BUG: twoPhaseEulerFoam: bug fixes and bounding improvements to interfacal models

applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/aspectRatioModels/TomiyamaAspectRatio/TomiyamaAspectRatio.C

@@ -52,7 +52,7 @@ Foam::aspectRatioModels::TomiyamaAspectRatio::TomiyamaAspectRatio
     const orderedPhasePair& pair
 )
 :
-    aspectRatioModel(dict, pair),
+    VakhrushevEfremov(dict, pair),
     yWall_(pair.phase1().mesh().lookupObject<volScalarField>("yWall"))
 {}
 
@@ -69,7 +69,7 @@ Foam::tmp<Foam::volScalarField>
 Foam::aspectRatioModels::TomiyamaAspectRatio::E() const
 {
     return
-        pair_.Eo()
+        VakhrushevEfremov::E()
        *max
        (
            scalar(1) - 0.35*yWall_/pair_.dispersed().d(),

applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/aspectRatioModels/TomiyamaAspectRatio/TomiyamaAspectRatio.H

@@ -44,7 +44,7 @@ SourceFiles
 #ifndef TomiyamaAspectRatio_H
 #define TomiyamaAspectRatio_H
 
-#include "aspectRatioModel.H"
+#include "VakhrushevEfremov.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -59,7 +59,7 @@ namespace aspectRatioModels
 
 class TomiyamaAspectRatio
 :
-    public aspectRatioModel
+    public VakhrushevEfremov
 {
 private:
 

applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/dragModels/TomiyamaAnalytic/TomiyamaAnalytic.C

@@ -69,17 +69,17 @@ Foam::dragModels::TomiyamaAnalytic::CdRe() const
     volScalarField Eo(max(pair_.Eo(), residualEo_));
     volScalarField E(max(pair_.E(), residualE_));
 
-    volScalarField OmEsq(max(scalar(1) - sqr(E), residualE_));
+    volScalarField OmEsq(max(scalar(1) - sqr(E), sqr(residualE_)));
     volScalarField rtOmEsq(sqrt(OmEsq));
 
-    volScalarField F((asin(rtOmEsq) - E*rtOmEsq)/OmEsq);
+    volScalarField F(max(asin(rtOmEsq) - E*rtOmEsq, residualE_)/OmEsq);
 
     return
         (8.0/3.0)
        *Eo
        /(
             Eo*pow(E, 2.0/3.0)/OmEsq
-          + 16*pow(E, 0.75)
+          + 16*pow(E, 4.0/3.0)
         )
        /sqr(F)
        *max(pair_.Re(), residualRe_);

applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/liftModels/Moraga/Moraga.C

@@ -62,15 +62,39 @@ Foam::liftModels::Moraga::~Moraga()
 
 Foam::tmp<Foam::volScalarField> Foam::liftModels::Moraga::Cl() const
 {
-    volScalarField ReSqrSr
+    volScalarField Re(pair_.Re());
+
+    volScalarField sqrSr
     (
-        pair_.Re()
-       *sqr(pair_.dispersed().d())
+        sqr(pair_.dispersed().d())
        /pair_.continuous().nu()
        *mag(fvc::grad(pair_.continuous().U()))
     );
 
-    return 0.2*exp(- ReSqrSr/3.6e5 - 0.12)*exp(ReSqrSr/3.0e7);
+    if
+    (
+        min(Re).value() < 1200.0
+     || max(Re).value() > 18800.0
+     || min(sqrSr).value() < 0.0016
+     || max(sqrSr).value() > 0.04
+    )
+    {
+        WarningIn
+        (
+            "Foam::tmp<Foam::volScalarField> "
+            "Foam::liftModels::Moraga::Cl() const"
+        )   << "Re and/or Sr are out of the range of applicability of the "
+            << "Moraga model. Clamping to range bounds"
+            << endl;
+    }
+
+    Re.min(1200.0);
+    Re.max(18800.0);
+
+    sqrSr.min(0.0016);
+    sqrSr.max(0.04);
+
+    return 0.2*exp(- Re*sqrSr/3.6e5 - 0.12)*exp(Re*sqrSr/3.0e7);
 }