[dvoosten]

Dear FOAM-ers,

I have been working with OpenFOAM for a few weeks now and I am coming to the point where I really have to check whether my results make any sense. The first thing I wanted to try was to see if the total energy in the system is conserved, so I modified the createFields.H file in rhoCentralFoam to write rhoE, so that I could externally integrate the total energy in the system. I checked that with the shockTube tutorial I indeed get a number that remains more or less constant. However, my own case it doesn't, in fact the energy rapidly doubles. My case is quite different from shockTube in the following way:

1) its use axi-symmetry using wedge boundary conditions
2) the initial temperature and pressure differences are quite a bit higher
3) the domain is not open on both ends (like in shocktube), but closed at one end

I don't understand the details of rhoCentralFoam well enough to judge how to control this problem. I tried varying the courant number, but this doesn't do much. Increasing the number of points does help, but only rather slowly. Am I doomed to keep increasing the number of point until the energy converges, or is there a shock-capturing solver that conserves energy a little better?

Thanks for your thoughts.

best,
Dries van Oosten

[dvoosten]

I've got some more details on a simple case that shows these problems.

Take the rhoCentralFoam/shockTube case and modify it such that the left side is a wall as opposed to a patch. Modify 0.orig/p, T and U accordingly. Run the case using Allrun. Afterwards, calculate the total mass and total energy for every timestep (using the included python scripts). You will find that they are constant to a high degree of accuracy.
time mass energy
0.001 1619.4241728 395803540.8
0.002 1619.4243168 395803509.12
0.003 1619.4243744 395803471.68
0.004 1619.4246048 395803465.92
0.005 1619.4246048 395803425.6
0.006 1619.42472 395803500.48
0.007 1619.4248928 395803396.8

Now, modify setFieldsDict such that the step pressure and temperature is at -4.5 (that is, closer to the left side, which we changed into a wall). Now run it and do the same analysis. We find
time mass energy
0.001 485.9774208 104401019.52
0.002 498.3053184 108470113.92
0.003 525.1805568 116935591.68
0.004 557.3084256 126861701.76
0.005 591.8136768 137447677.44
0.006 627.62688 148405492.8
0.007 664.2540576 159601962.24

A significant growth of both mass and energy!

Notice that this is at the pressure and temperature that is in the normal shockTube case. The only difference is the wall boundary. Moving the step closer to it makes the problem larger. Am I using the wrong boundary condition? I have a zeroGradient condition on that wall for p, T and U. Should I not be doing that?

As said, I have included the python code to calculate the total mass and energy. I have also include the Allrun script I use, which include a line to get the cell volume that the python code needs.