View Issue Details

Jump to NotesJump to History
IDProjectCategoryView StatusDate SubmittedLast Update
0003546OpenFOAMBugpublic2020-09-08 17:042020-09-11 12:28
Reporterjherb Assigned To 
PrioritynormalSeverityminorReproducibilityalways
Status newResolutionopen 
PlatformSUSE Linux Enterprise Server 11 OSOtherOS Version(please specify)
Summary0003546: Use of turbulentMixingLengthFrequencyInlet results in crash if velocity at patch is 0
DescriptionIf turbulentMixingLengthFrequencyBoundedInlet and turbulentIntensityKineticEnergyInlet are used for k and omega boundary conditions on a certain patch and the velocity at this patch is 0 (e. g. setting it by fixedValue to constant 0 over time or by using a table with time dependency with a certain period, where it is 0) then solution of the turbulence field equations will crash due to division by 0 on that patch.

The problem can be fixed by modifying the method turbulentMixingLengthFrequencyInletFvPatchScalarField::updateCoeffs() to be lower bounded by a non-zero value, e.g. by:
this->refValue() = max(sqrt(kp)/(Cmu25*mixingLength_), omegaMin_);
TagsNo tags attached.

Activities

tniemi

2020-09-09 06:43

reporter   ~0011499

The same issue also affects turbulentMixingLengthDissipationRateInletFvPatchScalarField, ie. k-epsilon combination. If U is zero in some face or at some time, k will be zero and hence epsilon leading to immediate crash.

will

2020-09-10 08:51

manager   ~0011500

There are a bunch of epsilon/k rates in the code, as well as the k^2/epsilon based calculation of the turbulent viscosity, so clipping epsilon and omega only really solves half the problem.

It would be better to change the turbulentIntensityKineticEnergyInlet condition to clip k given a minimum turbulent velocity fluctuation. That would fix all division issues, and a minimum/residual velocity magnitude is a more physically intuitive thing to ask users for than a value of epsilon or omega.

Better still would be to put a lower limit on the ratio of turbulent to laminar viscosity. That would be a nice non-dimensional setting so that it might be possible in some circumstances to set and forget it across multiple cases. The problem with doing this is that it couples the k and epsilon/omega boundary conditions more closely. To do it, the k condition would have to know the mixing length, which would make the BC less general and/or add a lot more logic.

Anyone got any better ideas as to how to apply this clip in as robust and user-friendly way as possible?

tniemi

2020-09-10 09:11

reporter   ~0011501

Residual value for k would be a simple and effective fix. Maybe it could have a default value of "small", as that is also used by default to bound the cell values of k, epsilon and omega. Could the bound-method be extended to cover also boundary values, as currently I think it only operates on cell values?

tniemi

2020-09-10 09:14

reporter   ~0011502

Hmm, my mistake, bound already does fix boundaryFieldRef. It has been a while since I encountered this issue, so I don't remeber where the division problem actually occured.

will

2020-09-10 09:56

manager   ~0011503

Last edited: 2020-09-10 09:56

 Yes, bound should sort it. Certainly after the model is corrected then the boundary fields should be clipped appropriately. It's possible that something in the evaluation sequence means that there are zeros hanging about for a bit causing this sort of error (like if something was doing a correctBoundaryConditions earlier on), but I can't identify what this would be, and I can't get a case to trigger it.

@jherb Can you give us a case that reproduces this problem?

tniemi

2020-09-10 10:12

reporter   ~0011504

I can also see if I can create a test case. When I had this issue, it was quite reproducible. Back then I "solved" the problem by simply forcing the inlet velocity to have a small positive value (it was a custom BC for velocity, with some faces possibly having pure 0). It did depend on the solver though, because if I remeber correctly, if I had p_rgh with fixedFluxPressure, the phi could be a tiny bit outwards, which turned the BCs into zeroGradient mode avoiding the issue.

jherb

2020-09-10 10:46

reporter   ~0011505

The problem can be reproduced in the $FOAM_TUT/heatTransfer/chtMultiRegionFoam/coolingSphere/ tutorial. Change the value of Uinlet in templates/0/fluid/U to (0 0 0) and the case will crash (tested with OpenFOAM-7):

Time = 0.0015


Solving for fluid region fluid
diagonal: Solving for rho, Initial residual = 0, Final residual = 0, No Iterations 0
DILUPBiCGStab: Solving for Ux, Initial residual = 0.17451096, Final residual = 9.4936014e-07, No Iterations 4
DILUPBiCGStab: Solving for Uy, Initial residual = 0.039102964, Final residual = 5.4949103e-07, No Iterations 4
DILUPBiCGStab: Solving for Uz, Initial residual = 0.35874008, Final residual = 7.5522014e-07, No Iterations 5
DILUPBiCGStab: Solving for h, Initial residual = 0.00090359875, Final residual = 9.720741e-07, No Iterations 1
Min/max T:295.99954 347.99059
GAMG: Solving for p_rgh, Initial residual = 0.52790977, Final residual = 0.0025859984, No Iterations 6
GAMG: Solving for p_rgh, Initial residual = 0.13235302, Final residual = 0.00048298206, No Iterations 4
diagonal: Solving for rho, Initial residual = 0, Final residual = 0, No Iterations 0
time step continuity errors : sum local = 5.1026755e-08, global = 2.5261954e-09, cumulative = 2.5261954e-09
Min/max rho:0.99883892 1.174286
GAMG: Solving for p_rgh, Initial residual = 0.10024547, Final residual = 0.00050102148, No Iterations 6
GAMG: Solving for p_rgh, Initial residual = 0.030951032, Final residual = 3.4537543e-08, No Iterations 13
diagonal: Solving for rho, Initial residual = 0, Final residual = 0, No Iterations 0
time step continuity errors : sum local = 5.9540449e-12, global = 2.4344848e-13, cumulative = 2.5264389e-09
Min/max rho:0.99883891 1.1742855
[1] #0 Foam::error::printStack(Foam::Ostream&) at ??:?
[1] #1 Foam::sigFpe::sigHandler(int) at ??:?
[1] #2 ? in "/lib64/libc.so.6"
[1] #3 Foam::divide(Foam::Field<double>&, Foam::UList<double> const&, Foam::UList<double> const&) at ??:?
[1] #4 Foam::tmp<Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh> > Foam::operator/<Foam::fvPatchField, Foam::volMesh>(Foam::tmp<Foam::GeometricField<double, Foam::fvPatchFi
eld, Foam::volMesh> > const&, Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh> const&) at ??:?
[1] #5 Foam::kOmegaSST<Foam::eddyViscosity<Foam::RASModel<Foam::EddyDiffusivity<Foam::ThermalDiffusivity<Foam::CompressibleTurbulenceModel<Foam::fluidThermo> > > > >, Foam::EddyDiffusivity
<Foam::ThermalDiffusivity<Foam::CompressibleTurbulenceModel<Foam::fluidThermo> > > >::correct() at ??:?
[1] #6 ? at ??:?
[1] #7 __libc_start_main in "/lib64/libc.so.6"
[1] #8 ? at /usr/src/packages/BUILD/glibc-2.11.3/csu/../sysdeps/x86_64/elf/start.S:116

Using a modified and bounded version of the turbulentMixingLengthFrequencyInlet boundary condition fixes the problem, i.e. the one line modification given in bug description above.

will

2020-09-10 11:24

manager   ~0011506

OK, thanks. The problem is that the BC isn't getting overridden by the bound operation because it is of mixed type and mixed BC-s aren't affected by assignment.

So, I reckon we're back at the idea of clipping the boundary condition itself.

It would be nice for the BC-s to clip using kMin/epsilonMin/omegaMin settings from the model; that would be logical and would save any additional input. I can't see how we'd access them, though. They are not part of the momentumTransport interface and nor should they be. We'd need another layer.

Lots of options, all of which imperfect.

tniemi

2020-09-11 10:47

reporter   ~0011507

I had time to investigate this a bit more and yes, the problem is that the bound does not work. In this case, as the BCs are derived from inletOutlet, they actually have a non-empty assignment operator, but it does not help because refValue() is zero and valueFraction() is one so the set value is ignored.

What if the k, epsilon, omega would have their own override of the operator=, which would just put the given value there?

tniemi

2020-09-11 11:13

reporter   ~0011508

I tested this minimal patch for just k and it lets the coolingSphere-tutorial to run.
patch.diff (1,909 bytes)    
diff --git a/src/finiteVolume/fields/fvPatchFields/derived/turbulentIntensityKineticEnergyInlet/turbulentIntensityKineticEnergyInletFvPatchScalarField.C b/src/finiteVolume/fields/fvPatchFields/derived/turbulentIntensityKineticEnergyInlet/turbulentIntensityKineticEnergyInletFvPatchScalarField.C
index 24ccbb5..48cde1f 100644
--- a/src/finiteVolume/fields/fvPatchFields/derived/turbulentIntensityKineticEnergyInlet/turbulentIntensityKineticEnergyInletFvPatchScalarField.C
+++ b/src/finiteVolume/fields/fvPatchFields/derived/turbulentIntensityKineticEnergyInlet/turbulentIntensityKineticEnergyInletFvPatchScalarField.C
@@ -155,6 +155,21 @@ void Foam::turbulentIntensityKineticEnergyInletFvPatchScalarField::write
 }
 
 
+// * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
+
+
+void Foam::turbulentIntensityKineticEnergyInletFvPatchScalarField::operator=
+(
+    const fvPatchScalarField& pvf
+)
+{
+    fvPatchScalarField::operator=
+    (
+        pvf
+    );
+}
+
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 namespace Foam
diff --git a/src/finiteVolume/fields/fvPatchFields/derived/turbulentIntensityKineticEnergyInlet/turbulentIntensityKineticEnergyInletFvPatchScalarField.H b/src/finiteVolume/fields/fvPatchFields/derived/turbulentIntensityKineticEnergyInlet/turbulentIntensityKineticEnergyInletFvPatchScalarField.H
index 61e8e33..72b182a 100644
--- a/src/finiteVolume/fields/fvPatchFields/derived/turbulentIntensityKineticEnergyInlet/turbulentIntensityKineticEnergyInletFvPatchScalarField.H
+++ b/src/finiteVolume/fields/fvPatchFields/derived/turbulentIntensityKineticEnergyInlet/turbulentIntensityKineticEnergyInletFvPatchScalarField.H
@@ -179,6 +179,11 @@ public:
 
         //- Write
         virtual void write(Ostream&) const;
+
+
+    // Member Operators
+
+        virtual void operator=(const fvPatchScalarField& pvf);
 };
patch.diff (1,909 bytes)   

will

2020-09-11 12:28

manager   ~0011509

Last edited: 2020-09-11 12:28

 Interesting. Thanks for the patch. I don't think we want to fix it that way, but it's helped me understand this some more.

None of this should be necessary. `phip` should be zero when `Up` is zero, so the `1 - pos0(phip)` thing that inletOutlet does for the value fraction should evaluate as zero, so the value should get overwritten on assignment. If you print out values for this case, though, `Up` is zero, but `phip` has some round off error in it. That's what's causing this.

This patch also fixes the problem. It makes sure the value fraction is zero when the velocity is zero. Again, though, ideally we wouldn't want to fix it this way. `phip` should be exactly zero here.

It's the buoyancy formulation. Something about it is accumulating error into `phip`. If you turn the gravity off, the problem also goes away.

bug3546_valueFractionFix.diff (3,967 bytes)    
diff --git a/src/MomentumTransportModels/momentumTransportModels/RAS/derivedFvPatchFields/turbulentMixingLengthDissipationRateInlet/turbulentMixingLengthDissipationRateInletFvPatchScalarField.C b/src/MomentumTransportModels/momentumTransportModels/RAS/derivedFvPatchFields/turbulentMixingLengthDissipationRateInlet/turbulentMixingLengthDissipationRateInletFvPatchScalarField.C
index f527511a1..38930b60d 100644
--- a/src/MomentumTransportModels/momentumTransportModels/RAS/derivedFvPatchFields/turbulentMixingLengthDissipationRateInlet/turbulentMixingLengthDissipationRateInletFvPatchScalarField.C
+++ b/src/MomentumTransportModels/momentumTransportModels/RAS/derivedFvPatchFields/turbulentMixingLengthDissipationRateInlet/turbulentMixingLengthDissipationRateInletFvPatchScalarField.C
@@ -144,13 +144,16 @@ void turbulentMixingLengthDissipationRateInletFvPatchScalarField::updateCoeffs()
     const fvPatchScalarField& kp =
         patch().lookupPatchField<volScalarField, scalar>(kName_);
 
+    const fvPatchVectorField& Up =
+        patch().lookupPatchField<volVectorField, vector>(/*UName_*/ "U");
+
     const fvsPatchScalarField& phip =
         patch().lookupPatchField<surfaceScalarField, scalar>(this->phiName_);
 
     this->refValue() = Cmu75*kp*sqrt(kp)/mixingLength_;
-    this->valueFraction() = 1.0 - pos0(phip);
+    this->valueFraction() = (1.0 - pos0(phip))*pos(magSqr(Up));
 
-    inletOutletFvPatchScalarField::updateCoeffs();
+    mixedFvPatchField<scalar>::updateCoeffs();
 }
 
 
diff --git a/src/MomentumTransportModels/momentumTransportModels/RAS/derivedFvPatchFields/turbulentMixingLengthFrequencyInlet/turbulentMixingLengthFrequencyInletFvPatchScalarField.C b/src/MomentumTransportModels/momentumTransportModels/RAS/derivedFvPatchFields/turbulentMixingLengthFrequencyInlet/turbulentMixingLengthFrequencyInletFvPatchScalarField.C
index 5fcfa17d3..d46ce02a3 100644
--- a/src/MomentumTransportModels/momentumTransportModels/RAS/derivedFvPatchFields/turbulentMixingLengthFrequencyInlet/turbulentMixingLengthFrequencyInletFvPatchScalarField.C
+++ b/src/MomentumTransportModels/momentumTransportModels/RAS/derivedFvPatchFields/turbulentMixingLengthFrequencyInlet/turbulentMixingLengthFrequencyInletFvPatchScalarField.C
@@ -140,13 +140,16 @@ void turbulentMixingLengthFrequencyInletFvPatchScalarField::updateCoeffs()
     const fvPatchScalarField& kp =
         patch().lookupPatchField<volScalarField, scalar>(kName_);
 
+    const fvPatchVectorField& Up =
+        patch().lookupPatchField<volVectorField, vector>(/*UName_*/ "U");
+
     const fvsPatchScalarField& phip =
         patch().lookupPatchField<surfaceScalarField, scalar>(this->phiName_);
 
     this->refValue() = sqrt(kp)/(Cmu25*mixingLength_);
-    this->valueFraction() = 1.0 - pos0(phip);
+    this->valueFraction() = (1.0 - pos0(phip))*pos(magSqr(Up));
 
-    inletOutletFvPatchScalarField::updateCoeffs();
+    mixedFvPatchField<scalar>::updateCoeffs();
 }
 
 
diff --git a/src/finiteVolume/fields/fvPatchFields/derived/turbulentIntensityKineticEnergyInlet/turbulentIntensityKineticEnergyInletFvPatchScalarField.C b/src/finiteVolume/fields/fvPatchFields/derived/turbulentIntensityKineticEnergyInlet/turbulentIntensityKineticEnergyInletFvPatchScalarField.C
index 24ccbb513..a43f85a23 100644
--- a/src/finiteVolume/fields/fvPatchFields/derived/turbulentIntensityKineticEnergyInlet/turbulentIntensityKineticEnergyInletFvPatchScalarField.C
+++ b/src/finiteVolume/fields/fvPatchFields/derived/turbulentIntensityKineticEnergyInlet/turbulentIntensityKineticEnergyInletFvPatchScalarField.C
@@ -136,9 +136,9 @@ updateCoeffs()
         patch().lookupPatchField<surfaceScalarField, scalar>(this->phiName_);
 
     this->refValue() = 1.5*sqr(intensity_)*magSqr(Up);
-    this->valueFraction() = 1.0 - pos0(phip);
+    this->valueFraction() = (1.0 - pos0(phip))*pos(magSqr(Up));
 
-    inletOutletFvPatchScalarField::updateCoeffs();
+    mixedFvPatchField<scalar>::updateCoeffs();
 }
bug3546_valueFractionFix.diff (3,967 bytes)   

Issue History

Date Modified Username Field Change
2020-09-08 17:04 jherb New Issue
2020-09-09 06:43 tniemi Note Added: 0011499
2020-09-10 08:51 will Note Added: 0011500
2020-09-10 09:11 tniemi Note Added: 0011501
2020-09-10 09:14 tniemi Note Added: 0011502
2020-09-10 09:56 will Note Added: 0011503
2020-09-10 09:56 will Note Edited: 0011503
2020-09-10 10:12 tniemi Note Added: 0011504
2020-09-10 10:46 jherb Note Added: 0011505
2020-09-10 11:24 will Note Added: 0011506
2020-09-11 10:47 tniemi Note Added: 0011507
2020-09-11 11:13 tniemi File Added: patch.diff
2020-09-11 11:13 tniemi Note Added: 0011508
2020-09-11 12:28 will File Added: bug3546_valueFractionFix.diff
2020-09-11 12:28 will Note Added: 0011509
2020-09-11 12:28 will Note Edited: 0011509