CFD Online Logo CFD Online URL
Home > Forums > OpenFOAM Running, Solving & CFD

conpressibleInterFoam for non-isothermal non-newtonian high-viscous flows

Register Blogs Members List Search Today's Posts Mark Forums Read

LinkBack Thread Tools Display Modes
Old   May 30, 2010, 21:37
Default conpressibleInterFoam for non-isothermal non-newtonian high-viscous flows
Jitao Liu
Join Date: Mar 2009
Location: Jinan , China
Posts: 64
Rep Power: 10
awacs is on a distinguished road
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:
Attached Images
File Type: jpeg Tait-equation.jpeg (17.0 KB, 159 views)
awacs is offline   Reply With Quote

Old   May 30, 2010, 21:44
Jitao Liu
Join Date: Mar 2009
Location: Jinan , China
Posts: 64
Rep Power: 10
awacs is on a distinguished road
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,
awacs is offline   Reply With Quote

Old   June 11, 2010, 14:30
Default Adding the temperature equation to compressibleInterFoam
New Member
Dan Gadensgaard
Join Date: Apr 2010
Posts: 13
Rep Power: 9
dgadensg is on a distinguished road
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!)
    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.
    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

dgadensg is offline   Reply With Quote

Old   June 23, 2010, 23:19
Jitao Liu
Join Date: Mar 2009
Location: Jinan , China
Posts: 64
Rep Power: 10
awacs is on a distinguished road
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.

awacs is offline   Reply With Quote

Old   January 24, 2011, 13:10
Senior Member
Illya Shevchuk
Join Date: Aug 2009
Location: Darmstadt, Germany
Posts: 176
Rep Power: 9
linch is on a distinguished road
Originally Posted by awacs View Post
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?

linch is offline   Reply With Quote

Old   March 18, 2011, 12:21
Senior Member
Ralph Moolenaar
Join Date: Aug 2010
Location: 's-Hertogenbosch, the Netherlands
Posts: 120
Rep Power: 8
Ralph M is on a distinguished road
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,

CFD for marine applications? Go to and join the OF Ship Hydromechanics Group:
Ralph M is offline   Reply With Quote

Old   January 8, 2012, 08:44
Senior Member
jeff osborne
Join Date: Mar 2010
Posts: 108
Rep Power: 9
ozzythewise is on a distinguished road
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
             + p*(fvm::div(U))
             -  fvm::laplacian(kT,T)
             - (tau && gradU)
I really have no idea, please can someone help?

ozzythewise is offline   Reply With Quote

Old   June 23, 2015, 05:22
Default Solving for enthalpy/internal energy instead of temperature in compressibleInterFoam
New Member
Quinn Reynolds
Join Date: Jun 2014
Posts: 6
Rep Power: 5
kittychunk is on a distinguished road
Originally Posted by ozzythewise View Post
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()
    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?

kittychunk is offline   Reply With Quote

Old   November 7, 2016, 02:52
Martin Aunskjaer
Join Date: Mar 2009
Location: Denmark
Posts: 52
Rep Power: 10
aunola is on a distinguished road
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 BCs 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.

Attached Files
File Type: gz cIEF.tar.gz (31.2 KB, 4 views)
aunola is offline   Reply With Quote


Thread Tools
Display Modes

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off
Trackbacks are On
Pingbacks are On
Refbacks are On

Similar Threads
Thread Thread Starter Forum Replies Last Post
y+ in high Pr flows vishyaroon Main CFD Forum 0 April 15, 2010 14:24
Oil (high Pr) flows: Near wall mesh resolution vishyaroon Main CFD Forum 0 April 13, 2010 09:37
Flows at High Angles of Attack Sarah C. FLUENT 2 March 18, 2009 08:39
High speed flows Raj FLUENT 4 January 10, 2006 05:14
Multicomponent fluid Andrea CFX 2 October 11, 2004 05:12

All times are GMT -4. The time now is 17:04.