From 1464c4ff5c6b568ed6bf28475dfc8ae79dff741e Mon Sep 17 00:00:00 2001
From: Henry <Henry>
Date: Fri, 27 Jul 2012 14:56:01 +0100
Subject: [PATCH] multiphase (VoF): Added support for general turbulence models
 Required the addition of the divDevRhoR function to all incompressible
 turbulence models

---
 .../solvers/multiphase/cavitatingFoam/UEqn.H  | 11 +----
 .../multiphase/compressibleInterFoam/UEqn.H   | 11 +----
 .../solvers/multiphase/interFoam/UEqn.H       | 11 +----
 .../multiphase/interPhaseChangeFoam/UEqn.H    | 11 +----
 .../multiphase/multiphaseInterFoam/UEqn.H     | 11 +----
 .../solvers/multiphase/settlingFoam/UEqn.H    |  1 -
 .../multiphase/twoLiquidMixingFoam/UEqn.H     | 11 +----
 .../functionObjects/forces/forces/forces.C    |  4 +-
 .../LES/GenEddyVisc/GenEddyVisc.C             | 23 ++++++++--
 .../LES/GenEddyVisc/GenEddyVisc.H             | 12 ++++-
 .../LES/GenSGSStress/GenSGSStress.C           | 36 ++++++++++++++-
 .../LES/GenSGSStress/GenSGSStress.H           | 16 +++++--
 .../incompressible/LES/LESModel/LESModel.H    | 20 ---------
 .../LES/Smagorinsky2/Smagorinsky2.C           | 24 +++++++++-
 .../LES/Smagorinsky2/Smagorinsky2.H           | 14 ++++--
 .../LES/SpalartAllmaras/SpalartAllmaras.C     | 23 ++++++++--
 .../LES/SpalartAllmaras/SpalartAllmaras.H     | 16 +++++--
 .../LES/kOmegaSSTSAS/kOmegaSSTSAS.C           | 23 ++++++++--
 .../LES/kOmegaSSTSAS/kOmegaSSTSAS.H           | 16 +++++--
 .../incompressible/LES/laminar/laminar.C      | 23 ++++++++--
 .../incompressible/LES/laminar/laminar.H      | 16 +++++--
 .../LES/mixedSmagorinsky/mixedSmagorinsky.C   | 26 ++++++++---
 .../LES/mixedSmagorinsky/mixedSmagorinsky.H   | 16 +++++--
 .../LES/scaleSimilarity/scaleSimilarity.C     | 16 +++++--
 .../LES/scaleSimilarity/scaleSimilarity.H     | 16 +++++--
 .../incompressible/RAS/LRR/LRR.C              | 38 ++++++++++++++++
 .../incompressible/RAS/LRR/LRR.H              |  7 +++
 .../RAS/LamBremhorstKE/LamBremhorstKE.C       | 16 +++++++
 .../RAS/LamBremhorstKE/LamBremhorstKE.H       |  7 +++
 .../RAS/LaunderGibsonRSTM/LaunderGibsonRSTM.C | 45 ++++++++++++++++++-
 .../RAS/LaunderGibsonRSTM/LaunderGibsonRSTM.H |  7 +++
 .../RAS/LaunderSharmaKE/LaunderSharmaKE.C     | 16 +++++++
 .../RAS/LaunderSharmaKE/LaunderSharmaKE.H     |  7 +++
 .../RAS/LienCubicKE/LienCubicKE.C             | 17 +++++++
 .../RAS/LienCubicKE/LienCubicKE.H             |  7 +++
 .../RAS/LienCubicKELowRe/LienCubicKELowRe.C   | 17 +++++++
 .../RAS/LienCubicKELowRe/LienCubicKELowRe.H   |  7 +++
 .../LienLeschzinerLowRe/LienLeschzinerLowRe.C | 17 ++++++-
 .../LienLeschzinerLowRe/LienLeschzinerLowRe.H |  7 +++
 .../RAS/NonlinearKEShih/NonlinearKEShih.C     | 17 +++++++
 .../RAS/NonlinearKEShih/NonlinearKEShih.H     |  7 +++
 .../RAS/RNGkEpsilon/RNGkEpsilon.C             | 16 +++++++
 .../RAS/RNGkEpsilon/RNGkEpsilon.H             |  7 +++
 .../RAS/SpalartAllmaras/SpalartAllmaras.C     | 16 +++++++
 .../RAS/SpalartAllmaras/SpalartAllmaras.H     |  7 +++
 .../incompressible/RAS/kEpsilon/kEpsilon.C    | 16 +++++++
 .../incompressible/RAS/kEpsilon/kEpsilon.H    |  7 +++
 .../incompressible/RAS/kOmega/kOmega.C        | 16 +++++++
 .../incompressible/RAS/kOmega/kOmega.H        |  7 +++
 .../incompressible/RAS/kOmegaSST/kOmegaSST.C  | 16 +++++++
 .../incompressible/RAS/kOmegaSST/kOmegaSST.H  |  7 +++
 .../incompressible/RAS/kkLOmega/kkLOmega.C    | 16 +++++++
 .../incompressible/RAS/kkLOmega/kkLOmega.H    |  7 +++
 .../incompressible/RAS/laminar/laminar.C      | 16 +++++++
 .../incompressible/RAS/laminar/laminar.H      |  7 +++
 .../incompressible/RAS/qZeta/qZeta.C          | 16 +++++++
 .../incompressible/RAS/qZeta/qZeta.H          |  7 +++
 .../RAS/realizableKE/realizableKE.C           | 16 +++++++
 .../RAS/realizableKE/realizableKE.H           |  7 +++
 .../turbulenceModel/laminar/laminar.C         | 16 +++++++
 .../turbulenceModel/laminar/laminar.H         |  7 +++
 .../turbulenceModel/turbulenceModel.H         |  7 +++
 .../LTSInterFoam/wigleyHull/system/fvSchemes  |  1 +
 .../mixerVessel2D/system/fvSchemes            |  1 +
 .../mixerVessel2D/system/fvSchemes            |  1 +
 .../les/throttle/system/fvSchemes             |  1 +
 .../les/throttle3D/system/fvSchemes           |  1 +
 .../ras/throttle/system/fvSchemes             |  1 +
 .../les/depthCharge2D/system/fvSchemes        |  2 +-
 .../les/depthCharge3D/system/fvSchemes        |  2 +-
 .../ras/damBreakWithObstacle/system/fvSchemes |  1 +
 .../ras/floatingObject/system/fvSchemes       |  1 +
 .../ras/sloshingTank2D/system/fvSchemes       |  1 +
 .../ras/sloshingTank2D3DoF/system/fvSchemes   |  1 +
 .../ras/sloshingTank3D/system/fvSchemes       |  1 +
 .../ras/sloshingTank3D3DoF/system/fvSchemes   |  1 +
 .../ras/sloshingTank3D6DoF/system/fvSchemes   |  1 +
 .../ras/testTubeMixer/system/fvSchemes        |  1 +
 .../laminar/capillaryRise/system/fvSchemes    |  1 +
 .../laminar/damBreak/system/fvSchemes         |  1 +
 .../les/nozzleFlow2D/system/fvSchemes         |  2 +-
 .../interFoam/ras/damBreak/system/fvSchemes   |  2 +-
 .../ras/damBreakPorousBaffle/system/fvSchemes |  2 +-
 .../LTSInterFoam/system/fvSchemes             |  1 +
 .../ras/waterChannel/system/fvSchemes         |  1 +
 .../ras/weirOverflow/system/fvSchemes         |  2 +-
 .../laminar/damBreak/system/fvSchemes         |  1 +
 .../cavitatingBullet/system/fvSchemes         |  2 +-
 .../laminar/damBreak4phase/system/fvSchemes   |  1 +
 .../damBreak4phaseFine/system/fvSchemes       |  1 +
 .../lockExchange/system/fvSchemes             |  2 +-
 91 files changed, 774 insertions(+), 146 deletions(-)

diff --git a/applications/solvers/multiphase/cavitatingFoam/UEqn.H b/applications/solvers/multiphase/cavitatingFoam/UEqn.H
index 9a5761b59f3..1ca0b9f0a5f 100644
--- a/applications/solvers/multiphase/cavitatingFoam/UEqn.H
+++ b/applications/solvers/multiphase/cavitatingFoam/UEqn.H
@@ -1,17 +1,8 @@
-    surfaceScalarField muEff
-    (
-        "muEff",
-        twoPhaseProperties.muf()
-      + fvc::interpolate(rho*turbulence->nut())
-    );
-
     fvVectorMatrix UEqn
     (
         fvm::ddt(rho, U)
       + fvm::div(phi, U)
-      - fvm::laplacian(muEff, U)
-    //- (fvc::grad(U) & fvc::grad(muf))
-      - fvc::div(muEff*(fvc::interpolate(dev(fvc::grad(U))) & mesh.Sf()))
+      + turbulence->divDevRhoReff(rho, U)
     );
 
     UEqn.relax();
diff --git a/applications/solvers/multiphase/compressibleInterFoam/UEqn.H b/applications/solvers/multiphase/compressibleInterFoam/UEqn.H
index 257f6d48b53..7cc250a66a2 100644
--- a/applications/solvers/multiphase/compressibleInterFoam/UEqn.H
+++ b/applications/solvers/multiphase/compressibleInterFoam/UEqn.H
@@ -1,17 +1,8 @@
-    surfaceScalarField muEff
-    (
-        "muEff",
-        twoPhaseProperties.muf()
-      + fvc::interpolate(rho*turbulence->nut())
-    );
-
     fvVectorMatrix UEqn
     (
         fvm::ddt(rho, U)
       + fvm::div(rhoPhi, U)
-      - fvm::laplacian(muEff, U)
-      - (fvc::grad(U) & fvc::grad(muEff))
-    //- fvc::div(muEff*(fvc::interpolate(dev(fvc::grad(U))) & mesh.Sf()))
+      + turbulence->divDevRhoReff(rho, U)
     );
 
     UEqn.relax();
diff --git a/applications/solvers/multiphase/interFoam/UEqn.H b/applications/solvers/multiphase/interFoam/UEqn.H
index 257f6d48b53..7cc250a66a2 100644
--- a/applications/solvers/multiphase/interFoam/UEqn.H
+++ b/applications/solvers/multiphase/interFoam/UEqn.H
@@ -1,17 +1,8 @@
-    surfaceScalarField muEff
-    (
-        "muEff",
-        twoPhaseProperties.muf()
-      + fvc::interpolate(rho*turbulence->nut())
-    );
-
     fvVectorMatrix UEqn
     (
         fvm::ddt(rho, U)
       + fvm::div(rhoPhi, U)
-      - fvm::laplacian(muEff, U)
-      - (fvc::grad(U) & fvc::grad(muEff))
-    //- fvc::div(muEff*(fvc::interpolate(dev(fvc::grad(U))) & mesh.Sf()))
+      + turbulence->divDevRhoReff(rho, U)
     );
 
     UEqn.relax();
diff --git a/applications/solvers/multiphase/interPhaseChangeFoam/UEqn.H b/applications/solvers/multiphase/interPhaseChangeFoam/UEqn.H
index b8fe82ff142..062c5523c9d 100644
--- a/applications/solvers/multiphase/interPhaseChangeFoam/UEqn.H
+++ b/applications/solvers/multiphase/interPhaseChangeFoam/UEqn.H
@@ -1,18 +1,9 @@
-    surfaceScalarField muEff
-    (
-        "muEff",
-        twoPhaseProperties->muf()
-      + fvc::interpolate(rho*turbulence->nut())
-    );
-
     fvVectorMatrix UEqn
     (
         fvm::ddt(rho, U)
       + fvm::div(rhoPhi, U)
       - fvm::Sp(fvc::ddt(rho) + fvc::div(rhoPhi), U)
-      - fvm::laplacian(muEff, U)
-      - (fvc::grad(U) & fvc::grad(muEff))
-    //- fvc::div(muEff*(fvc::interpolate(dev2(fvc::grad(U))) & mesh.Sf()))
+      + turbulence->divDevRhoReff(rho, U)
     );
 
     UEqn.relax();
diff --git a/applications/solvers/multiphase/multiphaseInterFoam/UEqn.H b/applications/solvers/multiphase/multiphaseInterFoam/UEqn.H
index 43ee1e3eb94..fbcba7db72d 100644
--- a/applications/solvers/multiphase/multiphaseInterFoam/UEqn.H
+++ b/applications/solvers/multiphase/multiphaseInterFoam/UEqn.H
@@ -1,17 +1,8 @@
-    surfaceScalarField muEff
-    (
-        "muEff",
-        mixture.muf()
-      + fvc::interpolate(rho*turbulence->nut())
-    );
-
     fvVectorMatrix UEqn
     (
         fvm::ddt(rho, U)
       + fvm::div(mixture.rhoPhi(), U)
-      - fvm::laplacian(muEff, U)
-      - (fvc::grad(U) & fvc::grad(muEff))
-    //- fvc::div(muEff*(fvc::interpolate(dev(fvc::grad(U))) & mesh.Sf()))
+      + turbulence->divDevRhoReff(rho, U)
     );
 
     UEqn.relax();
diff --git a/applications/solvers/multiphase/settlingFoam/UEqn.H b/applications/solvers/multiphase/settlingFoam/UEqn.H
index d6232da309d..ff6323bca27 100644
--- a/applications/solvers/multiphase/settlingFoam/UEqn.H
+++ b/applications/solvers/multiphase/settlingFoam/UEqn.H
@@ -11,7 +11,6 @@
         )
       - fvm::laplacian(muEff, U, "laplacian(muEff,U)")
       - (fvc::grad(U) & fvc::grad(muEff))
-    //- fvc::div(muEff*dev2(T(fvc::grad(U))))
     );
 
     UEqn.relax();
diff --git a/applications/solvers/multiphase/twoLiquidMixingFoam/UEqn.H b/applications/solvers/multiphase/twoLiquidMixingFoam/UEqn.H
index 3dbb3747743..7afd323dd0d 100644
--- a/applications/solvers/multiphase/twoLiquidMixingFoam/UEqn.H
+++ b/applications/solvers/multiphase/twoLiquidMixingFoam/UEqn.H
@@ -1,17 +1,8 @@
-    surfaceScalarField muEff
-    (
-        "muEff",
-        twoPhaseProperties.muf()
-      + fvc::interpolate(rho*turbulence->nut())
-    );
-
     fvVectorMatrix UEqn
     (
         fvm::ddt(rho, U)
       + fvm::div(rhoPhi, U)
-      - fvm::laplacian(muEff, U)
-      - (fvc::grad(U) & fvc::grad(muEff))
-    //- fvc::div(muEff*(fvc::interpolate(dev(fvc::grad(U))) & mesh.Sf()))
+      + turbulence->divDevRhoReff(rho, U)
     );
 
     UEqn.relax();
diff --git a/src/postProcessing/functionObjects/forces/forces/forces.C b/src/postProcessing/functionObjects/forces/forces/forces.C
index 1dd26906b74..ee8b8bb1b38 100644
--- a/src/postProcessing/functionObjects/forces/forces/forces.C
+++ b/src/postProcessing/functionObjects/forces/forces/forces.C
@@ -66,14 +66,14 @@ Foam::tmp<Foam::volSymmTensorField> Foam::forces::devRhoReff() const
         const compressible::LESModel& les =
         obr_.lookupObject<compressible::LESModel>("LESProperties");
 
-        return les.devRhoBeff();
+        return les.devRhoReff();
     }
     else if (obr_.foundObject<incompressible::LESModel>("LESProperties"))
     {
         const incompressible::LESModel& les
             = obr_.lookupObject<incompressible::LESModel>("LESProperties");
 
-        return rho()*les.devBeff();
+        return rho()*les.devReff();
     }
     else if (obr_.foundObject<basicThermo>("thermophysicalProperties"))
     {
diff --git a/src/turbulenceModels/incompressible/LES/GenEddyVisc/GenEddyVisc.C b/src/turbulenceModels/incompressible/LES/GenEddyVisc/GenEddyVisc.C
index 0c905a91ced..e4bde9ce7b6 100644
--- a/src/turbulenceModels/incompressible/LES/GenEddyVisc/GenEddyVisc.C
+++ b/src/turbulenceModels/incompressible/LES/GenEddyVisc/GenEddyVisc.C
@@ -87,17 +87,34 @@ tmp<volSymmTensorField> GenEddyVisc::B() const
 }
 
 
-tmp<volSymmTensorField> GenEddyVisc::devBeff() const
+tmp<volSymmTensorField> GenEddyVisc::devReff() const
 {
     return -nuEff()*dev(twoSymm(fvc::grad(U())));
 }
 
 
-tmp<fvVectorMatrix> GenEddyVisc::divDevBeff(volVectorField& U) const
+tmp<fvVectorMatrix> GenEddyVisc::divDevReff(volVectorField& U) const
 {
     return
     (
-      - fvm::laplacian(nuEff(), U) - fvc::div(nuEff()*dev(T(fvc::grad(U))))
+      - fvm::laplacian(nuEff(), U)
+      - fvc::div(nuEff()*dev(T(fvc::grad(U))))
+    );
+}
+
+
+tmp<fvVectorMatrix> GenEddyVisc::divDevRhoReff
+(
+    const volScalarField& rho,
+    volVectorField& U
+) const
+{
+    volScalarField muEff("muEff", rho*nuEff());
+
+    return
+    (
+      - fvm::laplacian(muEff, U)
+      - fvc::div(muEff*dev(T(fvc::grad(U))))
     );
 }
 
diff --git a/src/turbulenceModels/incompressible/LES/GenEddyVisc/GenEddyVisc.H b/src/turbulenceModels/incompressible/LES/GenEddyVisc/GenEddyVisc.H
index 0746377eb79..105f62c6647 100644
--- a/src/turbulenceModels/incompressible/LES/GenEddyVisc/GenEddyVisc.H
+++ b/src/turbulenceModels/incompressible/LES/GenEddyVisc/GenEddyVisc.H
@@ -118,11 +118,19 @@ public:
 
         //- Return the effective sub-grid turbulence stress tensor
         //  including the laminar stress
-        virtual tmp<volSymmTensorField> devBeff() const;
+        virtual tmp<volSymmTensorField> devReff() const;
 
         //- Return the deviatoric part of the effective sub-grid
         //  turbulence stress tensor including the laminar stress
-        virtual tmp<fvVectorMatrix> divDevBeff(volVectorField& U) const;
+        virtual tmp<fvVectorMatrix> divDevReff(volVectorField& U) const;
+
+        //- Return the deviatoric part of the effective sub-grid
+        //  turbulence stress tensor including the laminar stress
+        virtual tmp<fvVectorMatrix> divDevRhoReff
+        (
+            const volScalarField& rho,
+            volVectorField& U
+        ) const;
 
         //- Correct Eddy-Viscosity and related properties
         virtual void correct(const tmp<volTensorField>& gradU);
diff --git a/src/turbulenceModels/incompressible/LES/GenSGSStress/GenSGSStress.C b/src/turbulenceModels/incompressible/LES/GenSGSStress/GenSGSStress.C
index 9d385c21eb1..cbc803da9dd 100644
--- a/src/turbulenceModels/incompressible/LES/GenSGSStress/GenSGSStress.C
+++ b/src/turbulenceModels/incompressible/LES/GenSGSStress/GenSGSStress.C
@@ -115,7 +115,7 @@ GenSGSStress::GenSGSStress
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-tmp<volSymmTensorField> GenSGSStress::devBeff() const
+tmp<volSymmTensorField> GenSGSStress::devReff() const
 {
     return tmp<volSymmTensorField>
     (
@@ -135,7 +135,7 @@ tmp<volSymmTensorField> GenSGSStress::devBeff() const
 }
 
 
-tmp<fvVectorMatrix> GenSGSStress::divDevBeff
+tmp<fvVectorMatrix> GenSGSStress::divDevReff
 (
     volVectorField& U
 ) const
@@ -164,6 +164,38 @@ tmp<fvVectorMatrix> GenSGSStress::divDevBeff
 }
 
 
+tmp<fvVectorMatrix> GenSGSStress::divDevRhoReff
+(
+    const volScalarField& rho,
+    volVectorField& U
+) const
+{
+    volScalarField muEff("muEff", rho*nuEff());
+
+    if (couplingFactor_.value() > 0.0)
+    {
+        return
+        (
+            fvc::div(rho*B_ + couplingFactor_*rho*nuSgs_*fvc::grad(U))
+          + fvc::laplacian
+            (
+                (1.0 - couplingFactor_)*rho*nuSgs_, U, "laplacian(muEff,U)"
+            )
+          - fvm::laplacian(muEff, U)
+        );
+    }
+    else
+    {
+        return
+        (
+            fvc::div(rho*B_)
+          + fvc::laplacian(rho*nuSgs_, U, "laplacian(muEff,U)")
+          - fvm::laplacian(muEff, U)
+        );
+    }
+}
+
+
 bool GenSGSStress::read()
 {
     if (LESModel::read())
diff --git a/src/turbulenceModels/incompressible/LES/GenSGSStress/GenSGSStress.H b/src/turbulenceModels/incompressible/LES/GenSGSStress/GenSGSStress.H
index 36d84f9c771..6b805f70c2e 100644
--- a/src/turbulenceModels/incompressible/LES/GenSGSStress/GenSGSStress.H
+++ b/src/turbulenceModels/incompressible/LES/GenSGSStress/GenSGSStress.H
@@ -128,11 +128,19 @@ public:
 
         //- Return the effective sub-grid turbulence stress tensor
         //  including the laminar stress
-        virtual tmp<volSymmTensorField> devBeff() const;
+        virtual tmp<volSymmTensorField> devReff() const;
 
-        //- Returns div(B).
-        // This is the additional term due to the filtering of the NSE.
-        virtual tmp<fvVectorMatrix> divDevBeff(volVectorField& U) const;
+        //- Return the deviatoric part of the effective sub-grid
+        //  turbulence stress tensor including the laminar stress
+        virtual tmp<fvVectorMatrix> divDevReff(volVectorField& U) const;
+
+        //- Return the deviatoric part of the effective sub-grid
+        //  turbulence stress tensor including the laminar stress
+        virtual tmp<fvVectorMatrix> divDevRhoReff
+        (
+            const volScalarField& rho,
+            volVectorField& U
+        ) const;
 
         //- Read LESProperties dictionary
         virtual bool read();
diff --git a/src/turbulenceModels/incompressible/LES/LESModel/LESModel.H b/src/turbulenceModels/incompressible/LES/LESModel/LESModel.H
index deca26e7936..3688cb8dec9 100644
--- a/src/turbulenceModels/incompressible/LES/LESModel/LESModel.H
+++ b/src/turbulenceModels/incompressible/LES/LESModel/LESModel.H
@@ -202,14 +202,6 @@ public:
         //- Return the sub-grid stress tensor.
         virtual tmp<volSymmTensorField> B() const = 0;
 
-        //- Return the deviatoric part of the effective sub-grid
-        //  turbulence stress tensor including the laminar stress
-        virtual tmp<volSymmTensorField> devBeff() const = 0;
-
-        //- Returns div(dev(Beff)).
-        //  This is the additional term due to the filtering of the NSE.
-        virtual tmp<fvVectorMatrix> divDevBeff(volVectorField& U) const = 0;
-
 
         // RAS compatibility functions for the turbulenceModel base class
 
@@ -225,18 +217,6 @@ public:
                 return B();
             }
 
-            //- Return the effective stress tensor including the laminar stress
-            virtual tmp<volSymmTensorField> devReff() const
-            {
-                return devBeff();
-            }
-
-            //- Return the source term for the momentum equation
-            virtual tmp<fvVectorMatrix> divDevReff(volVectorField& U) const
-            {
-                return divDevBeff(U);
-            }
-
 
         //- Correct Eddy-Viscosity and related properties.
         //  This calls correct(const tmp<volTensorField>& gradU) by supplying
diff --git a/src/turbulenceModels/incompressible/LES/Smagorinsky2/Smagorinsky2.C b/src/turbulenceModels/incompressible/LES/Smagorinsky2/Smagorinsky2.C
index 97d32bce9dd..f53cf08282d 100644
--- a/src/turbulenceModels/incompressible/LES/Smagorinsky2/Smagorinsky2.C
+++ b/src/turbulenceModels/incompressible/LES/Smagorinsky2/Smagorinsky2.C
@@ -81,7 +81,7 @@ tmp<volSymmTensorField> Smagorinsky2::B() const
 }
 
 
-tmp<fvVectorMatrix> Smagorinsky2::divDevBeff
+tmp<fvVectorMatrix> Smagorinsky2::divDevReff
 (
     volVectorField& U
 ) const
@@ -101,6 +101,28 @@ tmp<fvVectorMatrix> Smagorinsky2::divDevBeff
 }
 
 
+tmp<fvVectorMatrix> Smagorinsky2::divDevRhoReff
+(
+    const volScalarField& rho,
+    volVectorField& U
+) const
+{
+    volTensorField gradU(fvc::grad(U));
+
+    volSymmTensorField aniMuEff
+    (
+        "muEff",
+        I*(rho*nuEff()) + (cD2_*rho*delta())*symm(gradU)
+    );
+
+    return
+    (
+      - fvm::laplacian(aniMuEff, U)
+      - fvc::div(rho*nuEff()*dev(T(fvc::grad(U))))
+    );
+}
+
+
 bool Smagorinsky2::read()
 {
     if (Smagorinsky::read())
diff --git a/src/turbulenceModels/incompressible/LES/Smagorinsky2/Smagorinsky2.H b/src/turbulenceModels/incompressible/LES/Smagorinsky2/Smagorinsky2.H
index 89c77912cde..bcce2151bf3 100644
--- a/src/turbulenceModels/incompressible/LES/Smagorinsky2/Smagorinsky2.H
+++ b/src/turbulenceModels/incompressible/LES/Smagorinsky2/Smagorinsky2.H
@@ -109,9 +109,17 @@ public:
         //- Return B.
         virtual tmp<volSymmTensorField> B() const;
 
-        //- Returns div(B).
-        // This is the additional term due to the filtering of the NSE.
-        virtual tmp<fvVectorMatrix> divDevBeff(volVectorField& U) const;
+        //- Return the deviatoric part of the effective sub-grid
+        //  turbulence stress tensor including the laminar stress
+        virtual tmp<fvVectorMatrix> divDevReff(volVectorField& U) const;
+
+        //- Return the deviatoric part of the effective sub-grid
+        //  turbulence stress tensor including the laminar stress
+        virtual tmp<fvVectorMatrix> divDevRhoReff
+        (
+            const volScalarField& rho,
+            volVectorField& U
+        ) const;
 
         //- Read LESProperties dictionary
         virtual bool read();
diff --git a/src/turbulenceModels/incompressible/LES/SpalartAllmaras/SpalartAllmaras.C b/src/turbulenceModels/incompressible/LES/SpalartAllmaras/SpalartAllmaras.C
index 6a34f90892e..dd16e0ee93b 100644
--- a/src/turbulenceModels/incompressible/LES/SpalartAllmaras/SpalartAllmaras.C
+++ b/src/turbulenceModels/incompressible/LES/SpalartAllmaras/SpalartAllmaras.C
@@ -338,17 +338,34 @@ tmp<volSymmTensorField> SpalartAllmaras::B() const
 }
 
 
-tmp<volSymmTensorField> SpalartAllmaras::devBeff() const
+tmp<volSymmTensorField> SpalartAllmaras::devReff() const
 {
     return -nuEff()*dev(twoSymm(fvc::grad(U())));
 }
 
 
-tmp<fvVectorMatrix> SpalartAllmaras::divDevBeff(volVectorField& U) const
+tmp<fvVectorMatrix> SpalartAllmaras::divDevReff(volVectorField& U) const
 {
     return
     (
-      - fvm::laplacian(nuEff(), U) - fvc::div(nuEff()*dev(T(fvc::grad(U))))
+      - fvm::laplacian(nuEff(), U)
+      - fvc::div(nuEff()*dev(T(fvc::grad(U))))
+    );
+}
+
+
+tmp<fvVectorMatrix> SpalartAllmaras::divDevRhoReff
+(
+    const volScalarField& rho,
+    volVectorField& U
+) const
+{
+    volScalarField muEff("muEff", rho*nuEff());
+
+    return
+    (
+      - fvm::laplacian(muEff, U)
+      - fvc::div(muEff*dev(T(fvc::grad(U))))
     );
 }
 
diff --git a/src/turbulenceModels/incompressible/LES/SpalartAllmaras/SpalartAllmaras.H b/src/turbulenceModels/incompressible/LES/SpalartAllmaras/SpalartAllmaras.H
index 7ed507712a2..6fd50a2772d 100644
--- a/src/turbulenceModels/incompressible/LES/SpalartAllmaras/SpalartAllmaras.H
+++ b/src/turbulenceModels/incompressible/LES/SpalartAllmaras/SpalartAllmaras.H
@@ -172,11 +172,19 @@ public:
 
         //- Return the effective sub-grid turbulence stress tensor
         //  including the laminar stress
-        virtual tmp<volSymmTensorField> devBeff() const;
+        virtual tmp<volSymmTensorField> devReff() const;
 
-        //- Return the deviatoric part of the divergence of Beff
-        //  i.e. the additional term in the filtered NSE.
-        virtual tmp<fvVectorMatrix> divDevBeff(volVectorField& U) const;
+        //- Return the deviatoric part of the effective sub-grid
+        //  turbulence stress tensor including the laminar stress
+        virtual tmp<fvVectorMatrix> divDevReff(volVectorField& U) const;
+
+        //- Return the deviatoric part of the effective sub-grid
+        //  turbulence stress tensor including the laminar stress
+        virtual tmp<fvVectorMatrix> divDevRhoReff
+        (
+            const volScalarField& rho,
+            volVectorField& U
+        ) const;
 
         //- Correct nuTilda and related properties
         virtual void correct(const tmp<volTensorField>& gradU);
diff --git a/src/turbulenceModels/incompressible/LES/kOmegaSSTSAS/kOmegaSSTSAS.C b/src/turbulenceModels/incompressible/LES/kOmegaSSTSAS/kOmegaSSTSAS.C
index 2ab771ac737..e676731654d 100644
--- a/src/turbulenceModels/incompressible/LES/kOmegaSSTSAS/kOmegaSSTSAS.C
+++ b/src/turbulenceModels/incompressible/LES/kOmegaSSTSAS/kOmegaSSTSAS.C
@@ -425,17 +425,34 @@ tmp<volSymmTensorField> kOmegaSSTSAS::B() const
 }
 
 
-tmp<volSymmTensorField> kOmegaSSTSAS::devBeff() const
+tmp<volSymmTensorField> kOmegaSSTSAS::devReff() const
 {
     return -nuEff()*dev(twoSymm(fvc::grad(U())));
 }
 
 
-tmp<fvVectorMatrix> kOmegaSSTSAS::divDevBeff(volVectorField& U) const
+tmp<fvVectorMatrix> kOmegaSSTSAS::divDevReff(volVectorField& U) const
 {
     return
     (
-      - fvm::laplacian(nuEff(), U) - fvc::div(nuEff()*dev(T(fvc::grad(U))))
+      - fvm::laplacian(nuEff(), U)
+      - fvc::div(nuEff()*dev(T(fvc::grad(U))))
+    );
+}
+
+
+tmp<fvVectorMatrix> kOmegaSSTSAS::divDevRhoReff
+(
+    const volScalarField& rho,
+    volVectorField& U
+) const
+{
+    volScalarField muEff("muEff", rho*nuEff());
+
+    return
+    (
+      - fvm::laplacian(muEff, U)
+      - fvc::div(muEff*dev(T(fvc::grad(U))))
     );
 }
 
diff --git a/src/turbulenceModels/incompressible/LES/kOmegaSSTSAS/kOmegaSSTSAS.H b/src/turbulenceModels/incompressible/LES/kOmegaSSTSAS/kOmegaSSTSAS.H
index 2ca98380afa..1b319429e59 100644
--- a/src/turbulenceModels/incompressible/LES/kOmegaSSTSAS/kOmegaSSTSAS.H
+++ b/src/turbulenceModels/incompressible/LES/kOmegaSSTSAS/kOmegaSSTSAS.H
@@ -249,11 +249,19 @@ public:
 
         //- Return the effective sub-grid turbulence stress tensor
         //  including the laminar stress
-        virtual tmp<volSymmTensorField> devBeff() const;
+        virtual tmp<volSymmTensorField> devReff() const;
 
-        //- Return the deviatoric part of the divergence of Beff
-        //  i.e. the additional term in the filtered NSE.
-        virtual tmp<fvVectorMatrix> divDevBeff(volVectorField& U) const;
+        //- Return the deviatoric part of the effective sub-grid
+        //  turbulence stress tensor including the laminar stress
+        virtual tmp<fvVectorMatrix> divDevReff(volVectorField& U) const;
+
+        //- Return the deviatoric part of the effective sub-grid
+        //  turbulence stress tensor including the laminar stress
+        virtual tmp<fvVectorMatrix> divDevRhoReff
+        (
+            const volScalarField& rho,
+            volVectorField& U
+        ) const;
 
         //- Solve the turbulence equations (k-w) and correct the turbulence
         //  viscosity
diff --git a/src/turbulenceModels/incompressible/LES/laminar/laminar.C b/src/turbulenceModels/incompressible/LES/laminar/laminar.C
index c88ff32beba..8b3ec0474e5 100644
--- a/src/turbulenceModels/incompressible/LES/laminar/laminar.C
+++ b/src/turbulenceModels/incompressible/LES/laminar/laminar.C
@@ -136,17 +136,34 @@ tmp<volSymmTensorField> laminar::B() const
 }
 
 
-tmp<volSymmTensorField> laminar::devBeff() const
+tmp<volSymmTensorField> laminar::devReff() const
 {
     return -nu()*dev(twoSymm(fvc::grad(U())));
 }
 
 
-tmp<fvVectorMatrix> laminar::divDevBeff(volVectorField& U) const
+tmp<fvVectorMatrix> laminar::divDevReff(volVectorField& U) const
 {
     return
     (
-      - fvm::laplacian(nu(), U) - fvc::div(nu()*dev(T(fvc::grad(U))))
+      - fvm::laplacian(nu(), U)
+      - fvc::div(nu()*dev(T(fvc::grad(U))))
+    );
+}
+
+
+tmp<fvVectorMatrix> laminar::divDevRhoReff
+(
+    const volScalarField& rho,
+    volVectorField& U
+) const
+{
+    volScalarField muEff("muEff", rho*nuEff());
+
+    return
+    (
+      - fvm::laplacian(muEff, U)
+      - fvc::div(muEff*dev(T(fvc::grad(U))))
     );
 }
 
diff --git a/src/turbulenceModels/incompressible/LES/laminar/laminar.H b/src/turbulenceModels/incompressible/LES/laminar/laminar.H
index 14de7a142ac..d7fa00c5557 100644
--- a/src/turbulenceModels/incompressible/LES/laminar/laminar.H
+++ b/src/turbulenceModels/incompressible/LES/laminar/laminar.H
@@ -104,13 +104,21 @@ public:
         //- Return the sub-grid stress tensor B.
         virtual tmp<volSymmTensorField> B() const;
 
+        //- Return the effective sub-grid turbulence stress tensor
+        //  including the laminar stress
+        virtual tmp<volSymmTensorField> devReff() const;
+
         //- Return the deviatoric part of the effective sub-grid
         //  turbulence stress tensor including the laminar stress
-        virtual tmp<volSymmTensorField> devBeff() const;
+        virtual tmp<fvVectorMatrix> divDevReff(volVectorField& U) const;
 
-        //- Return the deviatoric part of the divergence of Beff
-        //  i.e. the additional term in the filtered NSE.
-        virtual tmp<fvVectorMatrix> divDevBeff(volVectorField& U) const;
+        //- Return the deviatoric part of the effective sub-grid
+        //  turbulence stress tensor including the laminar stress
+        virtual tmp<fvVectorMatrix> divDevRhoReff
+        (
+            const volScalarField& rho,
+            volVectorField& U
+        ) const;
 
         //- Read LESProperties dictionary
         bool read();
diff --git a/src/turbulenceModels/incompressible/LES/mixedSmagorinsky/mixedSmagorinsky.C b/src/turbulenceModels/incompressible/LES/mixedSmagorinsky/mixedSmagorinsky.C
index 4f2804577fa..ed0e66e29cb 100644
--- a/src/turbulenceModels/incompressible/LES/mixedSmagorinsky/mixedSmagorinsky.C
+++ b/src/turbulenceModels/incompressible/LES/mixedSmagorinsky/mixedSmagorinsky.C
@@ -91,25 +91,39 @@ tmp<volSymmTensorField> mixedSmagorinsky::B() const
 }
 
 
-tmp<volSymmTensorField> mixedSmagorinsky::devBeff() const
+tmp<volSymmTensorField> mixedSmagorinsky::devReff() const
 {
     return
     (
-        scaleSimilarity::devBeff()
-      + Smagorinsky::devBeff()
+        scaleSimilarity::devReff()
+      + Smagorinsky::devReff()
     );
 }
 
 
-tmp<fvVectorMatrix> mixedSmagorinsky::divDevBeff
+tmp<fvVectorMatrix> mixedSmagorinsky::divDevReff
 (
     volVectorField& U
 ) const
 {
     return
     (
-        scaleSimilarity::divDevBeff(U)
-      + Smagorinsky::divDevBeff(U)
+        scaleSimilarity::divDevReff(U)
+      + Smagorinsky::divDevReff(U)
+    );
+}
+
+
+tmp<fvVectorMatrix> mixedSmagorinsky::divDevRhoReff
+(
+    const volScalarField& rho,
+    volVectorField& U
+) const
+{
+    return
+    (
+        scaleSimilarity::divDevRhoReff(rho, U)
+      + Smagorinsky::divDevRhoReff(rho, U)
     );
 }
 
diff --git a/src/turbulenceModels/incompressible/LES/mixedSmagorinsky/mixedSmagorinsky.H b/src/turbulenceModels/incompressible/LES/mixedSmagorinsky/mixedSmagorinsky.H
index c8557db2391..6d685f232ff 100644
--- a/src/turbulenceModels/incompressible/LES/mixedSmagorinsky/mixedSmagorinsky.H
+++ b/src/turbulenceModels/incompressible/LES/mixedSmagorinsky/mixedSmagorinsky.H
@@ -129,11 +129,19 @@ public:
 
         //- Return the effective sub-grid turbulence stress tensor
         //  including the laminar stress
-        virtual tmp<volSymmTensorField> devBeff() const;
+        virtual tmp<volSymmTensorField> devReff() const;
 
-        //- Implementation of div(B). This is necessary to override
-        // (and include) the div(B) terms from both the parent classes.
-        virtual tmp<fvVectorMatrix> divDevBeff(volVectorField& U) const;
+        //- Return the deviatoric part of the effective sub-grid
+        //  turbulence stress tensor including the laminar stress
+        virtual tmp<fvVectorMatrix> divDevReff(volVectorField& U) const;
+
+        //- Return the deviatoric part of the effective sub-grid
+        //  turbulence stress tensor including the laminar stress
+        virtual tmp<fvVectorMatrix> divDevRhoReff
+        (
+            const volScalarField& rho,
+            volVectorField& U
+        ) const;
 
         //- Correct Eddy-Viscosity and related properties
         virtual void correct(const tmp<volTensorField>& gradU);
diff --git a/src/turbulenceModels/incompressible/LES/scaleSimilarity/scaleSimilarity.C b/src/turbulenceModels/incompressible/LES/scaleSimilarity/scaleSimilarity.C
index b253573be13..c213f8432d6 100644
--- a/src/turbulenceModels/incompressible/LES/scaleSimilarity/scaleSimilarity.C
+++ b/src/turbulenceModels/incompressible/LES/scaleSimilarity/scaleSimilarity.C
@@ -79,15 +79,25 @@ tmp<volSymmTensorField> scaleSimilarity::B() const
 }
 
 
-tmp<volSymmTensorField> scaleSimilarity::devBeff() const
+tmp<volSymmTensorField> scaleSimilarity::devReff() const
 {
     return dev(B());
 }
 
 
-tmp<fvVectorMatrix> scaleSimilarity::divDevBeff(volVectorField& U) const
+tmp<fvVectorMatrix> scaleSimilarity::divDevReff(volVectorField& U) const
 {
-    return fvm::Su(fvc::div(devBeff()), U);
+    return fvm::Su(fvc::div(devReff()), U);
+}
+
+
+tmp<fvVectorMatrix> scaleSimilarity::divDevRhoReff
+(
+    const volScalarField& rho,
+    volVectorField& U
+) const
+{
+    return fvm::Su(fvc::div(rho*devReff()), U);
 }
 
 
diff --git a/src/turbulenceModels/incompressible/LES/scaleSimilarity/scaleSimilarity.H b/src/turbulenceModels/incompressible/LES/scaleSimilarity/scaleSimilarity.H
index 0dff31d87e3..4707d6b0cbf 100644
--- a/src/turbulenceModels/incompressible/LES/scaleSimilarity/scaleSimilarity.H
+++ b/src/turbulenceModels/incompressible/LES/scaleSimilarity/scaleSimilarity.H
@@ -106,13 +106,21 @@ public:
         //- Return the sub-grid stress tensor.
         virtual tmp<volSymmTensorField> B() const;
 
+        //- Return the effective sub-grid turbulence stress tensor
+        //  including the laminar stress
+        virtual tmp<volSymmTensorField> devReff() const;
+
         //- Return the deviatoric part of the effective sub-grid
         //  turbulence stress tensor including the laminar stress
-        virtual tmp<volSymmTensorField> devBeff() const;
+        virtual tmp<fvVectorMatrix> divDevReff(volVectorField& U) const;
 
-        //- Return the deviatoric part of the divergence of Beff
-        //  i.e. the additional term in the filtered NSE.
-        virtual tmp<fvVectorMatrix> divDevBeff(volVectorField& U) const;
+        //- Return the deviatoric part of the effective sub-grid
+        //  turbulence stress tensor including the laminar stress
+        virtual tmp<fvVectorMatrix> divDevRhoReff
+        (
+            const volScalarField& rho,
+            volVectorField& U
+        ) const;
 
         //- Correct Eddy-Viscosity and related properties
         virtual void correct(const tmp<volTensorField>&);
diff --git a/src/turbulenceModels/incompressible/RAS/LRR/LRR.C b/src/turbulenceModels/incompressible/RAS/LRR/LRR.C
index 30afa63ce99..b5a64a3f0c4 100644
--- a/src/turbulenceModels/incompressible/RAS/LRR/LRR.C
+++ b/src/turbulenceModels/incompressible/RAS/LRR/LRR.C
@@ -259,6 +259,44 @@ tmp<fvVectorMatrix> LRR::divDevReff(volVectorField& U) const
 }
 
 
+tmp<fvVectorMatrix> LRR::divDevRhoReff
+(
+    const volScalarField& rho,
+    volVectorField& U
+) const
+{
+    volScalarField muEff("muEff", rho*nuEff());
+
+    if (couplingFactor_.value() > 0.0)
+    {
+        return
+        (
+            fvc::div
+            (
+                rho*R_ + couplingFactor_*(rho*nut_)*fvc::grad(U),
+                "div((rho*R))"
+            )
+          + fvc::laplacian
+            (
+                (1.0 - couplingFactor_)*rho*nut_,
+                U,
+                "laplacian(muEff,U)"
+            )
+          - fvm::laplacian(muEff, U)
+        );
+    }
+    else
+    {
+        return
+        (
+            fvc::div(rho*R_)
+          + fvc::laplacian(rho*nut_, U, "laplacian(muEff,U)")
+          - fvm::laplacian(muEff, U)
+        );
+    }
+}
+
+
 bool LRR::read()
 {
     if (RASModel::read())
diff --git a/src/turbulenceModels/incompressible/RAS/LRR/LRR.H b/src/turbulenceModels/incompressible/RAS/LRR/LRR.H
index 2f1a119b9ca..dea191cfcee 100644
--- a/src/turbulenceModels/incompressible/RAS/LRR/LRR.H
+++ b/src/turbulenceModels/incompressible/RAS/LRR/LRR.H
@@ -173,6 +173,13 @@ public:
         //- Return the source term for the momentum equation
         virtual tmp<fvVectorMatrix> divDevReff(volVectorField& U) const;
 
+        //- Return the source term for the momentum equation
+        virtual tmp<fvVectorMatrix> divDevRhoReff
+        (
+            const volScalarField& rho,
+            volVectorField& U
+        ) const;
+
         //- Solve the turbulence equations and correct the turbulence viscosity
         virtual void correct();
 
diff --git a/src/turbulenceModels/incompressible/RAS/LamBremhorstKE/LamBremhorstKE.C b/src/turbulenceModels/incompressible/RAS/LamBremhorstKE/LamBremhorstKE.C
index 1671c625a08..7dab09b7400 100644
--- a/src/turbulenceModels/incompressible/RAS/LamBremhorstKE/LamBremhorstKE.C
+++ b/src/turbulenceModels/incompressible/RAS/LamBremhorstKE/LamBremhorstKE.C
@@ -204,6 +204,22 @@ tmp<fvVectorMatrix> LamBremhorstKE::divDevReff(volVectorField& U) const
 }
 
 
+tmp<fvVectorMatrix> LamBremhorstKE::divDevRhoReff
+(
+    const volScalarField& rho,
+    volVectorField& U
+) const
+{
+    volScalarField muEff("muEff", rho*nuEff());
+
+    return
+    (
+      - fvm::laplacian(muEff, U)
+      - fvc::div(muEff*dev(T(fvc::grad(U))))
+    );
+}
+
+
 bool LamBremhorstKE::read()
 {
     if (RASModel::read())
diff --git a/src/turbulenceModels/incompressible/RAS/LamBremhorstKE/LamBremhorstKE.H b/src/turbulenceModels/incompressible/RAS/LamBremhorstKE/LamBremhorstKE.H
index 449055f4c5d..ba618487126 100644
--- a/src/turbulenceModels/incompressible/RAS/LamBremhorstKE/LamBremhorstKE.H
+++ b/src/turbulenceModels/incompressible/RAS/LamBremhorstKE/LamBremhorstKE.H
@@ -147,6 +147,13 @@ public:
         //- Return the source term for the momentum equation
         virtual tmp<fvVectorMatrix> divDevReff(volVectorField& U) const;
 
+        //- Return the source term for the momentum equation
+        virtual tmp<fvVectorMatrix> divDevRhoReff
+        (
+            const volScalarField& rho,
+            volVectorField& U
+        ) const;
+
         //- Solve the turbulence equations and correct the turbulence viscosity
         virtual void correct();
 
diff --git a/src/turbulenceModels/incompressible/RAS/LaunderGibsonRSTM/LaunderGibsonRSTM.C b/src/turbulenceModels/incompressible/RAS/LaunderGibsonRSTM/LaunderGibsonRSTM.C
index aafa62af192..8fb4310bb31 100644
--- a/src/turbulenceModels/incompressible/RAS/LaunderGibsonRSTM/LaunderGibsonRSTM.C
+++ b/src/turbulenceModels/incompressible/RAS/LaunderGibsonRSTM/LaunderGibsonRSTM.C
@@ -276,7 +276,12 @@ tmp<fvVectorMatrix> LaunderGibsonRSTM::divDevReff(volVectorField& U) const
         return
         (
             fvc::div(R_ + couplingFactor_*nut_*fvc::grad(U), "div(R)")
-          + fvc::laplacian((1.0-couplingFactor_)*nut_, U, "laplacian(nuEff,U)")
+          + fvc::laplacian
+            (
+                (1.0 - couplingFactor_)*nut_,
+                U,
+                "laplacian(nuEff,U)"
+            )
           - fvm::laplacian(nuEff(), U)
         );
     }
@@ -292,6 +297,44 @@ tmp<fvVectorMatrix> LaunderGibsonRSTM::divDevReff(volVectorField& U) const
 }
 
 
+tmp<fvVectorMatrix> LaunderGibsonRSTM::divDevRhoReff
+(
+    const volScalarField& rho,
+    volVectorField& U
+) const
+{
+    volScalarField muEff("muEff", rho*nuEff());
+
+    if (couplingFactor_.value() > 0.0)
+    {
+        return
+        (
+            fvc::div
+            (
+                rho*R_ + couplingFactor_*(rho*nut_)*fvc::grad(U),
+                "div((rho*R))"
+            )
+          + fvc::laplacian
+            (
+                (1.0 - couplingFactor_)*rho*nut_,
+                U,
+                "laplacian(muEff,U)"
+            )
+          - fvm::laplacian(muEff, U)
+        );
+    }
+    else
+    {
+        return
+        (
+            fvc::div(rho*R_)
+          + fvc::laplacian(rho*nut_, U, "laplacian(muEff,U)")
+          - fvm::laplacian(muEff, U)
+        );
+    }
+}
+
+
 bool LaunderGibsonRSTM::read()
 {
     if (RASModel::read())
diff --git a/src/turbulenceModels/incompressible/RAS/LaunderGibsonRSTM/LaunderGibsonRSTM.H b/src/turbulenceModels/incompressible/RAS/LaunderGibsonRSTM/LaunderGibsonRSTM.H
index 0b6fc91b4bb..c6f325d8548 100644
--- a/src/turbulenceModels/incompressible/RAS/LaunderGibsonRSTM/LaunderGibsonRSTM.H
+++ b/src/turbulenceModels/incompressible/RAS/LaunderGibsonRSTM/LaunderGibsonRSTM.H
@@ -184,6 +184,13 @@ public:
         //- Return the source term for the momentum equation
         virtual tmp<fvVectorMatrix> divDevReff(volVectorField& U) const;
 
+        //- Return the source term for the momentum equation
+        virtual tmp<fvVectorMatrix> divDevRhoReff
+        (
+            const volScalarField& rho,
+            volVectorField& U
+        ) const;
+
         //- Solve the turbulence equations and correct the turbulence viscosity
         virtual void correct();
 
diff --git a/src/turbulenceModels/incompressible/RAS/LaunderSharmaKE/LaunderSharmaKE.C b/src/turbulenceModels/incompressible/RAS/LaunderSharmaKE/LaunderSharmaKE.C
index f0fa6e60c8d..e20299c29f8 100644
--- a/src/turbulenceModels/incompressible/RAS/LaunderSharmaKE/LaunderSharmaKE.C
+++ b/src/turbulenceModels/incompressible/RAS/LaunderSharmaKE/LaunderSharmaKE.C
@@ -210,6 +210,22 @@ tmp<fvVectorMatrix> LaunderSharmaKE::divDevReff(volVectorField& U) const
 }
 
 
+tmp<fvVectorMatrix> LaunderSharmaKE::divDevRhoReff
+(
+    const volScalarField& rho,
+    volVectorField& U
+) const
+{
+    volScalarField muEff("muEff", rho*nuEff());
+
+    return
+    (
+      - fvm::laplacian(muEff, U)
+      - fvc::div(muEff*dev(T(fvc::grad(U))))
+    );
+}
+
+
 bool LaunderSharmaKE::read()
 {
     if (RASModel::read())
diff --git a/src/turbulenceModels/incompressible/RAS/LaunderSharmaKE/LaunderSharmaKE.H b/src/turbulenceModels/incompressible/RAS/LaunderSharmaKE/LaunderSharmaKE.H
index 3f2338e18e8..625de1dcb8b 100644
--- a/src/turbulenceModels/incompressible/RAS/LaunderSharmaKE/LaunderSharmaKE.H
+++ b/src/turbulenceModels/incompressible/RAS/LaunderSharmaKE/LaunderSharmaKE.H
@@ -164,6 +164,13 @@ public:
         //- Return the source term for the momentum equation
         virtual tmp<fvVectorMatrix> divDevReff(volVectorField& U) const;
 
+        //- Return the source term for the momentum equation
+        virtual tmp<fvVectorMatrix> divDevRhoReff
+        (
+            const volScalarField& rho,
+            volVectorField& U
+        ) const;
+
         //- Solve the turbulence equations and correct the turbulence viscosity
         virtual void correct();
 
diff --git a/src/turbulenceModels/incompressible/RAS/LienCubicKE/LienCubicKE.C b/src/turbulenceModels/incompressible/RAS/LienCubicKE/LienCubicKE.C
index 03e1e5aa1ac..20d8dee9514 100644
--- a/src/turbulenceModels/incompressible/RAS/LienCubicKE/LienCubicKE.C
+++ b/src/turbulenceModels/incompressible/RAS/LienCubicKE/LienCubicKE.C
@@ -298,6 +298,23 @@ tmp<fvVectorMatrix> LienCubicKE::divDevReff(volVectorField& U) const
 }
 
 
+tmp<fvVectorMatrix> LienCubicKE::divDevRhoReff
+(
+    const volScalarField& rho,
+    volVectorField& U
+) const
+{
+    volScalarField muEff("muEff", rho*nuEff());
+
+    return
+    (
+        fvc::div(rho*nonlinearStress_)
+      - fvm::laplacian(muEff, U)
+      - fvc::div(muEff*dev(T(fvc::grad(U))))
+    );
+}
+
+
 bool LienCubicKE::read()
 {
     if (RASModel::read())
diff --git a/src/turbulenceModels/incompressible/RAS/LienCubicKE/LienCubicKE.H b/src/turbulenceModels/incompressible/RAS/LienCubicKE/LienCubicKE.H
index ebe7c373d8e..56696f42539 100644
--- a/src/turbulenceModels/incompressible/RAS/LienCubicKE/LienCubicKE.H
+++ b/src/turbulenceModels/incompressible/RAS/LienCubicKE/LienCubicKE.H
@@ -159,6 +159,13 @@ public:
         //- Return the source term for the momentum equation
         virtual tmp<fvVectorMatrix> divDevReff(volVectorField& U) const;
 
+        //- Return the source term for the momentum equation
+        virtual tmp<fvVectorMatrix> divDevRhoReff
+        (
+            const volScalarField& rho,
+            volVectorField& U
+        ) const;
+
         //- Solve the turbulence equations and correct the turbulence viscosity
         virtual void correct();
 
diff --git a/src/turbulenceModels/incompressible/RAS/LienCubicKELowRe/LienCubicKELowRe.C b/src/turbulenceModels/incompressible/RAS/LienCubicKELowRe/LienCubicKELowRe.C
index e0bd0f0df9f..da81b66eddd 100644
--- a/src/turbulenceModels/incompressible/RAS/LienCubicKELowRe/LienCubicKELowRe.C
+++ b/src/turbulenceModels/incompressible/RAS/LienCubicKELowRe/LienCubicKELowRe.C
@@ -365,6 +365,23 @@ tmp<fvVectorMatrix> LienCubicKELowRe::divDevReff(volVectorField& U) const
 }
 
 
+tmp<fvVectorMatrix> LienCubicKELowRe::divDevRhoReff
+(
+    const volScalarField& rho,
+    volVectorField& U
+) const
+{
+    volScalarField muEff("muEff", rho*nuEff());
+
+    return
+    (
+        fvc::div(rho*nonlinearStress_)
+      - fvm::laplacian(muEff, U)
+      - fvc::div(muEff*dev(T(fvc::grad(U))))
+    );
+}
+
+
 bool LienCubicKELowRe::read()
 {
     if (RASModel::read())
diff --git a/src/turbulenceModels/incompressible/RAS/LienCubicKELowRe/LienCubicKELowRe.H b/src/turbulenceModels/incompressible/RAS/LienCubicKELowRe/LienCubicKELowRe.H
index 0efd63a454c..8794af1fc4b 100644
--- a/src/turbulenceModels/incompressible/RAS/LienCubicKELowRe/LienCubicKELowRe.H
+++ b/src/turbulenceModels/incompressible/RAS/LienCubicKELowRe/LienCubicKELowRe.H
@@ -186,6 +186,13 @@ public:
         //- Return the source term for the momentum equation
         virtual tmp<fvVectorMatrix> divDevReff(volVectorField& U) const;
 
+        //- Return the source term for the momentum equation
+        virtual tmp<fvVectorMatrix> divDevRhoReff
+        (
+            const volScalarField& rho,
+            volVectorField& U
+        ) const;
+
         //- Solve the turbulence equations and correct the turbulence viscosity
         virtual void correct();
 
diff --git a/src/turbulenceModels/incompressible/RAS/LienLeschzinerLowRe/LienLeschzinerLowRe.C b/src/turbulenceModels/incompressible/RAS/LienLeschzinerLowRe/LienLeschzinerLowRe.C
index 91e288f1042..de4ce7be751 100644
--- a/src/turbulenceModels/incompressible/RAS/LienLeschzinerLowRe/LienLeschzinerLowRe.C
+++ b/src/turbulenceModels/incompressible/RAS/LienLeschzinerLowRe/LienLeschzinerLowRe.C
@@ -243,12 +243,27 @@ tmp<fvVectorMatrix> LienLeschzinerLowRe::divDevReff(volVectorField& U) const
     return
     (
       - fvm::laplacian(nuEff(), U)
-    //- (fvc::grad(U) & fvc::grad(nuEff()))
       - fvc::div(nuEff()*T(fvc::grad(U)))
     );
 }
 
 
+tmp<fvVectorMatrix> LienLeschzinerLowRe::divDevRhoReff
+(
+    const volScalarField& rho,
+    volVectorField& U
+) const
+{
+    volScalarField muEff("muEff", rho*nuEff());
+
+    return
+    (
+      - fvm::laplacian(muEff, U)
+      - fvc::div(muEff*dev(T(fvc::grad(U))))
+    );
+}
+
+
 bool LienLeschzinerLowRe::read()
 {
     if (RASModel::read())
diff --git a/src/turbulenceModels/incompressible/RAS/LienLeschzinerLowRe/LienLeschzinerLowRe.H b/src/turbulenceModels/incompressible/RAS/LienLeschzinerLowRe/LienLeschzinerLowRe.H
index 3fa1c790fe0..b99fffdcdf6 100644
--- a/src/turbulenceModels/incompressible/RAS/LienLeschzinerLowRe/LienLeschzinerLowRe.H
+++ b/src/turbulenceModels/incompressible/RAS/LienLeschzinerLowRe/LienLeschzinerLowRe.H
@@ -156,6 +156,13 @@ public:
         //- Return the source term for the momentum equation
         virtual tmp<fvVectorMatrix> divDevReff(volVectorField& U) const;
 
+        //- Return the source term for the momentum equation
+        virtual tmp<fvVectorMatrix> divDevRhoReff
+        (
+            const volScalarField& rho,
+            volVectorField& U
+        ) const;
+
         //- Solve the turbulence equations and correct the turbulence viscosity
         virtual void correct();
 
diff --git a/src/turbulenceModels/incompressible/RAS/NonlinearKEShih/NonlinearKEShih.C b/src/turbulenceModels/incompressible/RAS/NonlinearKEShih/NonlinearKEShih.C
index 33cc03edac8..4dbb7b7b683 100644
--- a/src/turbulenceModels/incompressible/RAS/NonlinearKEShih/NonlinearKEShih.C
+++ b/src/turbulenceModels/incompressible/RAS/NonlinearKEShih/NonlinearKEShih.C
@@ -287,6 +287,23 @@ tmp<fvVectorMatrix> NonlinearKEShih::divDevReff(volVectorField& U) const
 }
 
 
+tmp<fvVectorMatrix> NonlinearKEShih::divDevRhoReff
+(
+    const volScalarField& rho,
+    volVectorField& U
+) const
+{
+    volScalarField muEff("muEff", rho*nuEff());
+
+    return
+    (
+        fvc::div(rho*nonlinearStress_)
+      - fvm::laplacian(muEff, U)
+      - fvc::div(muEff*dev(T(fvc::grad(U))))
+    );
+}
+
+
 bool NonlinearKEShih::read()
 {
     if (RASModel::read())
diff --git a/src/turbulenceModels/incompressible/RAS/NonlinearKEShih/NonlinearKEShih.H b/src/turbulenceModels/incompressible/RAS/NonlinearKEShih/NonlinearKEShih.H
index c7ad8a37131..f63f0c6565c 100644
--- a/src/turbulenceModels/incompressible/RAS/NonlinearKEShih/NonlinearKEShih.H
+++ b/src/turbulenceModels/incompressible/RAS/NonlinearKEShih/NonlinearKEShih.H
@@ -162,6 +162,13 @@ public:
         //- Return the source term for the momentum equation
         virtual tmp<fvVectorMatrix> divDevReff(volVectorField& U) const;
 
+        //- Return the source term for the momentum equation
+        virtual tmp<fvVectorMatrix> divDevRhoReff
+        (
+            const volScalarField& rho,
+            volVectorField& U
+        ) const;
+
         //- Solve the turbulence equations and correct the turbulence viscosity
         virtual void correct();
 
diff --git a/src/turbulenceModels/incompressible/RAS/RNGkEpsilon/RNGkEpsilon.C b/src/turbulenceModels/incompressible/RAS/RNGkEpsilon/RNGkEpsilon.C
index 794d9f9b07f..99c7cc4c773 100644
--- a/src/turbulenceModels/incompressible/RAS/RNGkEpsilon/RNGkEpsilon.C
+++ b/src/turbulenceModels/incompressible/RAS/RNGkEpsilon/RNGkEpsilon.C
@@ -220,6 +220,22 @@ tmp<fvVectorMatrix> RNGkEpsilon::divDevReff(volVectorField& U) const
 }
 
 
+tmp<fvVectorMatrix> RNGkEpsilon::divDevRhoReff
+(
+    const volScalarField& rho,
+    volVectorField& U
+) const
+{
+    volScalarField muEff("muEff", rho*nuEff());
+
+    return
+    (
+      - fvm::laplacian(muEff, U)
+      - fvc::div(muEff*dev(T(fvc::grad(U))))
+    );
+}
+
+
 bool RNGkEpsilon::read()
 {
     if (RASModel::read())
diff --git a/src/turbulenceModels/incompressible/RAS/RNGkEpsilon/RNGkEpsilon.H b/src/turbulenceModels/incompressible/RAS/RNGkEpsilon/RNGkEpsilon.H
index 62b60b98046..077acc70f88 100644
--- a/src/turbulenceModels/incompressible/RAS/RNGkEpsilon/RNGkEpsilon.H
+++ b/src/turbulenceModels/incompressible/RAS/RNGkEpsilon/RNGkEpsilon.H
@@ -161,6 +161,13 @@ public:
         //- Return the source term for the momentum equation
         virtual tmp<fvVectorMatrix> divDevReff(volVectorField& U) const;
 
+        //- Return the source term for the momentum equation
+        virtual tmp<fvVectorMatrix> divDevRhoReff
+        (
+            const volScalarField& rho,
+            volVectorField& U
+        ) const;
+
         //- Solve the turbulence equations and correct the turbulence viscosity
         virtual void correct();
 
diff --git a/src/turbulenceModels/incompressible/RAS/SpalartAllmaras/SpalartAllmaras.C b/src/turbulenceModels/incompressible/RAS/SpalartAllmaras/SpalartAllmaras.C
index 1804ba5b010..4105cb1876e 100644
--- a/src/turbulenceModels/incompressible/RAS/SpalartAllmaras/SpalartAllmaras.C
+++ b/src/turbulenceModels/incompressible/RAS/SpalartAllmaras/SpalartAllmaras.C
@@ -337,6 +337,22 @@ tmp<fvVectorMatrix> SpalartAllmaras::divDevReff(volVectorField& U) const
 }
 
 
+tmp<fvVectorMatrix> SpalartAllmaras::divDevRhoReff
+(
+    const volScalarField& rho,
+    volVectorField& U
+) const
+{
+    volScalarField muEff("muEff", rho*nuEff());
+
+    return
+    (
+      - fvm::laplacian(muEff, U)
+      - fvc::div(muEff*dev(T(fvc::grad(U))))
+    );
+}
+
+
 bool SpalartAllmaras::read()
 {
     if (RASModel::read())
diff --git a/src/turbulenceModels/incompressible/RAS/SpalartAllmaras/SpalartAllmaras.H b/src/turbulenceModels/incompressible/RAS/SpalartAllmaras/SpalartAllmaras.H
index fbd12481810..2e2626534d2 100644
--- a/src/turbulenceModels/incompressible/RAS/SpalartAllmaras/SpalartAllmaras.H
+++ b/src/turbulenceModels/incompressible/RAS/SpalartAllmaras/SpalartAllmaras.H
@@ -184,6 +184,13 @@ public:
         //- Return the source term for the momentum equation
         virtual tmp<fvVectorMatrix> divDevReff(volVectorField& U) const;
 
+        //- Return the source term for the momentum equation
+        virtual tmp<fvVectorMatrix> divDevRhoReff
+        (
+            const volScalarField& rho,
+            volVectorField& U
+        ) const;
+
         //- Solve the turbulence equations and correct the turbulence viscosity
         virtual void correct();
 
diff --git a/src/turbulenceModels/incompressible/RAS/kEpsilon/kEpsilon.C b/src/turbulenceModels/incompressible/RAS/kEpsilon/kEpsilon.C
index 99b4188cc01..6560ff4db6f 100644
--- a/src/turbulenceModels/incompressible/RAS/kEpsilon/kEpsilon.C
+++ b/src/turbulenceModels/incompressible/RAS/kEpsilon/kEpsilon.C
@@ -192,6 +192,22 @@ tmp<fvVectorMatrix> kEpsilon::divDevReff(volVectorField& U) const
 }
 
 
+tmp<fvVectorMatrix> kEpsilon::divDevRhoReff
+(
+    const volScalarField& rho,
+    volVectorField& U
+) const
+{
+    volScalarField muEff("muEff", rho*nuEff());
+
+    return
+    (
+      - fvm::laplacian(muEff, U)
+      - fvc::div(muEff*dev(T(fvc::grad(U))))
+    );
+}
+
+
 bool kEpsilon::read()
 {
     if (RASModel::read())
diff --git a/src/turbulenceModels/incompressible/RAS/kEpsilon/kEpsilon.H b/src/turbulenceModels/incompressible/RAS/kEpsilon/kEpsilon.H
index a3baf205134..df9b15d4999 100644
--- a/src/turbulenceModels/incompressible/RAS/kEpsilon/kEpsilon.H
+++ b/src/turbulenceModels/incompressible/RAS/kEpsilon/kEpsilon.H
@@ -155,6 +155,13 @@ public:
         //- Return the source term for the momentum equation
         virtual tmp<fvVectorMatrix> divDevReff(volVectorField& U) const;
 
+        //- Return the source term for the momentum equation
+        virtual tmp<fvVectorMatrix> divDevRhoReff
+        (
+            const volScalarField& rho,
+            volVectorField& U
+        ) const;
+
         //- Solve the turbulence equations and correct the turbulence viscosity
         virtual void correct();
 
diff --git a/src/turbulenceModels/incompressible/RAS/kOmega/kOmega.C b/src/turbulenceModels/incompressible/RAS/kOmega/kOmega.C
index 7106c0ea2cd..05eb253f5fd 100644
--- a/src/turbulenceModels/incompressible/RAS/kOmega/kOmega.C
+++ b/src/turbulenceModels/incompressible/RAS/kOmega/kOmega.C
@@ -201,6 +201,22 @@ tmp<fvVectorMatrix> kOmega::divDevReff(volVectorField& U) const
 }
 
 
+tmp<fvVectorMatrix> kOmega::divDevRhoReff
+(
+    const volScalarField& rho,
+    volVectorField& U
+) const
+{
+    volScalarField muEff("muEff", rho*nuEff());
+
+    return
+    (
+      - fvm::laplacian(muEff, U)
+      - fvc::div(muEff*dev(T(fvc::grad(U))))
+    );
+}
+
+
 bool kOmega::read()
 {
     if (RASModel::read())
diff --git a/src/turbulenceModels/incompressible/RAS/kOmega/kOmega.H b/src/turbulenceModels/incompressible/RAS/kOmega/kOmega.H
index 10729f334e7..712aa1a193c 100644
--- a/src/turbulenceModels/incompressible/RAS/kOmega/kOmega.H
+++ b/src/turbulenceModels/incompressible/RAS/kOmega/kOmega.H
@@ -188,6 +188,13 @@ public:
         //- Return the source term for the momentum equation
         virtual tmp<fvVectorMatrix> divDevReff(volVectorField& U) const;
 
+        //- Return the source term for the momentum equation
+        virtual tmp<fvVectorMatrix> divDevRhoReff
+        (
+            const volScalarField& rho,
+            volVectorField& U
+        ) const;
+
         //- Solve the turbulence equations and correct the turbulence viscosity
         virtual void correct();
 
diff --git a/src/turbulenceModels/incompressible/RAS/kOmegaSST/kOmegaSST.C b/src/turbulenceModels/incompressible/RAS/kOmegaSST/kOmegaSST.C
index c025076baa9..3481bc5dbaf 100644
--- a/src/turbulenceModels/incompressible/RAS/kOmegaSST/kOmegaSST.C
+++ b/src/turbulenceModels/incompressible/RAS/kOmegaSST/kOmegaSST.C
@@ -308,6 +308,22 @@ tmp<fvVectorMatrix> kOmegaSST::divDevReff(volVectorField& U) const
 }
 
 
+tmp<fvVectorMatrix> kOmegaSST::divDevRhoReff
+(
+    const volScalarField& rho,
+    volVectorField& U
+) const
+{
+    volScalarField muEff("muEff", rho*nuEff());
+
+    return
+    (
+      - fvm::laplacian(muEff, U)
+      - fvc::div(muEff*dev(T(fvc::grad(U))))
+    );
+}
+
+
 bool kOmegaSST::read()
 {
     if (RASModel::read())
diff --git a/src/turbulenceModels/incompressible/RAS/kOmegaSST/kOmegaSST.H b/src/turbulenceModels/incompressible/RAS/kOmegaSST/kOmegaSST.H
index bd4b3742b5e..94438d6825b 100644
--- a/src/turbulenceModels/incompressible/RAS/kOmegaSST/kOmegaSST.H
+++ b/src/turbulenceModels/incompressible/RAS/kOmegaSST/kOmegaSST.H
@@ -258,6 +258,13 @@ public:
         //- Return the source term for the momentum equation
         virtual tmp<fvVectorMatrix> divDevReff(volVectorField& U) const;
 
+        //- Return the source term for the momentum equation
+        virtual tmp<fvVectorMatrix> divDevRhoReff
+        (
+            const volScalarField& rho,
+            volVectorField& U
+        ) const;
+
         //- Solve the turbulence equations and correct the turbulence viscosity
         virtual void correct();
 
diff --git a/src/turbulenceModels/incompressible/RAS/kkLOmega/kkLOmega.C b/src/turbulenceModels/incompressible/RAS/kkLOmega/kkLOmega.C
index abf9b29cde7..45f4f69d130 100644
--- a/src/turbulenceModels/incompressible/RAS/kkLOmega/kkLOmega.C
+++ b/src/turbulenceModels/incompressible/RAS/kkLOmega/kkLOmega.C
@@ -570,6 +570,22 @@ tmp<fvVectorMatrix> kkLOmega::divDevReff(volVectorField& U) const
 }
 
 
+tmp<fvVectorMatrix> kkLOmega::divDevRhoReff
+(
+    const volScalarField& rho,
+    volVectorField& U
+) const
+{
+    volScalarField muEff("muEff", rho*nuEff());
+
+    return
+    (
+      - fvm::laplacian(muEff, U)
+      - fvc::div(muEff*dev(T(fvc::grad(U))))
+    );
+}
+
+
 bool kkLOmega::read()
 {
     if (RASModel::read())
diff --git a/src/turbulenceModels/incompressible/RAS/kkLOmega/kkLOmega.H b/src/turbulenceModels/incompressible/RAS/kkLOmega/kkLOmega.H
index 83db5f82d5c..51ddfb7f04d 100644
--- a/src/turbulenceModels/incompressible/RAS/kkLOmega/kkLOmega.H
+++ b/src/turbulenceModels/incompressible/RAS/kkLOmega/kkLOmega.H
@@ -276,6 +276,13 @@ public:
         //- Return the source term for the momentum equation
         virtual tmp<fvVectorMatrix> divDevReff(volVectorField& U) const;
 
+        //- Return the source term for the momentum equation
+        virtual tmp<fvVectorMatrix> divDevRhoReff
+        (
+            const volScalarField& rho,
+            volVectorField& U
+        ) const;
+
         //- Solve the turbulence equations and correct the turbulence viscosity
         virtual void correct();
 
diff --git a/src/turbulenceModels/incompressible/RAS/laminar/laminar.C b/src/turbulenceModels/incompressible/RAS/laminar/laminar.C
index 68adca77b3f..3c4560c4dc3 100644
--- a/src/turbulenceModels/incompressible/RAS/laminar/laminar.C
+++ b/src/turbulenceModels/incompressible/RAS/laminar/laminar.C
@@ -177,6 +177,22 @@ tmp<fvVectorMatrix> laminar::divDevReff(volVectorField& U) const
 }
 
 
+tmp<fvVectorMatrix> laminar::divDevRhoReff
+(
+    const volScalarField& rho,
+    volVectorField& U
+) const
+{
+    volScalarField muEff("muEff", rho*nuEff());
+
+    return
+    (
+      - fvm::laplacian(muEff, U)
+      - fvc::div(muEff*dev(T(fvc::grad(U))))
+    );
+}
+
+
 bool laminar::read()
 {
     return RASModel::read();
diff --git a/src/turbulenceModels/incompressible/RAS/laminar/laminar.H b/src/turbulenceModels/incompressible/RAS/laminar/laminar.H
index 90afadff132..ce51d01c003 100644
--- a/src/turbulenceModels/incompressible/RAS/laminar/laminar.H
+++ b/src/turbulenceModels/incompressible/RAS/laminar/laminar.H
@@ -105,6 +105,13 @@ public:
         //- Return the source term for the momentum equation
         virtual tmp<fvVectorMatrix> divDevReff(volVectorField& U) const;
 
+        //- Return the source term for the momentum equation
+        virtual tmp<fvVectorMatrix> divDevRhoReff
+        (
+            const volScalarField& rho,
+            volVectorField& U
+        ) const;
+
         //- Correct the laminar viscosity
         virtual void correct();
 
diff --git a/src/turbulenceModels/incompressible/RAS/qZeta/qZeta.C b/src/turbulenceModels/incompressible/RAS/qZeta/qZeta.C
index 25940abb9ef..722bff88b18 100644
--- a/src/turbulenceModels/incompressible/RAS/qZeta/qZeta.C
+++ b/src/turbulenceModels/incompressible/RAS/qZeta/qZeta.C
@@ -262,6 +262,22 @@ tmp<fvVectorMatrix> qZeta::divDevReff(volVectorField& U) const
 }
 
 
+tmp<fvVectorMatrix> qZeta::divDevRhoReff
+(
+    const volScalarField& rho,
+    volVectorField& U
+) const
+{
+    volScalarField muEff("muEff", rho*nuEff());
+
+    return
+    (
+      - fvm::laplacian(muEff, U)
+      - fvc::div(muEff*dev(T(fvc::grad(U))))
+    );
+}
+
+
 bool qZeta::read()
 {
     if (RASModel::read())
diff --git a/src/turbulenceModels/incompressible/RAS/qZeta/qZeta.H b/src/turbulenceModels/incompressible/RAS/qZeta/qZeta.H
index ba8dc1b413f..c9c272975da 100644
--- a/src/turbulenceModels/incompressible/RAS/qZeta/qZeta.H
+++ b/src/turbulenceModels/incompressible/RAS/qZeta/qZeta.H
@@ -211,6 +211,13 @@ public:
         //- Return the source term for the momentum equation
         virtual tmp<fvVectorMatrix> divDevReff(volVectorField& U) const;
 
+        //- Return the source term for the momentum equation
+        virtual tmp<fvVectorMatrix> divDevRhoReff
+        (
+            const volScalarField& rho,
+            volVectorField& U
+        ) const;
+
         //- Solve the turbulence equations and correct the turbulence viscosity
         virtual void correct();
 
diff --git a/src/turbulenceModels/incompressible/RAS/realizableKE/realizableKE.C b/src/turbulenceModels/incompressible/RAS/realizableKE/realizableKE.C
index 71058caa7d9..2595b3f294b 100644
--- a/src/turbulenceModels/incompressible/RAS/realizableKE/realizableKE.C
+++ b/src/turbulenceModels/incompressible/RAS/realizableKE/realizableKE.C
@@ -246,6 +246,22 @@ tmp<fvVectorMatrix> realizableKE::divDevReff(volVectorField& U) const
 }
 
 
+tmp<fvVectorMatrix> realizableKE::divDevRhoReff
+(
+    const volScalarField& rho,
+    volVectorField& U
+) const
+{
+    volScalarField muEff("muEff", rho*nuEff());
+
+    return
+    (
+      - fvm::laplacian(muEff, U)
+      - fvc::div(muEff*dev(T(fvc::grad(U))))
+    );
+}
+
+
 bool realizableKE::read()
 {
     if (RASModel::read())
diff --git a/src/turbulenceModels/incompressible/RAS/realizableKE/realizableKE.H b/src/turbulenceModels/incompressible/RAS/realizableKE/realizableKE.H
index 045f3333908..147febcb508 100644
--- a/src/turbulenceModels/incompressible/RAS/realizableKE/realizableKE.H
+++ b/src/turbulenceModels/incompressible/RAS/realizableKE/realizableKE.H
@@ -180,6 +180,13 @@ public:
         //- Return the source term for the momentum equation
         virtual tmp<fvVectorMatrix> divDevReff(volVectorField& U) const;
 
+        //- Return the source term for the momentum equation
+        virtual tmp<fvVectorMatrix> divDevRhoReff
+        (
+            const volScalarField& rho,
+            volVectorField& U
+        ) const;
+
         //- Solve the turbulence equations and correct the turbulence viscosity
         virtual void correct();
 
diff --git a/src/turbulenceModels/incompressible/turbulenceModel/laminar/laminar.C b/src/turbulenceModels/incompressible/turbulenceModel/laminar/laminar.C
index c102b52cb1e..9d1e77d5abe 100644
--- a/src/turbulenceModels/incompressible/turbulenceModel/laminar/laminar.C
+++ b/src/turbulenceModels/incompressible/turbulenceModel/laminar/laminar.C
@@ -203,6 +203,22 @@ tmp<fvVectorMatrix> laminar::divDevReff(volVectorField& U) const
 }
 
 
+tmp<fvVectorMatrix> laminar::divDevRhoReff
+(
+    const volScalarField& rho,
+    volVectorField& U
+) const
+{
+    volScalarField muEff("muEff", rho*nuEff());
+
+    return
+    (
+      - fvm::laplacian(muEff, U)
+      - fvc::div(muEff*dev(T(fvc::grad(U))))
+    );
+}
+
+
 void laminar::correct()
 {
     turbulenceModel::correct();
diff --git a/src/turbulenceModels/incompressible/turbulenceModel/laminar/laminar.H b/src/turbulenceModels/incompressible/turbulenceModel/laminar/laminar.H
index 72994e90fb8..09cdf4faa7e 100644
--- a/src/turbulenceModels/incompressible/turbulenceModel/laminar/laminar.H
+++ b/src/turbulenceModels/incompressible/turbulenceModel/laminar/laminar.H
@@ -111,6 +111,13 @@ public:
         //- Return the source term for the momentum equation
         virtual tmp<fvVectorMatrix> divDevReff(volVectorField& U) const;
 
+        //- Return the source term for the momentum equation
+        virtual tmp<fvVectorMatrix> divDevRhoReff
+        (
+            const volScalarField& rho,
+            volVectorField& U
+        ) const;
+
         //- Correct the laminar viscosity
         virtual void correct();
 
diff --git a/src/turbulenceModels/incompressible/turbulenceModel/turbulenceModel.H b/src/turbulenceModels/incompressible/turbulenceModel/turbulenceModel.H
index 405cfef0bfa..c2eacdaf392 100644
--- a/src/turbulenceModels/incompressible/turbulenceModel/turbulenceModel.H
+++ b/src/turbulenceModels/incompressible/turbulenceModel/turbulenceModel.H
@@ -204,6 +204,13 @@ public:
         //- Return the source term for the momentum equation
         virtual tmp<fvVectorMatrix> divDevReff(volVectorField& U) const = 0;
 
+        //- Return the source term for the momentum equation
+        virtual tmp<fvVectorMatrix> divDevRhoReff
+        (
+            const volScalarField& rho,
+            volVectorField& U
+        ) const = 0;
+
         //- Solve the turbulence equations and correct the turbulence viscosity
         virtual void correct() = 0;
 
diff --git a/tutorials/multiphase/LTSInterFoam/wigleyHull/system/fvSchemes b/tutorials/multiphase/LTSInterFoam/wigleyHull/system/fvSchemes
index 22366f437b7..7018c5eaa73 100644
--- a/tutorials/multiphase/LTSInterFoam/wigleyHull/system/fvSchemes
+++ b/tutorials/multiphase/LTSInterFoam/wigleyHull/system/fvSchemes
@@ -32,6 +32,7 @@ divSchemes
     div(phirb,alpha) Gauss interfaceCompression;
     div(phi,k)      Gauss upwind;
     div(phi,omega)  Gauss upwind;
+    div((muEff*dev(T(grad(U))))) Gauss linear;
 }
 
 laplacianSchemes
diff --git a/tutorials/multiphase/MRFInterFoam/mixerVessel2D/system/fvSchemes b/tutorials/multiphase/MRFInterFoam/mixerVessel2D/system/fvSchemes
index 41159c5be64..3827ac6548c 100644
--- a/tutorials/multiphase/MRFInterFoam/mixerVessel2D/system/fvSchemes
+++ b/tutorials/multiphase/MRFInterFoam/mixerVessel2D/system/fvSchemes
@@ -30,6 +30,7 @@ divSchemes
     div(rho*phi,U)  Gauss linear;
     div(phi,alpha)  Gauss vanLeer;
     div(phirb,alpha) Gauss interfaceCompression;
+    div((muEff*dev(T(grad(U))))) Gauss linear;
 }
 
 laplacianSchemes
diff --git a/tutorials/multiphase/MRFMultiphaseInterFoam/mixerVessel2D/system/fvSchemes b/tutorials/multiphase/MRFMultiphaseInterFoam/mixerVessel2D/system/fvSchemes
index b4ad45237d1..d1ed26cae3b 100644
--- a/tutorials/multiphase/MRFMultiphaseInterFoam/mixerVessel2D/system/fvSchemes
+++ b/tutorials/multiphase/MRFMultiphaseInterFoam/mixerVessel2D/system/fvSchemes
@@ -30,6 +30,7 @@ divSchemes
     div(rho*phi,U)  Gauss limitedLinearV 1;
     div(phi,alpha)  Gauss vanLeer;
     div(phirb,alpha) Gauss interfaceCompression;
+    div((muEff*dev(T(grad(U))))) Gauss linear;
 }
 
 laplacianSchemes
diff --git a/tutorials/multiphase/cavitatingFoam/les/throttle/system/fvSchemes b/tutorials/multiphase/cavitatingFoam/les/throttle/system/fvSchemes
index 0c8bf54bb51..3d506aa1d7f 100644
--- a/tutorials/multiphase/cavitatingFoam/les/throttle/system/fvSchemes
+++ b/tutorials/multiphase/cavitatingFoam/les/throttle/system/fvSchemes
@@ -31,6 +31,7 @@ divSchemes
     div(phiv,rho)   Gauss limitedLinear 0.2;
     div(phi,U)      Gauss filteredLinear2V 0.2 0;
     div(phiv,k)     Gauss filteredLinear2 0.2 0;
+    div((muEff*dev(T(grad(U))))) Gauss linear;
 }
 
 gradSchemes
diff --git a/tutorials/multiphase/cavitatingFoam/les/throttle3D/system/fvSchemes b/tutorials/multiphase/cavitatingFoam/les/throttle3D/system/fvSchemes
index 0c8bf54bb51..3d506aa1d7f 100644
--- a/tutorials/multiphase/cavitatingFoam/les/throttle3D/system/fvSchemes
+++ b/tutorials/multiphase/cavitatingFoam/les/throttle3D/system/fvSchemes
@@ -31,6 +31,7 @@ divSchemes
     div(phiv,rho)   Gauss limitedLinear 0.2;
     div(phi,U)      Gauss filteredLinear2V 0.2 0;
     div(phiv,k)     Gauss filteredLinear2 0.2 0;
+    div((muEff*dev(T(grad(U))))) Gauss linear;
 }
 
 gradSchemes
diff --git a/tutorials/multiphase/cavitatingFoam/ras/throttle/system/fvSchemes b/tutorials/multiphase/cavitatingFoam/ras/throttle/system/fvSchemes
index 633e4a749e3..32379fb2c2b 100644
--- a/tutorials/multiphase/cavitatingFoam/ras/throttle/system/fvSchemes
+++ b/tutorials/multiphase/cavitatingFoam/ras/throttle/system/fvSchemes
@@ -32,6 +32,7 @@ divSchemes
     div(phi,U)      Gauss limitedLinearV 1;
     div(phiv,omega) Gauss limitedLinear 1;
     div(phiv,k)     Gauss limitedLinear 1;
+    div((muEff*dev(T(grad(U))))) Gauss linear;
 }
 
 gradSchemes
diff --git a/tutorials/multiphase/compressibleInterFoam/les/depthCharge2D/system/fvSchemes b/tutorials/multiphase/compressibleInterFoam/les/depthCharge2D/system/fvSchemes
index c35af6cc4f1..903d94d30c6 100644
--- a/tutorials/multiphase/compressibleInterFoam/les/depthCharge2D/system/fvSchemes
+++ b/tutorials/multiphase/compressibleInterFoam/les/depthCharge2D/system/fvSchemes
@@ -34,7 +34,7 @@ divSchemes
     div(phid2,p_rgh) Gauss upwind;
     div(rho*phi,T)  Gauss upwind;
     div(phi,k)      Gauss vanLeer;
-    div((nuEff*dev(T(grad(U))))) Gauss linear;
+    div((muEff*dev(T(grad(U))))) Gauss linear;
 }
 
 laplacianSchemes
diff --git a/tutorials/multiphase/compressibleInterFoam/les/depthCharge3D/system/fvSchemes b/tutorials/multiphase/compressibleInterFoam/les/depthCharge3D/system/fvSchemes
index c35af6cc4f1..903d94d30c6 100644
--- a/tutorials/multiphase/compressibleInterFoam/les/depthCharge3D/system/fvSchemes
+++ b/tutorials/multiphase/compressibleInterFoam/les/depthCharge3D/system/fvSchemes
@@ -34,7 +34,7 @@ divSchemes
     div(phid2,p_rgh) Gauss upwind;
     div(rho*phi,T)  Gauss upwind;
     div(phi,k)      Gauss vanLeer;
-    div((nuEff*dev(T(grad(U))))) Gauss linear;
+    div((muEff*dev(T(grad(U))))) Gauss linear;
 }
 
 laplacianSchemes
diff --git a/tutorials/multiphase/interDyMFoam/ras/damBreakWithObstacle/system/fvSchemes b/tutorials/multiphase/interDyMFoam/ras/damBreakWithObstacle/system/fvSchemes
index fa6ef27c434..ba444728e19 100644
--- a/tutorials/multiphase/interDyMFoam/ras/damBreakWithObstacle/system/fvSchemes
+++ b/tutorials/multiphase/interDyMFoam/ras/damBreakWithObstacle/system/fvSchemes
@@ -30,6 +30,7 @@ divSchemes
     div(rho*phi,U)  Gauss upwind;
     div(phi,alpha)  Gauss vanLeer;
     div(phirb,alpha) Gauss interfaceCompression;
+    div((muEff*dev(T(grad(U))))) Gauss linear;
 }
 
 laplacianSchemes
diff --git a/tutorials/multiphase/interDyMFoam/ras/floatingObject/system/fvSchemes b/tutorials/multiphase/interDyMFoam/ras/floatingObject/system/fvSchemes
index 5ce6b8ee171..852f5e0c744 100644
--- a/tutorials/multiphase/interDyMFoam/ras/floatingObject/system/fvSchemes
+++ b/tutorials/multiphase/interDyMFoam/ras/floatingObject/system/fvSchemes
@@ -32,6 +32,7 @@ divSchemes
     div(phirb,alpha) Gauss interfaceCompression;
     div(phi,k)      Gauss upwind;
     div(phi,epsilon) Gauss upwind;
+    div((muEff*dev(T(grad(U))))) Gauss linear;
 }
 
 laplacianSchemes
diff --git a/tutorials/multiphase/interDyMFoam/ras/sloshingTank2D/system/fvSchemes b/tutorials/multiphase/interDyMFoam/ras/sloshingTank2D/system/fvSchemes
index a0a60315bd4..b2c0af769de 100644
--- a/tutorials/multiphase/interDyMFoam/ras/sloshingTank2D/system/fvSchemes
+++ b/tutorials/multiphase/interDyMFoam/ras/sloshingTank2D/system/fvSchemes
@@ -30,6 +30,7 @@ divSchemes
     div(rho*phi,U)  Gauss vanLeerV;
     div(phi,alpha)  Gauss vanLeer;
     div(phirb,alpha) Gauss vanLeer;
+    div((muEff*dev(T(grad(U))))) Gauss linear;
 }
 
 laplacianSchemes
diff --git a/tutorials/multiphase/interDyMFoam/ras/sloshingTank2D3DoF/system/fvSchemes b/tutorials/multiphase/interDyMFoam/ras/sloshingTank2D3DoF/system/fvSchemes
index a0a60315bd4..b2c0af769de 100644
--- a/tutorials/multiphase/interDyMFoam/ras/sloshingTank2D3DoF/system/fvSchemes
+++ b/tutorials/multiphase/interDyMFoam/ras/sloshingTank2D3DoF/system/fvSchemes
@@ -30,6 +30,7 @@ divSchemes
     div(rho*phi,U)  Gauss vanLeerV;
     div(phi,alpha)  Gauss vanLeer;
     div(phirb,alpha) Gauss vanLeer;
+    div((muEff*dev(T(grad(U))))) Gauss linear;
 }
 
 laplacianSchemes
diff --git a/tutorials/multiphase/interDyMFoam/ras/sloshingTank3D/system/fvSchemes b/tutorials/multiphase/interDyMFoam/ras/sloshingTank3D/system/fvSchemes
index a0a60315bd4..b2c0af769de 100644
--- a/tutorials/multiphase/interDyMFoam/ras/sloshingTank3D/system/fvSchemes
+++ b/tutorials/multiphase/interDyMFoam/ras/sloshingTank3D/system/fvSchemes
@@ -30,6 +30,7 @@ divSchemes
     div(rho*phi,U)  Gauss vanLeerV;
     div(phi,alpha)  Gauss vanLeer;
     div(phirb,alpha) Gauss vanLeer;
+    div((muEff*dev(T(grad(U))))) Gauss linear;
 }
 
 laplacianSchemes
diff --git a/tutorials/multiphase/interDyMFoam/ras/sloshingTank3D3DoF/system/fvSchemes b/tutorials/multiphase/interDyMFoam/ras/sloshingTank3D3DoF/system/fvSchemes
index a0a60315bd4..b2c0af769de 100644
--- a/tutorials/multiphase/interDyMFoam/ras/sloshingTank3D3DoF/system/fvSchemes
+++ b/tutorials/multiphase/interDyMFoam/ras/sloshingTank3D3DoF/system/fvSchemes
@@ -30,6 +30,7 @@ divSchemes
     div(rho*phi,U)  Gauss vanLeerV;
     div(phi,alpha)  Gauss vanLeer;
     div(phirb,alpha) Gauss vanLeer;
+    div((muEff*dev(T(grad(U))))) Gauss linear;
 }
 
 laplacianSchemes
diff --git a/tutorials/multiphase/interDyMFoam/ras/sloshingTank3D6DoF/system/fvSchemes b/tutorials/multiphase/interDyMFoam/ras/sloshingTank3D6DoF/system/fvSchemes
index a0a60315bd4..b2c0af769de 100644
--- a/tutorials/multiphase/interDyMFoam/ras/sloshingTank3D6DoF/system/fvSchemes
+++ b/tutorials/multiphase/interDyMFoam/ras/sloshingTank3D6DoF/system/fvSchemes
@@ -30,6 +30,7 @@ divSchemes
     div(rho*phi,U)  Gauss vanLeerV;
     div(phi,alpha)  Gauss vanLeer;
     div(phirb,alpha) Gauss vanLeer;
+    div((muEff*dev(T(grad(U))))) Gauss linear;
 }
 
 laplacianSchemes
diff --git a/tutorials/multiphase/interDyMFoam/ras/testTubeMixer/system/fvSchemes b/tutorials/multiphase/interDyMFoam/ras/testTubeMixer/system/fvSchemes
index a0a60315bd4..b2c0af769de 100644
--- a/tutorials/multiphase/interDyMFoam/ras/testTubeMixer/system/fvSchemes
+++ b/tutorials/multiphase/interDyMFoam/ras/testTubeMixer/system/fvSchemes
@@ -30,6 +30,7 @@ divSchemes
     div(rho*phi,U)  Gauss vanLeerV;
     div(phi,alpha)  Gauss vanLeer;
     div(phirb,alpha) Gauss vanLeer;
+    div((muEff*dev(T(grad(U))))) Gauss linear;
 }
 
 laplacianSchemes
diff --git a/tutorials/multiphase/interFoam/laminar/capillaryRise/system/fvSchemes b/tutorials/multiphase/interFoam/laminar/capillaryRise/system/fvSchemes
index c2c563142ca..2a580ed9eb9 100644
--- a/tutorials/multiphase/interFoam/laminar/capillaryRise/system/fvSchemes
+++ b/tutorials/multiphase/interFoam/laminar/capillaryRise/system/fvSchemes
@@ -30,6 +30,7 @@ divSchemes
     div(rho*phi,U)  Gauss upwind;
     div(phi,alpha)  Gauss vanLeer;
     div(phirb,alpha) Gauss interfaceCompression;
+    div((muEff*dev(T(grad(U))))) Gauss linear;
 }
 
 laplacianSchemes
diff --git a/tutorials/multiphase/interFoam/laminar/damBreak/system/fvSchemes b/tutorials/multiphase/interFoam/laminar/damBreak/system/fvSchemes
index 7bb78c013cf..986a20eebeb 100644
--- a/tutorials/multiphase/interFoam/laminar/damBreak/system/fvSchemes
+++ b/tutorials/multiphase/interFoam/laminar/damBreak/system/fvSchemes
@@ -30,6 +30,7 @@ divSchemes
     div(rho*phi,U)  Gauss limitedLinearV 1;
     div(phi,alpha)  Gauss vanLeer;
     div(phirb,alpha) Gauss interfaceCompression;
+    div((muEff*dev(T(grad(U))))) Gauss linear;
 }
 
 laplacianSchemes
diff --git a/tutorials/multiphase/interFoam/les/nozzleFlow2D/system/fvSchemes b/tutorials/multiphase/interFoam/les/nozzleFlow2D/system/fvSchemes
index be6806ad732..26e06a656bd 100644
--- a/tutorials/multiphase/interFoam/les/nozzleFlow2D/system/fvSchemes
+++ b/tutorials/multiphase/interFoam/les/nozzleFlow2D/system/fvSchemes
@@ -34,7 +34,7 @@ divSchemes
     div(phi,B)      Gauss limitedLinear 1;
     div(B)          Gauss linear;
     div(phi,nuTilda) Gauss limitedLinear 1;
-    div((nuEff*dev(T(grad(U))))) Gauss linear;
+    div((muEff*dev(T(grad(U))))) Gauss linear;
 }
 
 laplacianSchemes
diff --git a/tutorials/multiphase/interFoam/ras/damBreak/system/fvSchemes b/tutorials/multiphase/interFoam/ras/damBreak/system/fvSchemes
index 28b3d7d20ae..8ed8e03cc08 100644
--- a/tutorials/multiphase/interFoam/ras/damBreak/system/fvSchemes
+++ b/tutorials/multiphase/interFoam/ras/damBreak/system/fvSchemes
@@ -35,7 +35,7 @@ divSchemes
     div(phi,R)      Gauss upwind;
     div(R)          Gauss linear;
     div(phi,nuTilda) Gauss upwind;
-    div((nuEff*dev(T(grad(U))))) Gauss linear;
+    div((muEff*dev(T(grad(U))))) Gauss linear;
 }
 
 laplacianSchemes
diff --git a/tutorials/multiphase/interFoam/ras/damBreakPorousBaffle/system/fvSchemes b/tutorials/multiphase/interFoam/ras/damBreakPorousBaffle/system/fvSchemes
index 28b3d7d20ae..8ed8e03cc08 100644
--- a/tutorials/multiphase/interFoam/ras/damBreakPorousBaffle/system/fvSchemes
+++ b/tutorials/multiphase/interFoam/ras/damBreakPorousBaffle/system/fvSchemes
@@ -35,7 +35,7 @@ divSchemes
     div(phi,R)      Gauss upwind;
     div(R)          Gauss linear;
     div(phi,nuTilda) Gauss upwind;
-    div((nuEff*dev(T(grad(U))))) Gauss linear;
+    div((muEff*dev(T(grad(U))))) Gauss linear;
 }
 
 laplacianSchemes
diff --git a/tutorials/multiphase/interFoam/ras/waterChannel/LTSInterFoam/system/fvSchemes b/tutorials/multiphase/interFoam/ras/waterChannel/LTSInterFoam/system/fvSchemes
index 22366f437b7..7018c5eaa73 100644
--- a/tutorials/multiphase/interFoam/ras/waterChannel/LTSInterFoam/system/fvSchemes
+++ b/tutorials/multiphase/interFoam/ras/waterChannel/LTSInterFoam/system/fvSchemes
@@ -32,6 +32,7 @@ divSchemes
     div(phirb,alpha) Gauss interfaceCompression;
     div(phi,k)      Gauss upwind;
     div(phi,omega)  Gauss upwind;
+    div((muEff*dev(T(grad(U))))) Gauss linear;
 }
 
 laplacianSchemes
diff --git a/tutorials/multiphase/interFoam/ras/waterChannel/system/fvSchemes b/tutorials/multiphase/interFoam/ras/waterChannel/system/fvSchemes
index 43847dd8414..bb536c017ce 100644
--- a/tutorials/multiphase/interFoam/ras/waterChannel/system/fvSchemes
+++ b/tutorials/multiphase/interFoam/ras/waterChannel/system/fvSchemes
@@ -33,6 +33,7 @@ divSchemes
 
     div(phi,k)      Gauss upwind;
     div(phi,omega)  $div(phi,k);
+    div((muEff*dev(T(grad(U))))) Gauss linear;
 }
 
 laplacianSchemes
diff --git a/tutorials/multiphase/interFoam/ras/weirOverflow/system/fvSchemes b/tutorials/multiphase/interFoam/ras/weirOverflow/system/fvSchemes
index 28b3d7d20ae..8ed8e03cc08 100644
--- a/tutorials/multiphase/interFoam/ras/weirOverflow/system/fvSchemes
+++ b/tutorials/multiphase/interFoam/ras/weirOverflow/system/fvSchemes
@@ -35,7 +35,7 @@ divSchemes
     div(phi,R)      Gauss upwind;
     div(R)          Gauss linear;
     div(phi,nuTilda) Gauss upwind;
-    div((nuEff*dev(T(grad(U))))) Gauss linear;
+    div((muEff*dev(T(grad(U))))) Gauss linear;
 }
 
 laplacianSchemes
diff --git a/tutorials/multiphase/interMixingFoam/laminar/damBreak/system/fvSchemes b/tutorials/multiphase/interMixingFoam/laminar/damBreak/system/fvSchemes
index ed937573ff9..c80ccb314ef 100644
--- a/tutorials/multiphase/interMixingFoam/laminar/damBreak/system/fvSchemes
+++ b/tutorials/multiphase/interMixingFoam/laminar/damBreak/system/fvSchemes
@@ -30,6 +30,7 @@ divSchemes
     div(rho*phi,U)  Gauss limitedLinearV 1;
     div(phi,alpha)  Gauss vanLeer;
     div(phirb,alpha) Gauss interfaceCompression;
+    div((muEff*dev(T(grad(U))))) Gauss linear;
 }
 
 laplacianSchemes
diff --git a/tutorials/multiphase/interPhaseChangeFoam/cavitatingBullet/system/fvSchemes b/tutorials/multiphase/interPhaseChangeFoam/cavitatingBullet/system/fvSchemes
index c1b7f712065..44313d502d0 100644
--- a/tutorials/multiphase/interPhaseChangeFoam/cavitatingBullet/system/fvSchemes
+++ b/tutorials/multiphase/interPhaseChangeFoam/cavitatingBullet/system/fvSchemes
@@ -30,9 +30,9 @@ divSchemes
     div(rhoPhi,U)        Gauss linearUpwind grad(U);
     div(phi,omega)       Gauss linearUpwind grad(omega);
     div(phi,k)           Gauss linearUpwind grad(k);
-
     div(phi,alpha)       Gauss vanLeer;
     div(phirb,alpha)     Gauss interfaceCompression;
+    div((muEff*dev(T(grad(U))))) Gauss linear;
 }
 
 gradSchemes
diff --git a/tutorials/multiphase/multiphaseInterFoam/laminar/damBreak4phase/system/fvSchemes b/tutorials/multiphase/multiphaseInterFoam/laminar/damBreak4phase/system/fvSchemes
index 6e3d6f1b848..e1436932da8 100644
--- a/tutorials/multiphase/multiphaseInterFoam/laminar/damBreak4phase/system/fvSchemes
+++ b/tutorials/multiphase/multiphaseInterFoam/laminar/damBreak4phase/system/fvSchemes
@@ -32,6 +32,7 @@ divSchemes
     div(rho*phi,U)  Gauss upwind;
     div(phi,alpha)  Gauss vanLeer;
     div(phirb,alpha) Gauss interfaceCompression;
+    div((muEff*dev(T(grad(U))))) Gauss linear;
 }
 
 laplacianSchemes
diff --git a/tutorials/multiphase/multiphaseInterFoam/laminar/damBreak4phaseFine/system/fvSchemes b/tutorials/multiphase/multiphaseInterFoam/laminar/damBreak4phaseFine/system/fvSchemes
index 6e3d6f1b848..e1436932da8 100644
--- a/tutorials/multiphase/multiphaseInterFoam/laminar/damBreak4phaseFine/system/fvSchemes
+++ b/tutorials/multiphase/multiphaseInterFoam/laminar/damBreak4phaseFine/system/fvSchemes
@@ -32,6 +32,7 @@ divSchemes
     div(rho*phi,U)  Gauss upwind;
     div(phi,alpha)  Gauss vanLeer;
     div(phirb,alpha) Gauss interfaceCompression;
+    div((muEff*dev(T(grad(U))))) Gauss linear;
 }
 
 laplacianSchemes
diff --git a/tutorials/multiphase/twoLiquidMixingFoam/lockExchange/system/fvSchemes b/tutorials/multiphase/twoLiquidMixingFoam/lockExchange/system/fvSchemes
index 4bb0230588c..35930b92037 100644
--- a/tutorials/multiphase/twoLiquidMixingFoam/lockExchange/system/fvSchemes
+++ b/tutorials/multiphase/twoLiquidMixingFoam/lockExchange/system/fvSchemes
@@ -33,7 +33,7 @@ divSchemes
     div(rho*phi,U)  Gauss linear;
     div(phi,alpha1) Gauss vanLeer;
     div(phi,k)      Gauss limitedLinear 1;
-    div(((rho*nuEff)*dev(grad(U).T()))) Gauss linear;
+    div((muEff*dev(T(grad(U))))) Gauss linear;
 }
 
 laplacianSchemes
-- 
GitLab