Turbulence models (kEpsilon, kOmegaSSTSato) TwoPhaseSystem - preconditioner (#1789) · Issues · Development / openfoam

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; }

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; }

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)

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);

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);

alpha1.max(SMALL); alpha1.min(1 - SMALL);

<!--
  *** Please read this first! ***

Before opening a new issue, make sure to search for keywords in the issues
  filtered by the "bug" label and check to see if it has already been reported

You can see how your report will be rendered on the platform by using the
  "preview" tab above
-->

### 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;                                                                                                                              
}`

`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;                                                                                                                              
}`

`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

```#0  Foam::error::printStack(Foam::Ostream&) at ~/OpenFOAM/OpenFOAM-v1906-debug/src/OSspecific/POSIX/printStack/printStack.C:236
#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.
<!--
  If possible, please create a SMALL example and attach it to your report

If you are using an older version of OpenFOAM this will also determine
  whether the bug has been fixed in a more recent version
-->

### 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);`

`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);
`

`
alpha1.max(SMALL);
alpha1.min(1 - SMALL);
`

<!--
  If you can, link to the line of code that might be responsible for the
  problem

The "/label ~bug" text is a gitlab flag that will add the "bug" label to this
  issue
-->

Edited Jul 26, 2020 by Robin Kamenicky

To upload designs, you'll need to enable LFS. More information