Problem in internal energy calculation, rhoCentralFoam version > 2006
Summary
As told to @mark, we encountered in a problem related to calculation of internal energy e
for rhoCentralFoam in version > 2006
Steps to reproduce
Run one of the tutorials of rhoCentralFoam
, for instance, forwardStep and display the residuals (which are zero, direct solver) and the internal energy (that starts with negative values).
Example case
Load OpenFOAM V2112;
[spissoi@dvlogin02 forwardStep]$ module load openfoam/v2112/GccOpt
[spissoi@dvlogin02 forwardStep]$ module li
Currently Loaded Modulefiles:
1) lapack/gcc/64/3.9.0 2) blas/gcc/64/3.8.0 3) fftw3/openmpi/gcc/64/3.3.8 4) openmpi/mlnx/gcc/64/4.0.4rc3 5) gcc/10.2.0 6) boost/1.66.0/gnu/10.2 7) openfoam/v2112/GccOpt
cp -r solvers/compressible/rhoCentralFoam/rhoCentralFoam/ $FOAM_USER_APPBIN
Now modify the rhoCentralFoam.C
file in myrhoCentralFoam.C
to add a line to monitor initial value of e
at line 83
vi myrhoCentralFoam.C
82
83 Info<< "e\n" << e << endl;
84
85 Info<< "\nStarting time loop\n" << endl;
Adjust the relevant path in myrhoCentralFoam/ dir, renaming rhoCentralFoam.C
in myrhoCentralFoam.C
, and recompile:
./Allclean
./Allwmake
Now copy one of the tutorials, in that case forwardStep
cp $FOAM_TUTORIALS/compressible/rhoCentralFoam/forwardStep/ $FOAM_RUN
create the mesh with
blockMesh
And set only one-time step in system/controlDict
application rhoCentralFoam;
startFrom startTime;
startTime 0;
stopAt endTime;
endTime 0.002;
deltaT 0.002;
What is the current bug behaviour?
By running the instrumented solver
myrhoCentralFoam > log &
Once the simulation is done is it possible to find, in the log file, the info about the internal energy which is negative with a value of -743.589
If the standard set-up is run with standard rhoCentralFoam solver:
myrhoCentralFoam > log &
[spissoi@dvlogin02 forwardStep]$ vim rho.log
type hePsiThermo;
mixture pureMixture;
transport const;
thermo hConst;
equationOfState perfectGas;
specie specie;
energy sensibleInternalEnergy;
}
Reading field U
Creating turbulence model
Selecting turbulence model type laminar
Selecting laminar stress model Stokes
fluxScheme: Kurganov
Starting time loop
deltaT = 0.00238095
Mean and max Courant Numbers = 0.66643 0.666666
Time = 0.00238095
diagonal: Solving for rho, Initial residual = 0, Final residual = 0, No Iterations 0
diagonal: Solving for rhoUx, Initial residual = 0, Final residual = 0, No Iterations 0
diagonal: Solving for rhoUy, Initial residual = 0, Final residual = 0, No Iterations 0
diagonal: Solving for rhoE, Initial residual = 0, Final residual = 0, No Iterations 0
ExecutionTime = 0.11 s ClockTime = 0 s
deltaT = 0.000712548
Mean and max Courant Numbers = 0.199446 0.200515
Time = 0.0030935
diagonal: Solving for rho, Initial residual = 0, Final residual = 0, No Iterations 0
diagonal: Solving for rhoUx, Initial residual = 0, Final residual = 0, No Iterations 0
diagonal: Solving for rhoUy, Initial residual = 0, Final residual = 0, No Iterations 0
diagonal: Solving for rhoE, Initial residual = 0, Final residual = 0, No Iterations 0
ExecutionTime = 0.12 s ClockTime = 0 s
deltaT = 0.000712548
Mean and max Courant Numbers = 0.199447 0.201005
Time = 0.00380605
diagonal: Solving for rho, Initial residual = 0, Final residual = 0, No Iterations 0
diagonal: Solving for rhoUx, Initial residual = 0, Final residual = 0, No Iterations 0
diagonal: Solving for rhoUy, Initial residual = 0, Final residual = 0, No Iterations 0
diagonal: Solving for rhoE, Initial residual = 0, Final residual = 0, No Iterations 0
ExecutionTime = 0.13 s ClockTime = 0 s
Even if, by displaying the fields, it seems that the solution advance and produces the expected shock waves, see attached figure.
What is the expected correct behavior?
If openfoam v2006 is run, normal behavior is reproduced. By reproducing the same step as before, once the simulation is done it is possible to find, in the log file, the info about the internal energy which is positive with a value of 1.78572. And a normal residual behaviuor
Relevant logs and/or images
U field version openfoam v2112
Pressure field version openfoam v2112
Environment information
Providing details of your set-up can help us identify any issues, e.g. OpenFOAM version : v2112|v2006| Operating system : ubuntu laptop|centos HPC etc Hardware info : local installation on laptop and HPC installation Compiler : gcc
Possible fixes
Maybe it is due a non correct initialization of energy in the
constant/thermophysicalProperties
or in the creation of the initial field in
createFields.H
Thanks in advance
Ivan