diff --git a/applications/solvers/multiphase/cavitatingFoam/UEqn.H b/applications/solvers/multiphase/cavitatingFoam/UEqn.H
index 9a5761b59f3af176c34fde4f0b366573cf510cfe..1ca0b9f0a5ffd09df362b53148b291968f9ea4f3 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 257f6d48b53176b612a8732887d32a4f1f767905..7cc250a66a2ab15a22680ab352f321e731b17c17 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 257f6d48b53176b612a8732887d32a4f1f767905..7cc250a66a2ab15a22680ab352f321e731b17c17 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 b8fe82ff1423ed793c82bdcb0a5d43d226be0bc4..062c5523c9d4c94a469e49e202b65e25ad1b285d 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 43ee1e3eb9452061200f3908241358d3f3215a24..fbcba7db72d751ac5b48d8055789b54eb8877b73 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 d6232da309d69b324c0eb327ac7f11287df7fb57..ff6323bca2776b6e06cef19d4ad56209f022c015 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 3dbb3747743536d11d95ce2021c81cab1641ea3f..7afd323dd0d0b373d36b96a1ecc20df73920401b 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 1dd26906b74487d46da744e23d9f1c2a68dc67cd..ff9ba3a4e3e749c46189eb80334f68df95f95d53 100644
--- a/src/postProcessing/functionObjects/forces/forces/forces.C
+++ b/src/postProcessing/functionObjects/forces/forces/forces.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -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 0c905a91cedb52c7953accfd971beeb04b320dcb..69795c56d7af1d6edaebddf0361e1c8f66063ae2 100644
--- a/src/turbulenceModels/incompressible/LES/GenEddyVisc/GenEddyVisc.C
+++ b/src/turbulenceModels/incompressible/LES/GenEddyVisc/GenEddyVisc.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -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 0746377eb79ece413a09ad6ce6229f45109f6779..62c216ee3ec06fd13b84ef69c3baff5d2991d977 100644
--- a/src/turbulenceModels/incompressible/LES/GenEddyVisc/GenEddyVisc.H
+++ b/src/turbulenceModels/incompressible/LES/GenEddyVisc/GenEddyVisc.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -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 9d385c21eb1687c6e72d5e43b965716c746deb95..1618fb7f3be0a3f9b0e96e2dff12c4bd08855207 100644
--- a/src/turbulenceModels/incompressible/LES/GenSGSStress/GenSGSStress.C
+++ b/src/turbulenceModels/incompressible/LES/GenSGSStress/GenSGSStress.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -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 36d84f9c771a1d6f3aa23c4a107201acbf8b0e5f..0456d69a7178193ffa6b0f05cea10f521796c247 100644
--- a/src/turbulenceModels/incompressible/LES/GenSGSStress/GenSGSStress.H
+++ b/src/turbulenceModels/incompressible/LES/GenSGSStress/GenSGSStress.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -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 deca26e7936f4a692adc8d0066cf08ba27d1a4ba..3688cb8dec91781b9b33c7f28a960630f16116a3 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 97d32bce9dd2ac39bc2fb3184487085be42132e9..84fc1fc7be5cfdf412cb9838089552668770af87 100644
--- a/src/turbulenceModels/incompressible/LES/Smagorinsky2/Smagorinsky2.C
+++ b/src/turbulenceModels/incompressible/LES/Smagorinsky2/Smagorinsky2.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -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 89c77912cde1bccf4610dff930ff5f6224d7babd..c9133e4e8732cf7e265552fedaf276f65de7a7ae 100644
--- a/src/turbulenceModels/incompressible/LES/Smagorinsky2/Smagorinsky2.H
+++ b/src/turbulenceModels/incompressible/LES/Smagorinsky2/Smagorinsky2.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -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 6a34f90892e0ae075f022959b87b516c3b9e8037..da8ab55e07b45492e216296381b0f4cd89a7d78d 100644
--- a/src/turbulenceModels/incompressible/LES/SpalartAllmaras/SpalartAllmaras.C
+++ b/src/turbulenceModels/incompressible/LES/SpalartAllmaras/SpalartAllmaras.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -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 7ed507712a2eafc893ca10cdbb2abd1faaac8b88..68e5ef4c73bd60cc5d5f54afeff2b467d30c07e4 100644
--- a/src/turbulenceModels/incompressible/LES/SpalartAllmaras/SpalartAllmaras.H
+++ b/src/turbulenceModels/incompressible/LES/SpalartAllmaras/SpalartAllmaras.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -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 2ab771ac73730b236d71ae1158aeb8ed17665f26..12a6f1b3d111b051302587ca94d4709c24f0077a 100644
--- a/src/turbulenceModels/incompressible/LES/kOmegaSSTSAS/kOmegaSSTSAS.C
+++ b/src/turbulenceModels/incompressible/LES/kOmegaSSTSAS/kOmegaSSTSAS.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -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 2ca98380afad07e6d762abd6ba539e23e5aa5a9d..e61156f9c9081626d85018a744f76549497587ad 100644
--- a/src/turbulenceModels/incompressible/LES/kOmegaSSTSAS/kOmegaSSTSAS.H
+++ b/src/turbulenceModels/incompressible/LES/kOmegaSSTSAS/kOmegaSSTSAS.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -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 c88ff32beba4648761bf68b8f8d669b816e17257..5a78c213ee3f78201d77366f2e2bf54b38590c5c 100644
--- a/src/turbulenceModels/incompressible/LES/laminar/laminar.C
+++ b/src/turbulenceModels/incompressible/LES/laminar/laminar.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -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 14de7a142ac90302938c5196c5d04ab16ebf7a5e..777a9228a73e61fd0fe09724297142515347de1e 100644
--- a/src/turbulenceModels/incompressible/LES/laminar/laminar.H
+++ b/src/turbulenceModels/incompressible/LES/laminar/laminar.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -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 4f2804577fa4020e007e222bdc946badee3e3715..b17f4ea572a2a0ab80ecc67e567f851ddfdbcc62 100644
--- a/src/turbulenceModels/incompressible/LES/mixedSmagorinsky/mixedSmagorinsky.C
+++ b/src/turbulenceModels/incompressible/LES/mixedSmagorinsky/mixedSmagorinsky.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -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 c8557db2391af5e0f45fece323d0c75ae55b8ff3..72dca5d35484ab6929c7635e3b09f70b08df3f58 100644
--- a/src/turbulenceModels/incompressible/LES/mixedSmagorinsky/mixedSmagorinsky.H
+++ b/src/turbulenceModels/incompressible/LES/mixedSmagorinsky/mixedSmagorinsky.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -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 b253573be13a387b48239fae3bd0359e2060b968..61322de356f08e3401d7454eeb442ec26b7a4d70 100644
--- a/src/turbulenceModels/incompressible/LES/scaleSimilarity/scaleSimilarity.C
+++ b/src/turbulenceModels/incompressible/LES/scaleSimilarity/scaleSimilarity.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -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 0dff31d87e302b3970b00efc51fa39bc78f033a0..c757a8165e2f266c1c0f9e01407ed62be8c59652 100644
--- a/src/turbulenceModels/incompressible/LES/scaleSimilarity/scaleSimilarity.H
+++ b/src/turbulenceModels/incompressible/LES/scaleSimilarity/scaleSimilarity.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -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 30afa63ce9964921dff160aeffa0f2fc38f2b00a..c2c7e5efbe456224c622afbfb168183dfb430dae 100644
--- a/src/turbulenceModels/incompressible/RAS/LRR/LRR.C
+++ b/src/turbulenceModels/incompressible/RAS/LRR/LRR.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -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 2f1a119b9ca061942c806b74b566fed8c48aa251..e7873560081a0c7ad1755f69b5a2dafd383fea49 100644
--- a/src/turbulenceModels/incompressible/RAS/LRR/LRR.H
+++ b/src/turbulenceModels/incompressible/RAS/LRR/LRR.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -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 1671c625a083b67609baccbb85d59592f925706a..10de46f58703089b06356d9265b17ccde0efdc8b 100644
--- a/src/turbulenceModels/incompressible/RAS/LamBremhorstKE/LamBremhorstKE.C
+++ b/src/turbulenceModels/incompressible/RAS/LamBremhorstKE/LamBremhorstKE.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -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 449055f4c5dde75b5f4e63ab52eda72b9cc89ee7..3ae7cfda29753bc9f6d960bb213ddee40360f35c 100644
--- a/src/turbulenceModels/incompressible/RAS/LamBremhorstKE/LamBremhorstKE.H
+++ b/src/turbulenceModels/incompressible/RAS/LamBremhorstKE/LamBremhorstKE.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -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 aafa62af19292d35f01cae9828208b453614932a..49326b1ff50c27f6453223ec8f489c590b13cce3 100644
--- a/src/turbulenceModels/incompressible/RAS/LaunderGibsonRSTM/LaunderGibsonRSTM.C
+++ b/src/turbulenceModels/incompressible/RAS/LaunderGibsonRSTM/LaunderGibsonRSTM.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -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 0b6fc91b4bb16c7736af1e1b637617b67286677b..528a7dd5d6d0c0556b32958600da3050e9a857b8 100644
--- a/src/turbulenceModels/incompressible/RAS/LaunderGibsonRSTM/LaunderGibsonRSTM.H
+++ b/src/turbulenceModels/incompressible/RAS/LaunderGibsonRSTM/LaunderGibsonRSTM.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -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 f0fa6e60c8d3abc3fd363025c6f45d39e9e85afb..c514a378e20c4b268b229e1960fb4084d3e777b4 100644
--- a/src/turbulenceModels/incompressible/RAS/LaunderSharmaKE/LaunderSharmaKE.C
+++ b/src/turbulenceModels/incompressible/RAS/LaunderSharmaKE/LaunderSharmaKE.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -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 3f2338e18e837675fecca62cd9bc3c47ed1755de..5e4ae36f5ae2f7103536f5056563ea5ecc9ece02 100644
--- a/src/turbulenceModels/incompressible/RAS/LaunderSharmaKE/LaunderSharmaKE.H
+++ b/src/turbulenceModels/incompressible/RAS/LaunderSharmaKE/LaunderSharmaKE.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -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 03e1e5aa1acb2673e463975a1b8483cf47cab590..20d8dee95140a3c513c4c3f1729266176f63c844 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 ebe7c373d8e9f0d9b96bfcc4e8077f16aa3c07b0..56696f42539287ccac73766159ee3a2a02e81c60 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 e0bd0f0df9f76afb0b15821460b37463f405903f..da81b66edddd428183e706efb5f1357fa58de2f4 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 0efd63a454c5f5cbe04e0df94c86ae877cc60271..8794af1fc4bec51c8f73f42a5ba7d331ed2146e8 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 91e288f104258375661e2d1fda755d9e13c4b36b..15296cdf5bb03f0b80b576d58ffab1629baceca5 100644
--- a/src/turbulenceModels/incompressible/RAS/LienLeschzinerLowRe/LienLeschzinerLowRe.C
+++ b/src/turbulenceModels/incompressible/RAS/LienLeschzinerLowRe/LienLeschzinerLowRe.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -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 3fa1c790fe0cf159a72768d10ee94fad02c60669..5aebe9ea50a30f3aaa84b8a54e9c3f5ae59eb728 100644
--- a/src/turbulenceModels/incompressible/RAS/LienLeschzinerLowRe/LienLeschzinerLowRe.H
+++ b/src/turbulenceModels/incompressible/RAS/LienLeschzinerLowRe/LienLeschzinerLowRe.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -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 33cc03edac8927bb0aa577b1bfef6b92870ad940..4dbb7b7b6832155c1d6f4d85829e2dd6235a01c5 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 c7ad8a371313005e06ff2f9499ebbc803d9dd994..f63f0c6565c3300fd96378f7294a8c9af81a0712 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 794d9f9b07f10a4888eb2b27ec7b671eaa8d0a4f..460a6609acdd01f2497333776058deaf808cf90a 100644
--- a/src/turbulenceModels/incompressible/RAS/RNGkEpsilon/RNGkEpsilon.C
+++ b/src/turbulenceModels/incompressible/RAS/RNGkEpsilon/RNGkEpsilon.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -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 62b60b98046b1265bf0e7aae1426a349c3c6ee36..2d06ee9ec1a804e2cb68793698e974aa12158348 100644
--- a/src/turbulenceModels/incompressible/RAS/RNGkEpsilon/RNGkEpsilon.H
+++ b/src/turbulenceModels/incompressible/RAS/RNGkEpsilon/RNGkEpsilon.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -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 1804ba5b01007b741f4ad25a14c612032c0940d4..4105cb1876e1664fd74006e7193ccd895f01b951 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 fbd1248181036cdf5116a3ab4d6d1cdd2a103039..4190fc7d8090fe1429a02bf9917258350abae941 100644
--- a/src/turbulenceModels/incompressible/RAS/SpalartAllmaras/SpalartAllmaras.H
+++ b/src/turbulenceModels/incompressible/RAS/SpalartAllmaras/SpalartAllmaras.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -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 99b4188cc016f2343d582fcfa798b602151f0c00..029b9f9db790f4a44b555f4e01182741e16c3fe4 100644
--- a/src/turbulenceModels/incompressible/RAS/kEpsilon/kEpsilon.C
+++ b/src/turbulenceModels/incompressible/RAS/kEpsilon/kEpsilon.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -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 a3baf20513454efadb9130ec30f12cf5aad35353..89c94fe42faff1c5a885688d05a33d41c2ad325d 100644
--- a/src/turbulenceModels/incompressible/RAS/kEpsilon/kEpsilon.H
+++ b/src/turbulenceModels/incompressible/RAS/kEpsilon/kEpsilon.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -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 7106c0ea2cd234b78401e0fccdbed220365043ff..51e1fbd491963020544af6770dc5c282d574cee8 100644
--- a/src/turbulenceModels/incompressible/RAS/kOmega/kOmega.C
+++ b/src/turbulenceModels/incompressible/RAS/kOmega/kOmega.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -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 10729f334e714713cfb9f710530415bb6a2e2554..1e63fd2eca384a027fbbdf879a7d1d8fbc2aea6d 100644
--- a/src/turbulenceModels/incompressible/RAS/kOmega/kOmega.H
+++ b/src/turbulenceModels/incompressible/RAS/kOmega/kOmega.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -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 c025076baa9936cfc78b97ff206c7a8d1f468b5b..e6c49db1c86724808fed61940c3d9ac16246383b 100644
--- a/src/turbulenceModels/incompressible/RAS/kOmegaSST/kOmegaSST.C
+++ b/src/turbulenceModels/incompressible/RAS/kOmegaSST/kOmegaSST.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -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 bd4b3742b5e4ee48b042f12bd1a8da8f129f3f7c..60f2775b1b0d618665c41d2b880e689ef0f54183 100644
--- a/src/turbulenceModels/incompressible/RAS/kOmegaSST/kOmegaSST.H
+++ b/src/turbulenceModels/incompressible/RAS/kOmegaSST/kOmegaSST.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -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 abf9b29cde745690e802753df899e10c57595a14..45f4f69d1302091cc5bec564af3c8d1f3680b97a 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 83db5f82d5cb36a8682dc2c44c3985ba86960226..dbd5ce0832a09f23c5ab3411933f9eb383c3f721 100644
--- a/src/turbulenceModels/incompressible/RAS/kkLOmega/kkLOmega.H
+++ b/src/turbulenceModels/incompressible/RAS/kkLOmega/kkLOmega.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -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 68adca77b3f79374264bf36cd4316b78543aa62d..f0597b58da5f3fd4fe7c692f391b857a5072b238 100644
--- a/src/turbulenceModels/incompressible/RAS/laminar/laminar.C
+++ b/src/turbulenceModels/incompressible/RAS/laminar/laminar.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -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 90afadff1329a512ea35ca91886e96635cf0f925..ccebe774ed171d0b24216327cc2e735bb676dea5 100644
--- a/src/turbulenceModels/incompressible/RAS/laminar/laminar.H
+++ b/src/turbulenceModels/incompressible/RAS/laminar/laminar.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -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 25940abb9ef0c81be459ea2effc2f19dffc3f8b3..1fd761b1d7e9e23e5b334684add0937e4cdaa7f3 100644
--- a/src/turbulenceModels/incompressible/RAS/qZeta/qZeta.C
+++ b/src/turbulenceModels/incompressible/RAS/qZeta/qZeta.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -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 ba8dc1b413f76bf329855600808da3e17c2da946..9347998a2e94f23b7fa50d4ac21924789d9accd9 100644
--- a/src/turbulenceModels/incompressible/RAS/qZeta/qZeta.H
+++ b/src/turbulenceModels/incompressible/RAS/qZeta/qZeta.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -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 71058caa7d9472fb2667c09183e1d4e4349e7bcb..495b8930dfc6b0df2b18c350e3725af06a65fe0d 100644
--- a/src/turbulenceModels/incompressible/RAS/realizableKE/realizableKE.C
+++ b/src/turbulenceModels/incompressible/RAS/realizableKE/realizableKE.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -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 045f33339088825e8caa8d670d5a12834f31c50e..ed2bc4841e1e79cfa7eb05dfc3d3be5ef7142d2f 100644
--- a/src/turbulenceModels/incompressible/RAS/realizableKE/realizableKE.H
+++ b/src/turbulenceModels/incompressible/RAS/realizableKE/realizableKE.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -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 c102b52cb1e10e2ae7a5525b590be169360b9a1c..a633c7a70e30df29bdc8c7feb1b33e4fa86d3ffc 100644
--- a/src/turbulenceModels/incompressible/turbulenceModel/laminar/laminar.C
+++ b/src/turbulenceModels/incompressible/turbulenceModel/laminar/laminar.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -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 72994e90fb84908cb4f6a29e93a702a8cfc0afb6..5b79ecd48538d6226854350183a0746e7ab9b5ca 100644
--- a/src/turbulenceModels/incompressible/turbulenceModel/laminar/laminar.H
+++ b/src/turbulenceModels/incompressible/turbulenceModel/laminar/laminar.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -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 405cfef0bfa47b20b1e34765820b0deb0ba7598f..c2eacdaf39201458b69422646d363383a7fde7f3 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 22366f437b7bcaa5224b1effab5efdaeb66d0ba5..7018c5eaa736966eedef0d3358cf514187949c90 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 41159c5be641e0b096ffce563f9b4a105ed2a7a7..3827ac6548c856001186bea04d81494e454c8392 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 b4ad45237d10b2f82fe10d6a8e7f4eadd2356b5d..d1ed26cae3bca94b43c279390ae992152d8c1682 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 0c8bf54bb5147c6e3efe5a0212ddbf9d31f23449..3d506aa1d7f2e631efb8635916edcd7f8dcc0ab1 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 0c8bf54bb5147c6e3efe5a0212ddbf9d31f23449..3d506aa1d7f2e631efb8635916edcd7f8dcc0ab1 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 633e4a749e31ff2499e39cf4bc8b0ef6b26c2dae..32379fb2c2bd82a8a82c7ce40705b1c31dba608f 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 c35af6cc4f1f5d6fee1dbd9232c03ad306fbbec1..903d94d30c6e6ac6ecc85e0d048a668f7d3cb220 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 c35af6cc4f1f5d6fee1dbd9232c03ad306fbbec1..903d94d30c6e6ac6ecc85e0d048a668f7d3cb220 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 fa6ef27c434ef90089fa8009d6f23362c604aaf9..ba444728e19c244ad788ce3878e4c65a2609b292 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 5ce6b8ee1718544482f1496e5c24cc8275cb6e34..852f5e0c74425f84b7f2d62fe4dda0a2ab1fec4a 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 a0a60315bd4849d92792125419fdb9570653849d..b2c0af769de60be4fa5704d4973f1d0e6ae62d45 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 a0a60315bd4849d92792125419fdb9570653849d..b2c0af769de60be4fa5704d4973f1d0e6ae62d45 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 a0a60315bd4849d92792125419fdb9570653849d..b2c0af769de60be4fa5704d4973f1d0e6ae62d45 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 a0a60315bd4849d92792125419fdb9570653849d..b2c0af769de60be4fa5704d4973f1d0e6ae62d45 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 a0a60315bd4849d92792125419fdb9570653849d..b2c0af769de60be4fa5704d4973f1d0e6ae62d45 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 a0a60315bd4849d92792125419fdb9570653849d..b2c0af769de60be4fa5704d4973f1d0e6ae62d45 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 c2c563142ca051a8a9bcb924df1baa06b8264939..2a580ed9eb90ab5088bf09e7caf69fd7ba0ac138 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 7bb78c013cfdfd415ca1868448b5a7739cf2814b..986a20eebeb97fe96e535bdd4b6ba3a632f1a4b3 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 be6806ad7320924d7a785e8e826a48d3f45e549c..26e06a656bd68191f6c494cdb2f9a44f7b785205 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 28b3d7d20aeb00f17348671a1a85ba2700c6f874..8ed8e03cc081883369e850bae8f758e5823c4b5b 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 28b3d7d20aeb00f17348671a1a85ba2700c6f874..8ed8e03cc081883369e850bae8f758e5823c4b5b 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 22366f437b7bcaa5224b1effab5efdaeb66d0ba5..7018c5eaa736966eedef0d3358cf514187949c90 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 43847dd8414adf639354202fa7159ef5620866ee..bb536c017ce4cef960716bc13d34c040d9f6d17a 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 28b3d7d20aeb00f17348671a1a85ba2700c6f874..8ed8e03cc081883369e850bae8f758e5823c4b5b 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 ed937573ff92bb26dc2313204024c0bd69473172..c80ccb314ef4160a628f59e7e4ffc6d97caa3dd8 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 c1b7f7120651910e965db385aebc9c9757227b01..44313d502d0f1656f7c94ed9109b2673d16b600c 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 6e3d6f1b848066c7a0ccfb94e97c541fcea03398..e1436932da80a557362978007f606cf262f31dd5 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 6e3d6f1b848066c7a0ccfb94e97c541fcea03398..e1436932da80a557362978007f606cf262f31dd5 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 4bb0230588cf9b264ce877886a7bcec8f595b18e..35930b92037ebce04ce8e97605ce4f145759b445 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