Large time step continuity error in first time step with reactingFoam
Summary
Large time step continuity error in reactingFoam due to missing old time in thermo::psi field.
When running the tutorial tutorials/combustion/reactingFoam/laminar/counterFlowFlame2D, a large time step continuity error appears in the first iteration. This is independent if it starts from 0 or from a later time step.
The error could be traced back to the thermo.correct()
function called in the energy equation. There the values of the fields, especially the compressibility psi, are overwritten with new values based on the newly calculated temperature and known pressure. Even though GeometricField<>::primitiveFieldRef()
is called the old time field is not updated because no current old time pointer exists, see also void Foam::GeometricField<Type, PatchField, GeoMesh>::storeOldTimes()
.
This can be prevented by forcing an oldTime pointer to exist by calling in thermo::correct() this->psi_.oldTime()
, as had been done in e.g. OpenFOAM 5.x. Adding this line to the thermo::correct() function of hePsiThermo.C would reduce the observed error!
Steps to reproduce
Run tutorials/combustion/reactingFoam/laminar/counterFlowFlame2D and observe the time continuity error in the first time step. Then add the line
this->psi_.oldTime
to hePsiThermo::correct() and the error is reduced. It is also possible to check by printing out the number of old time values stored of psi before the energy equation, after the energy equation and after the pressure equation. It clearly shows, that only after it is solved for in the pressure equation with fvm::ddt(psi,p)
it has the old time fields.
Relevant logs and/or images
Example output of the first time step
Starting time loop
Courant Number mean: 2.73409e-05 max: 0.00175681
deltaT = 1.19999e-06
Time = 1.19999e-06
diagonal: Solving for rho, Initial residual = 0, Final residual = 0, No Iterations 0
PIMPLE: iteration 1
DILUPBiCGStab: Solving for O2, Initial residual = 1, Final residual = 1.02704e-09, No Iterations 1
DILUPBiCGStab: Solving for H2O, Initial residual = 0, Final residual = 0, No Iterations 0
DILUPBiCGStab: Solving for CH4, Initial residual = 1, Final residual = 1.0339e-09, No Iterations 1
DILUPBiCGStab: Solving for CO2, Initial residual = 0, Final residual = 0, No Iterations 0
DILUPBiCGStab: Solving for h, Initial residual = 1, Final residual = 4.10484e-09, No Iterations 1
min/max(T) = 293, 2000
DICPCG: Solving for p, Initial residual = 1, Final residual = 0.0488686, No Iterations 5
diagonal: Solving for rho, Initial residual = 0, Final residual = 0, No Iterations 0
time step continuity errors : sum local = 0.000333765, global = -0.000330463, cumulative = -0.000330463
DICPCG: Solving for p, Initial residual = 0.0147831, Final residual = 6.72619e-07, No Iterations 12
diagonal: Solving for rho, Initial residual = 0, Final residual = 0, No Iterations 0
time step continuity errors : sum local = 0.000331263, global = -0.000331263, cumulative = -0.000661725
ExecutionTime = 0.04 s ClockTime = 0 s
Environment information
- OpenFOAM version : v2012
- Operating system : Ubuntu 22.04
- Hardware info : AMD Ryzen 7 2700X Eight-Core Processor
- Compiler : GCC 9.5.0
Possible fixes
Change hePsiThermo.C
to:
template<class BasicPsiThermo, class MixtureType>
void Foam::hePsiThermo<BasicPsiThermo, MixtureType>::correct()
{
DebugInFunction << endl;
// force old time
this->psi_.oldTime();
calculate
(
this->p_,
this->T_,
this->he_,
this->psi_,
this->mu_,
this->alpha_,
false // No need to update old times
);
DebugInFunction << "Finished" << endl;
}