diff --git a/src/TurbulenceModels/turbulenceModels/ReynoldsStress/ReynoldsStress.C b/src/TurbulenceModels/turbulenceModels/ReynoldsStress/ReynoldsStress.C
index 0164074765..b9f91bd180 100644
--- a/src/TurbulenceModels/turbulenceModels/ReynoldsStress/ReynoldsStress.C
+++ b/src/TurbulenceModels/turbulenceModels/ReynoldsStress/ReynoldsStress.C
@@ -227,9 +227,11 @@ Foam::ReynoldsStress<BasicTurbulenceModel>::devRhoReff() const
 
 
 template<class BasicTurbulenceModel>
+template<class RhoFieldType>
 Foam::tmp<Foam::fvVectorMatrix>
-Foam::ReynoldsStress<BasicTurbulenceModel>::divDevRhoReff
+Foam::ReynoldsStress<BasicTurbulenceModel>::DivDevRhoReff
 (
+    const RhoFieldType& rho,
     volVectorField& U
 ) const
 {
@@ -239,19 +241,19 @@ Foam::ReynoldsStress<BasicTurbulenceModel>::divDevRhoReff
         (
             fvc::laplacian
             (
-                (1.0 - couplingFactor_)*this->alpha_*this->rho_*this->nut(),
+                (1.0 - couplingFactor_)*this->alpha_*rho*this->nut(),
                 U,
                 "laplacian(nuEff,U)"
             )
           + fvc::div
             (
-                this->alpha_*this->rho_*R_
+                this->alpha_*rho*R_
               + couplingFactor_
-               *this->alpha_*this->rho_*this->nut()*fvc::grad(U),
+               *this->alpha_*rho*this->nut()*fvc::grad(U),
                 "div(devRhoReff)"
             )
-          - fvc::div(this->alpha_*this->rho_*this->nu()*dev2(T(fvc::grad(U))))
-          - fvm::laplacian(this->alpha_*this->rho_*this->nuEff(), U)
+          - fvc::div(this->alpha_*rho*this->nu()*dev2(T(fvc::grad(U))))
+          - fvm::laplacian(this->alpha_*rho*this->nuEff(), U)
         );
     }
     else
@@ -260,21 +262,26 @@ Foam::ReynoldsStress<BasicTurbulenceModel>::divDevRhoReff
         (
             fvc::laplacian
             (
-                this->alpha_*this->rho_*this->nut(),
+                this->alpha_*rho*this->nut(),
                 U,
                 "laplacian(nuEff,U)"
             )
-          + fvc::div(this->alpha_*this->rho_*R_)
-          - fvc::div(this->alpha_*this->rho_*this->nu()*dev2(T(fvc::grad(U))))
-          - fvm::laplacian(this->alpha_*this->rho_*this->nuEff(), U)
+          + fvc::div(this->alpha_*rho*R_)
+          - fvc::div(this->alpha_*rho*this->nu()*dev2(T(fvc::grad(U))))
+          - fvm::laplacian(this->alpha_*rho*this->nuEff(), U)
         );
     }
+}
 
-    return
-    (
-      - fvc::div((this->alpha_*this->rho_*this->nuEff())*dev2(T(fvc::grad(U))))
-      - fvm::laplacian(this->alpha_*this->rho_*this->nuEff(), U)
-    );
+
+template<class BasicTurbulenceModel>
+Foam::tmp<Foam::fvVectorMatrix>
+Foam::ReynoldsStress<BasicTurbulenceModel>::divDevRhoReff
+(
+    volVectorField& U
+) const
+{
+    return DivDevRhoReff(this->rho_, U);
 }
 
 
@@ -286,11 +293,7 @@ Foam::ReynoldsStress<BasicTurbulenceModel>::divDevRhoReff
     volVectorField& U
 ) const
 {
-    return
-    (
-      - fvc::div((this->alpha_*rho*this->nuEff())*dev2(T(fvc::grad(U))))
-      - fvm::laplacian(this->alpha_*rho*this->nuEff(), U)
-    );
+    return DivDevRhoReff(rho, U);
 }
 
 
diff --git a/src/TurbulenceModels/turbulenceModels/ReynoldsStress/ReynoldsStress.H b/src/TurbulenceModels/turbulenceModels/ReynoldsStress/ReynoldsStress.H
index e29c447459..084830b0d6 100644
--- a/src/TurbulenceModels/turbulenceModels/ReynoldsStress/ReynoldsStress.H
+++ b/src/TurbulenceModels/turbulenceModels/ReynoldsStress/ReynoldsStress.H
@@ -75,6 +75,14 @@ protected:
         //- Update the eddy-viscosity
         virtual void correctNut() = 0;
 
+        //- Return the source term for the momentum equation
+        template<class RhoFieldType>
+        tmp<fvVectorMatrix> DivDevRhoReff
+        (
+            const RhoFieldType& rho,
+            volVectorField& U
+        ) const;
+
 
 public: