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