CFD Online Discussion Forums

CFD Online Discussion Forums (
-   OpenFOAM Running, Solving & CFD (
-   -   OpenFoam vs CFX5 mass balance in OpenFoam (

tangd August 30, 2006 06:46

Hello, I ran my case with b

I ran my case with both OpenFoam (sonicTurbFoam) and CFX5 and compared
results by generating graphs of temperature, pressure, mass and mass flow rate
as well as some contour screenshots (see attached images). In principle,
the axisymmetric-2D-system is a simple volume with an inlet. For OpenFoam
a 5° sector was used with wedges.

I found that temperature and pressure of OpenFoam are lower than the corresponding
results of CFX-5. Additionally, in OpenFoam the overall mass in the closed volume
is lower than the integral value of the inlet mass flow rate. Thus, it looks as if
mass is lost during the calculation. Note, the massflow input (and also the
total temperature) at inlet is nearly the same as in CFX5. The mesh was the same,
of course.

I have checked a lot of things: turbulence model, boundary conditions, thermal
properties, etc, but haven't found the reason why in OpenForm the mass balance (and
consequently the energy balance) is not fulfilled. The schemes I used are:

================================================== ============================
default CrankNicholson 0.5;

default Gauss linear;

default none;
div(phi,U) Gauss limitedLinearV 1;
div(phi,k) Gauss upwind;
div(phi,epsilon) Gauss upwind;
div(phi,R) Gauss upwind;
div(R) Gauss linear;
div(phid,p) Gauss limitedLinear 1;
div(phiU,p) Gauss limitedLinear 1;
div(phi,h) Gauss limitedLinear 1;
div((muEff*dev2(grad(U).T()))) Gauss linear;

default none;
laplacian(muEff,U) Gauss linear corrected;
laplacian(DkEff,k) Gauss linear corrected;
laplacian(DREff,R) Gauss linear corrected;
laplacian(DepsilonEff,epsilon) Gauss linear corrected;
laplacian((rho*1|A(U)),p) Gauss linear corrected;
laplacian(alphaEff,h) Gauss linear corrected;

default linear;

default corrected;

default no;
================================================== ============================

I've also used Euler time discretisation, but nothing changed. So far, I didn't
change the spatial discretisation, it's still the default setting from prism.
The solution settings are:

================================================== ============================
p BICCG 1e-12 0;
rho ICCG 1e-12 0;
U BICCG 1e-12 0;
h BICCG 1e-12 0;
k BICCG 1e-12 0;
epsilon BICCG 1e-12 0;
R BICCG 1e-12 0;

momentumPredictor yes;
nCorrectors 2;
nNonOrthogonalCorrectors 0;
pRefCell 0;
pRefValue 0;
================================================== ============================

Thus, does anybody have an idea where the problem could be? The unconserevd
energy and mass balance make me uncomfortable ;-(

lr103476 August 30, 2006 07:13

Don't you need CrankNicholson
Don't you need CrankNicholson 1.0 for fully second order accuracy?

Regards, Frank

tangd August 30, 2006 07:26

I also used Euler here, it doe
I also used Euler here, it doesn't matter so much I think. Other ideas?

david_h August 30, 2006 08:11

How did you calculate the syst
How did you calculate the system mass, m, the dotted curve in the t vs. (m0 - m) plot ?

tangd August 30, 2006 08:24

The system mass is calculated
The system mass is calculated by using pV=mRT, p and T are the mean values within the volume. Then I get the mass(m) value for each saved step, subtract the initial mass m0.

david_h August 30, 2006 09:30

If I understood correctly you
If I understood correctly you are calculating the mass as m = V/R {p}/{T}. Where { } denotes the volume average. Have you tried calculating the mass as m = V/R {p/T} ?

tangd August 31, 2006 03:41

Yes, your idea turned out to b
Yes, your idea turned out to be correct! I have checked the (m-m0) again by using your proposal m = V/R {p/T}, the mass is now conserved. Thanks a lot! But any idea about the pressure and Temperature?

olwi August 31, 2006 03:42

If you really want to integrat
If you really want to integrate the influx, you should look at *exactly* what goes in. That is exactly what's on the boundary patch. Not what's in the cell next to it... In addition, your own way of integarting plays a role. Is *your* integration 100% identical with the OpenFOAM time integration? For example using mid-point values matching the CN scheme?

It's relatively easy to assure mass conservation in a fv code like OpenFOAM, and from the quality in detail I've seen so far in OpenFOAM it would be profoundly surprising if you found a lack of mass conservation. "Surprising" is too weak a word, actually.


tangd August 31, 2006 04:04

You are right! My incorrect po
You are right! My incorrect postprocessing calculation introduced the error. Now mass is conserved. But there are still problems during the simulation that make temperature and pressure smaller than that of CFX-5. any further suggestion?

olwi August 31, 2006 04:23

Then remains to compare exactl
Then remains to compare exactly how much influx of enthalpy you have in the two codes. If I understood correclty, there's no heat production inside? If the mass inventory is the same in both codes at time t, then you only need to look for the enthalpy influx. Same here: I doubt they've missed conservation in CFX5 or OpenFOAM... It's likley a matter of how sources or inlet properties are defined. Even if you think you do they same in both codes, the codes may have different conventions for how things are defined on the boundary.

tangd September 1, 2006 08:42

The finding is that the enthal
The finding is that the enthalpy flow of OpenFOAM is very close to that of CFX5.

Here, we have checked the static part, mdot*h_static, as well as the kinetic part mdot*uČ/2. Both are the same on the inlet for OpenFoam and CFX.

I inspected this by considering each cell at inlet then summing up them. Using the same method, I also checked the pressures and velocities at inlet from both tools. They are nearly the same. The inlet calculation seems to be no problem. But I still cannot figure out what kind of problem causes the different pressure and temperature values within the volume. Can you give me more suggestions?

david_h September 1, 2006 08:54

Did you run CFX with viscous d
Did you run CFX with viscous dissipation. The energy equation in sonicTurbFoam, "hEqn.H", does not appear to have a viscous dissipation term, which might explain a few degrees of temperature difference.

tangd September 1, 2006 09:14

In the setting of heat transfe
In the setting of heat transfer model "Total Energy" in CFX5, option "Incl. viscous Work Term" was _not_ activated.

olwi September 2, 2006 09:08

Have you done the same check o
Have you done the same check on energy/enthalpy conservation as you did for the mass? (total energy in the flow domain before/after in relation to the enthalpy influx.) If energy is conserved, then I guess there's only the consitutive relationship, etc., relating pressure to temperature and enthalpy...

tangd September 4, 2006 04:00

Now the situation is that if I
Now the situation is that if I subtract the mean temperature <t> of OpenFoam with 293 and multiply with gamma (gamma=cp/cv), then I get almost the same temperature as in CFX5. So I guess the wrong equation is being used in calculating temperature in OpenFOAM. Any suggestion?

olwi September 4, 2006 06:04

It's a bit hard to say OpenFOA
It's a bit hard to say OpenFOAM solves the "wrong" equation, isn't it? Sounds to me OF uses SI units for temperature (Kelvin), while CFX5 uses something else. The factor gamma can be explained if someone used cv instead of cp, or vice versa. Sounds to me it's more a question of reading in the manuals/case files to know what is actually calculated and ensuring you use the same physical properties and units in both codes. That's obviously the user's responsability. Should be easy enough in OF, just have a look in the source code.

tangd September 5, 2006 06:18

david_h September 5, 2006 13:31

One other thing you might take
One other thing you might take look at is the turbulence model. For example in:
the R() function includes a "(2/3)*k" term, however divRhoR() does not include this term.

The quickest fix is to add
-(2/3)*fvc::grad(rho*turbulence->k() )
after the fvc::grad(p) term when the momentum equation is solved.


tangd September 6, 2006 06:12

Thank you so much. I have chan
Thank you so much. I have changed the code by using your suggestion. But unfortunately, the result is nearly the same as before.

I'm still thinking the "gamma" plays the key role here. So please give me more hints based on my previous post. Thanks again!

chris1980 September 6, 2006 08:05

Can someone comment if David i
Can someone comment if David is right because this would be a potential bug!

All times are GMT -4. The time now is 09:01.