Merge branch 'master' of github.com:OpenFOAM/OpenFOAM-2.3.x

applications/solvers/multiphase/driftFluxFoam/driftFluxFoam.C

@@ -87,7 +87,6 @@ int main(int argc, char *argv[])
             #include "alphaEqnSubCycle.H"
 
             twoPhaseProperties.correct();
-            Info<< average(twoPhaseProperties.mu()) << endl;
 
             #include "UEqn.H"
 

applications/solvers/multiphase/driftFluxFoam/kEpsilon.H

@@ -10,8 +10,6 @@ if (turbulence)
     dimensionedScalar epsilon0("epsilon0", epsilon.dimensions(), 0);
     dimensionedScalar epsilonMin("epsilonMin", epsilon.dimensions(), SMALL);
 
-    volScalarField divU(fvc::div(phi));
-
     tmp<volTensorField> tgradU = fvc::grad(U);
     volScalarField G(mut*(tgradU() && dev(twoSymm(tgradU()))));
     tgradU.clear();
@@ -21,7 +19,7 @@ if (turbulence)
         Cmu*k/sigmak*(g & fvc::grad(rho))/(epsilon + epsilonMin)
     );
 
-    volScalarField muc(twoPhaseProperties.nucModel().nu()*rho2);
+    volScalarField mul(twoPhaseProperties.mu());
 
     #include "wallFunctions.H"
 
@@ -32,12 +30,12 @@ if (turbulence)
       + fvm::div(rhoPhi, epsilon)
       - fvm::laplacian
         (
-            mut/sigmaEps + muc, epsilon,
+            mut/sigmaEps + mul, epsilon,
             "laplacian(DepsilonEff,epsilon)"
         )
      ==
         C1*G*epsilon/(k + kMin)
-      - fvm::SuSp(C1*(1.0 - C3)*Gcoef + (2.0/3.0*C1)*rho*divU, epsilon)
+      - fvm::SuSp(C1*(1.0 - C3)*Gcoef, epsilon)
       - fvm::Sp(C2*rho*epsilon/(k + kMin), epsilon)
     );
 
@@ -56,12 +54,12 @@ if (turbulence)
       + fvm::div(rhoPhi, k)
       - fvm::laplacian
         (
-            mut/sigmak + muc, k,
+            mut/sigmak + mul, k,
             "laplacian(DkEff,k)"
         )
      ==
         G
-      - fvm::SuSp(Gcoef + 2.0/3.0*rho*divU, k)
+      - fvm::SuSp(Gcoef, k)
       - fvm::Sp(rho*epsilon/(k + kMin), k)
     );
 
@@ -75,6 +73,10 @@ if (turbulence)
     mut = rho*Cmu*sqr(k)/(epsilon + epsilonMin);
 
     #include "wallViscosity.H"
-}
 
-muEff = mut + twoPhaseProperties.mu();
+    muEff = mut + mul;
+}
+else
+{
+    muEff = mut + twoPhaseProperties.mu();
+}

applications/solvers/multiphase/driftFluxFoam/wallFunctions.H

@@ -33,7 +33,7 @@
         if (isA<wallFvPatch>(curPatch))
         {
             const scalarField& mutw = mut.boundaryField()[patchi];
-            const scalarField& mucw = muc.boundaryField()[patchi];
+            const scalarField& mulw = mul.boundaryField()[patchi];
 
             scalarField magFaceGradU
             (
@@ -55,7 +55,7 @@
                    /(kappa_*y[patchi][facei]);
 
                 G[faceCelli] +=
-                    (mutw[facei] + mucw[facei])
+                    (mutw[facei] + mulw[facei])
                    *magFaceGradU[facei]
                    *Cmu25*::sqrt(k[faceCelli])
                    /(kappa_*y[patchi][facei]);

applications/solvers/multiphase/driftFluxFoam/wallViscosity.H

@@ -12,7 +12,8 @@
         if (isA<wallFvPatch>(curPatch))
         {
             scalarField& mutw = mut.boundaryField()[patchi];
-            const scalarField& mucw = muc.boundaryField()[patchi];
+            const scalarField& mulw = mul.boundaryField()[patchi];
+            const scalarField& rhow = rho.boundaryField()[patchi];
 
             forAll(curPatch, facei)
             {
@@ -20,12 +21,12 @@
 
                 scalar yPlus =
                     Cmu25*y[patchi][facei]*::sqrt(k[faceCelli])
-                   /(mucw[facei]/rho2.value());
+                   /(mulw[facei]/rhow[facei]);
 
                 if (yPlus > 11.6)
                 {
                     mutw[facei] =
-                        mucw[facei]*(yPlus*kappa_/::log(E_*yPlus) - 1);
+                        mulw[facei]*(yPlus*kappa_/::log(E_*yPlus) - 1);
                 }
                 else
                 {

src/TurbulenceModels/phaseIncompressible/RAS/continuousGasKEpsilon/continuousGasKEpsilon.C

@@ -124,9 +124,17 @@ void continuousGasKEpsilon<BasicTurbulenceModel>::correctNut()
     const transportModel& liquid = fluid.otherPhase(gas);
 
     volScalarField thetal(liquidTurbulence.k()/liquidTurbulence.epsilon());
-    volScalarField thetag((1.0/(18*liquid.nu()))*sqr(gas.d()));
-    volScalarField expThetar(exp(min(thetal/thetag, scalar(50))));
-    volScalarField omega(sqr(expThetar - 1)/(sqr(expThetar) - 1));
+    volScalarField rhodv(gas.rho() + fluid.virtualMass(gas).Cvm()*liquid.rho());
+    volScalarField thetag((rhodv/(18*liquid.rho()*liquid.nu()))*sqr(gas.d()));
+    volScalarField expThetar
+    (
+        min
+        (
+            exp(min(thetal/thetag, scalar(50))),
+            scalar(1)
+        )
+    );
+    volScalarField omega((1 - expThetar)/(1 + expThetar));
 
     nutEff_ = omega*liquidTurbulence.nut();
 }