Turbulence models (kEpsilon, kOmegaSSTSato) TwoPhaseSystem - preconditioner
Summary
When kEpsolin, kOmegaSST, kOmegaSSTSato is used and the alpha fields are initialized via utility such as setField to value 0 or 1 the preconditioner gives floating point exception during usage of preconditioner class. This happens during first iteration.
Steps to reproduce
Use the tutorial $FOAM_TUTORIALS/heatTransfer/chtMultiRegionTwoPhaseEulerFoam/solidQuenching2D. Change the constant/water/turbulenceProperties.liquid from
simulationType RAS; RAS { RASModel mixtureKEpsilon;//kEpsilon; turbulence on; printCoeffs on; }
to
simulationType RAS; RAS { RASModel kEpsilon; turbulence on; printCoeffs on; }
and change file constant/water/turbulenceProperties.gas from
simulationType RAS; RAS { RASModel mixtureKEpsilon; turbulence on; printCoeffs on; }
to
simulationType RAS; RAS { RASModel kEpsilon; turbulence on; printCoeffs on; }
The same happens if we use for instance a combination for gas-liquid kEpsilon-laminar or kOmegaSST-kOmegaSSTSato.
The error is
#1 Foam::sigFpe::sigHandler(int) at ~/OpenFOAM/OpenFOAM-v1906-debug/src/OSspecific/POSIX/signals/sigFpe.C:130
#2 ? in /lib/x86_64-linux-gnu/libc.so.6
#3 Foam::DILUPreconditioner::calcReciprocalD(Foam::Field<double>&, Foam::lduMatrix const&) at ~/OpenFOAM/OpenFOAM-v1906-debug/src/OpenFOAM/matrices/lduMatrix/preconditioners/DILUPreconditioner/DILUPreconditioner.C:80 (discriminator 2)
It happens already at first construction of the class. It seems that division by 0 occurs on line 80 in DILUPreconditioner.C. Change of preconditioner or solver does not help. The only option is not to use preconditioner at all. When that is done, the simulation runs smoothly and the results seem to be fine.
Example case
Described in the steps to reproduce.
What is the current bug behaviour?
SIGFPE when preconditioner used with kEpsilon, kOmegaSST, kOmegaSSSato and internalField of phase volume fractions initialized to 0 or 1.
What is the expected correct behavior?
No SIGFPE.
Relevant logs and/or images
Environment information
- OpenFOAM version : v1906
- Operating system : Ubuntu 16.04.6 LTS
- Hardware info :
- Compiler : gcc
Possible fixes
Bound the phase fraction in $FOAM_SRC/phaseSystemModels/reactingEulerFoam/reactingTwoPhaseEulerFoam/twoPhaseSystem/twoPhaseSystem.C:342
from
alpha1.clip(0, 1);
to
alpha1.clip(SMALL, 1 - SMALL);
the same will probably happen also for twoPhaseEulerFoam. I have not tested that. the fix would be in $FOAM_APP/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/twoPhaseSystem.C
where lines 532, 533 would need to be changed from
alpha1.max(0); alpha1.min(1);
to
alpha1.max(SMALL); alpha1.min(1 - SMALL);