bubbleFoam, twoPhaseEulerFoam: rationalise kEpsilon

applications/solvers/multiphase/bubbleFoam/bubbleFoam.C

@@ -90,6 +90,7 @@ int main(int argc, char *argv[])
             if (pimple.turbCorr())
             {
                 #include "kEpsilon.H"
+                nuEffa = sqr(Ct)*nutb + nua;
             }
         }
 

applications/solvers/multiphase/bubbleFoam/kEpsilon.H

@@ -14,16 +14,17 @@ if (turbulence)
     // Dissipation equation
     fvScalarMatrix epsEqn
     (
-        fvm::ddt(beta, epsilon)
+        fvm::ddt(epsilon)
       + fvm::div(phib, epsilon)
+      - fvm::Sp(fvc::div(phib), epsilon)
       - fvm::laplacian
         (
             alphaEps*nuEffb, epsilon,
             "laplacian(DepsilonEff,epsilon)"
         )
       ==
-         C1*beta*G*epsilon/k
+         C1*G*epsilon/k
-       - fvm::Sp(C2*beta*epsilon/k, epsilon)
+       - fvm::Sp(C2*epsilon/k, epsilon)
     );
 
     #include "wallDissipation.H"
@@ -37,16 +38,17 @@ if (turbulence)
     // Turbulent kinetic energy equation
     fvScalarMatrix kEqn
     (
-        fvm::ddt(beta, k)
+        fvm::ddt(k)
       + fvm::div(phib, k)
+      - fvm::Sp(fvc::div(phib), k)
       - fvm::laplacian
         (
             alphak*nuEffb, k,
             "laplacian(DkEff,k)"
         )
       ==
-        beta*G
+        G
-      - fvm::Sp(beta*epsilon/k, k)
+      - fvm::Sp(epsilon/k, k)
     );
     kEqn.relax();
     kEqn.solve();
@@ -59,5 +61,4 @@ if (turbulence)
     #include "wallViscosity.H"
 }
 
-nuEffa = sqr(Ct)*nutb + nua;
 nuEffb = nutb + nub;

applications/solvers/multiphase/twoPhaseEulerFoam/kEpsilon.H

@@ -1,64 +0,0 @@
-if (turbulence)
-{
-    if (mesh.changing())
-    {
-        y.correct();
-    }
-
-    tmp<volTensorField> tgradUb = fvc::grad(Ub);
-    volScalarField G(2*nutb*(tgradUb() && dev(symm(tgradUb()))));
-    tgradUb.clear();
-
-    #include "wallFunctions.H"
-
-    // Dissipation equation
-    fvScalarMatrix epsEqn
-    (
-        fvm::ddt(epsilon)
-      + fvm::div(phib, epsilon)
-      - fvm::Sp(fvc::div(phib), epsilon)
-      - fvm::laplacian
-        (
-            alphaEps*nuEffb, epsilon,
-            "laplacian(DepsilonEff,epsilon)"
-        )
-      ==
-         C1*G*epsilon/k
-       - fvm::Sp(C2*epsilon/k, epsilon)
-    );
-
-    #include "wallDissipation.H"
-
-    epsEqn.relax();
-    epsEqn.solve();
-
-    epsilon.max(dimensionedScalar("zero", epsilon.dimensions(), 1.0e-15));
-
-
-    // Turbulent kinetic energy equation
-    fvScalarMatrix kEqn
-    (
-        fvm::ddt(k)
-      + fvm::div(phib, k)
-      - fvm::Sp(fvc::div(phib), k)
-      - fvm::laplacian
-        (
-            alphak*nuEffb, k,
-            "laplacian(DkEff,k)"
-        )
-      ==
-        G
-      - fvm::Sp(epsilon/k, k)
-    );
-    kEqn.relax();
-    kEqn.solve();
-
-    k.max(dimensionedScalar("zero", k.dimensions(), 1.0e-8));
-
-    //- Re-calculate turbulence viscosity
-    nutb = Cmu*sqr(k)/epsilon;
-
-    #include "wallViscosity.H"
-}
-
-nuEffb = nutb + nub;