[awacs]

Dear foamers,

I am trying to simulate two compressible, non-isothermal, non-newtonian (high viscous), laminar, and immiscible fluids flow. There is already a compressibleInterFoam solver for 2 compressible, isothermal immiscible fluids using a VOF (volume of fluid) phase-fraction based interface capturing approach. I want to create a new compressibleNonIsothermalInterFoam. Tait state equation is used to express the fuid dencity varing with temperature and pressure:

[awacs]

Is it possible to modify the compressibleInterFoam solver to take the temperature into account for non-isothermal compressible laminar flows? Any suggestions will be appreciated.

Best regards,
Jitao

[dgadensg]

Hi Jitao.

I have succesfully implemented the temperature equation in the interFoam solver, and the procedure should be somewhat the same for the compressible case.

1. First you have to add the appropriate constants and fields to the createFields.H file. This is best done by looking at step 3 in the guide provided here.
   
2. Next you should add the temperature equation to the compressibleInterFoam solver algorithm. I did this by making a new file "TEqn.H" and then included this just after the PISO loop. The contents of the "TEqn.H" i added is shown below:
   
   

   volScalarField mu //dynamic viscosity (at cell centers!)
       (
           "mu",
           twoPhaseProperties.mu()
       );
   
   volTensorField gradU = fvc::grad(U); //velocity gradient.
   
   volTensorField tau = mu * (gradU + gradU.T()); // viscous stress tensor
   
   //Temperature equation (with viscous dissipation):
   fvScalarMatrix TEqn
   	(
   		
   		rho * cp * (fvm::ddt(T) + fvm::div(phi,T))
   		- fvm::laplacian(kT,T) 
   		- ( tau && gradU ) //viscous dissipation heating.
   	);
   
   TEqn.solve();

   You should note that the above temperature equation is for an incompressible flow with constant thermal conductivity and specifik heat. Furthermore the equation takes the viscous dissipation heating into account.
   
   You should use a compressible version of this equation, which can be found in the litterature (e.g. Transport Phenomena by R. Byron Bird et al. - ISBN: 978-0-470-11539-8). Also note that the stress tensor should be corrected to the compressible case as well, but you should consider whether it is necessary for you to add the viscous dissipation term in your application.
   
3. When this is done you should be able to compile the new solver and have a nonisothermal version of the compressibleInterFoam. Before running a case, be sure to add the appropriate boundary conditions for the temperature, as well as to choose the schemes for the new terms in "fvSchemes" and a solver algorithm in "fvSolution".

Let me know if there is any further questions! Have you implemented the Tait state equation yet? If you have, i would like to know/see how you have done this?

With Kind Regards

Dan

[awacs]

Hi Dan,

I have added temperature field into compressibleInterFoam. The new solver works well in two-phase (polymer and air) flow simulation. Constant thermal conductivity and specifik heat are assumed in this silver.

The compressiblility is implemented in pEqn.H:
rho1 = rho10 + psi1*p;
rho2 = rho20 + psi2*p;
I am trying to take into account the effects of temperature on dencity. The aforementioned Tait state equation is such a model defining dencity as a function of temperature and pressure.

Sincerely,
Jitao

[linch]

Originally Posted by awacs

I am trying to take into account the effects of temperature on dencity. The aforementioned Tait state equation is such a model defining dencity as a function of temperature and pressure.

Hi Jitao,

have you succeeded with the implementation of the new equation of state with temperature and pressure depending density?

regards,
Ilya

[Ralph M]

Hello Dan and Jitao (and others),

I have some practical questions concerning the implementation of the heat equation.

First a remark: the laplacian in the code of Dan should be fvm::laplacian(DT,T) according to step 1..right?

I'm a but confused how DT and cp are calculated since both fluids have their own properties. Right now I'm adding read-terms in the "twoPhaseMixture"-files (both the .H and .C) to read these terms but I'm not sure whether this would be correct. these paramaters are then defined in the transportProperties-file?

Thanks for your support,

Ralph

[ozzythewise]

I realize that this thread is a little old, but I am in some desperate need of help. I'm trying to do the same thing that awacs was doing, that is implementing non-isothermal into the compressibleInterFoam solver, but am quite confused as to what to edit in the TEqn.H file. Isn't it when you account for compressibility you need to use enthalpy instead of temperature? That seems to be how all the built-in compressible solvers do for the rest of OF. I'm still trying to do it with TEqn.H, but I'm not sure how to edit it, would it be something like:

fvScalarMatrix TEqn
	(
		cp*(fvm::ddt(rho,T))
             + p*(fvm::div(U))
             -  fvm::laplacian(kT,T)
             - (tau && gradU)
	);

I really have no idea, please can someone help?

Thanks
Jeff

[kittychunk]

Originally Posted by ozzythewise

Isn't it when you account for compressibility you need to use enthalpy instead of temperature? That seems to be how all the built-in compressible solvers do for the rest of OF.

I'm having the same issue. Unless I'm missing something, the compressible multiphase flow solvers shipping with OF as of version 2.4.0 seem to assume that internal energy e = Cv*T (rather than the usual integral over T), implying that Cv is constant with temperature. In general this is not the case, and a formulation of the energy conservation equation should rather use enthalpy or internal energy directly as the raw field variable; this is the way it's done in OF's vanilla compressible flow solvers. It looks like a few others have noticed this, but I couldn't find any useful answers on to how to address the issue

It's not even possible to create a reference to (for example) the mixture.he() field and then write your own "EEqn.H" for compressibleInterFoam - the solver will happily compile if you do this, but when you run a case you'll find that it will throw a "Not Implemented" error. This is because for some reason he as a field reference is not implemented in the twoPhaseMixtureThermo and multiPhaseMixtureThermo libraries. For example in "multiphaseMixtureThermo.H":

//- Enthalpy/Internal energy [J/kg]
//  Non-const access allowed for transport equations
virtual volScalarField& he()
{
    notImplemented("multiphaseMixtureThermo::he()");
    return phases_[0]->thermo().he();
}

//- Enthalpy/Internal energy [J/kg]
virtual const volScalarField& he() const
{
    notImplemented("multiphaseMixtureThermo::he() const");
    return phases_[0]->thermo().he();
}

Short of writing one's own phase mixing thermo library, I'm wondering if is there any clever way to address this?

-q

[aunola]

The temperature equation in compressibleInterFoam (cIF) is derived under the assumption that the specific internal energy of each phase depends only on temperature and is given by e = c_v*T, where c_v is the constant isochoric specific heat capacity. Formally, this holds only for calorically perfect gases. It is not possible to use cIF in its current form for real gases or even thermally perfect gases (de = c_v(T)dT). In practice, this is enforced through the thermo library selection restricting what is available for the sensibleInternalEnergy energy variable required by cIF.

The use of enthalpy or energy as a primary unknown instead of temperature reflects its role as the conserved variable in the energy equation. This choice of primary unknown introduces complexity because many of the thermodynamic variables are given in terms of temperature. In particular, the functional form of the enthalpy/energy is he = he(p,T), which implicitly defines temperature as a function of pressure and enthalpy/energy. Temperature does not occur naturally in the conservation equations. For a compressible two-phase VOF-solver like cIF the defining equation for temperature is then rho*he(p,T) = alpha1*rho1*he1(p,T) + alpha2*rho2*he2(p,T). Temperature must be found by an iterative approach to solving this equation, given that the flow solver produces all the values except the he1/2.

Extending cIF to solve for the conserved variable in the energy equation (total enthalpy or energy) instead of temperature requires sorting out both implementation and numerical problems.

With a view towards the current thermo library implementation in OF at least the following implementation issues arise:

• A mixture enthalpy/energy equation is the sum of enthalpy/energy equations for each phase. This sum cannot be formed unless differences in the enthalpy/energy reference states of the fluids are accounted for. With the current assumption made by cIF as to the form of the phasic internal energy this is implicitly enforced, but in general it is not.
• As there is one thermo object per phase, there is no overall mixture knowledge in a thermo object. Consequently, either a new two-phase thermo structure must be created as a parallel to the existing single-phase thermo classes or it must be implemented outside of the current thermo structure.

The following numerical issues can be anticipated:

In a segregated solver with no coupling between the energy and volume fraction advection equations there is no reason to expect that the enthalpy/energy field is in exact agreement with the phase volume fraction and density fields. Consequently, iterating temperature from these uncoupled fields is likely to produce unphysical values in cells/faces containing the interface. This will be exacerbated when the enthalpy/energy levels of the phases are very different. Consider as an example a cell containing the interface between dry air and liquid water. The absolute specific enthalpy of liquid water is about a factor of four higher than that of dry air at room temperature (both referred to the triple point of water). Thus, if the enthalpy in the cell as calculated by the energy equation is not in exact agreement with the liquid volume fraction, the contribution from water to the cell enthalpy will either be significantly too low or too high. Solving for the temperature will not produce the correct result. It is also conceivable that different numerical schemes will have different effects on the inconsistency between the density, volume fraction and enthalpy/energy fields. Hence, in general, unphysical temperature fluctuations are to be expected in interface cells when temperature is derived from the conserved energy variable in a VOF-solver.

The energy BC’s should still work, I think, but it should at least be considered if a revision is necessary.

To demonstrate the above problems I made a quick modification to cIF which replaces the current temperature equation with an enthalpy equation and implements its own thermo library . I simply undocked all the necessary thermo functions from the OF thermo tree, placed them in the twoPhaseMixtureThermo.H/C, adapted as needed and created a number of pseudo-classes to allow objects of this new thermo class to be recognized as a thermo library. Running this solver on a test case (toy water rocket in 2D) immediately shows the problem of unphysical temperatures in interface cells. Eventually, other fields will become unphysical. The same case (with appropriate modifications to the thermophysicalProperties dictionary) runs fine with cIF.

Experimental solver and case attached. It is for OF-2.2.x but easily portable to later versions.

I stress that the solver is only meant to demonstrate the numerical problems – it does not address them.

How to address these issues? This is an interesting problem. Perhaps a coupled solver is needed to ensure consistency between the fields.

[chaohui]

[QUOTE=aunola;624441]The temperature equation in compressibleInterFoam (cIF) is derived under the assumption that the specific internal energy of each phase depends only on temperature and is given by e = c_v*T, where c_v is the constant isochoric specific heat capacity. Formally, this holds only for calorically perfect gases. It is not possible to use cIF in its current form for real gases or even thermally perfect gases (de = c_v(T)dT). In practice, this is enforced through the thermo library selection restricting what is available for the sensibleInternalEnergy energy variable required by cIF.

The use of enthalpy or energy as a primary unknown instead of temperature reflects its role as the conserved variable in the energy equation. This choice of primary unknown introduces complexity because many of the thermodynamic variables are given in terms of temperature. In particular, the functional form of the enthalpy/energy is he = he(p,T), which implicitly defines temperature as a function of pressure and enthalpy/energy. Temperature does not occur naturally in the conservation equations. For a compressible two-phase VOF-solver like cIF the defining equation for temperature is then rho*he(p,T) = alpha1*rho1*he1(p,T) + alpha2*rho2*he2(p,T). Temperature must be found by an iterative approach to solving this equation, given that the flow solver produces all the values except the he1/2.

Extending cIF to solve for the conserved variable in the energy equation (total enthalpy or energy) instead of temperature requires sorting out both implementation and numerical problems.

With a view towards the current thermo library implementation in OF at least the following implementation issues arise:

• A mixture enthalpy/energy equation is the sum of enthalpy/energy equations for each phase. This sum cannot be formed unless differences in the enthalpy/energy reference states of the fluids are accounted for. With the current assumption made by cIF as to the form of the phasic internal energy this is implicitly enforced, but in general it is not.
• As there is one thermo object per phase, there is no overall mixture knowledge in a thermo object. Consequently, either a new two-phase thermo structure must be created as a parallel to the existing single-phase thermo classes or it must be implemented outside of the current thermo structure.

The following numerical issues can be anticipated:

In a segregated solver with no coupling between the energy and volume fraction advection equations there is no reason to expect that the enthalpy/energy field is in exact agreement with the phase volume fraction and density fields. Consequently, iterating temperature from these uncoupled fields is likely to produce unphysical values in cells/faces containing the interface. This will be exacerbated when the enthalpy/energy levels of the phases are very different. Consider as an example a cell containing the interface between dry air and liquid water. The absolute specific enthalpy of liquid water is about a factor of four higher than that of dry air at room temperature (both referred to the triple point of water). Thus, if the enthalpy in the cell as calculated by the energy equation is not in exact agreement with the liquid volume fraction, the contribution from water to the cell enthalpy will either be significantly too low or too high. Solving for the temperature will not produce the correct result. It is also conceivable that different numerical schemes will have different effects on the inconsistency between the density, volume fraction and enthalpy/energy fields. Hence, in general, unphysical temperature fluctuations are to be expected in interface cells when temperature is derived from the conserved energy variable in a VOF-solver.

The energy BC’s should still work, I think, but it should at least be considered if a revision is necessary.

To demonstrate the above problems I made a quick modification to cIF which replaces the current temperature equation with an enthalpy equation and implements its own thermo library . I simply undocked all the necessary thermo functions from the OF thermo tree, placed them in the twoPhaseMixtureThermo.H/C, adapted as needed and created a number of pseudo-classes to allow objects of this new thermo class to be recognized as a thermo library. Running this solver on a test case (toy water rocket in 2D) immediately shows the problem of unphysical temperatures in interface cells. Eventually, other fields will become unphysical. The same case (with appropriate modifications to the thermophysicalProperties dictionary) runs fine with cIF.

Experimental solver and case attached. It is for OF-2.2.x but easily portable to later versions.

I stress that the solver is only meant to demonstrate the numerical problems – it does not address them.

How to address these issues? This is an interesting problem. Perhaps a coupled solver is needed to ensure consistency between the fields.

[/QUOTE
I have tried to install in of version 2.2, but the compiler is still a problem