diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/turbulentDispersionModels/Burns/Burns.C b/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/turbulentDispersionModels/Burns/Burns.C
index 7fa4064980ef75eed6f6aed9b24eefb3971e16f4..1de7c87de598346824d8f14de52f4007e1c1d7a2 100644
--- a/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/turbulentDispersionModels/Burns/Burns.C
+++ b/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/turbulentDispersionModels/Burns/Burns.C
@@ -25,7 +25,6 @@ License
 
 #include "Burns.H"
 #include "phasePair.H"
-#include "fvc.H"
 #include "PhaseCompressibleTurbulenceModel.H"
 #include "addToRunTimeSelectionTable.H"
 
@@ -70,8 +69,8 @@ Foam::turbulentDispersionModels::Burns::~Burns()
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-Foam::tmp<Foam::volVectorField>
-Foam::turbulentDispersionModels::Burns::F() const
+Foam::tmp<Foam::volScalarField>
+Foam::turbulentDispersionModels::Burns::Fprime() const
 {
     const fvMesh& mesh(pair_.phase1().mesh());
     const dragModel&
@@ -93,7 +92,6 @@ Foam::turbulentDispersionModels::Burns::F() const
            *sqr(pair_.dispersed().d())
         )
        *pair_.continuous().rho()
-       *fvc::grad(pair_.continuous())
        *(1.0 + pair_.dispersed()/max(pair_.continuous(), residualAlpha_));
 }
 
diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/turbulentDispersionModels/Burns/Burns.H b/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/turbulentDispersionModels/Burns/Burns.H
index edd33dc5990a076228e93841a13a9cbadfe23575..0d414bd5f3deaf92698c38b5aeff3190355aa3fb 100644
--- a/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/turbulentDispersionModels/Burns/Burns.H
+++ b/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/turbulentDispersionModels/Burns/Burns.H
@@ -103,8 +103,9 @@ public:
 
     // Member Functions
 
-        //- Turbulent dispersion force
-        virtual tmp<volVectorField> F() const;
+        //- Turbulent dispersion force coefficient
+        //  multiplying the gradient of the phase-fraction
+        virtual tmp<volScalarField> Fprime() const;
 };
 
 
diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/turbulentDispersionModels/Gosman/Gosman.C b/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/turbulentDispersionModels/Gosman/Gosman.C
index 7f7cc56526c95930947414dd06f2136fdc0c345c..005ef7c1a59ae9d5c12e3b6f8a6da4b08ff3b288 100644
--- a/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/turbulentDispersionModels/Gosman/Gosman.C
+++ b/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/turbulentDispersionModels/Gosman/Gosman.C
@@ -25,7 +25,6 @@ License
 
 #include "Gosman.H"
 #include "phasePair.H"
-#include "fvc.H"
 #include "PhaseCompressibleTurbulenceModel.H"
 #include "addToRunTimeSelectionTable.H"
 
@@ -69,8 +68,8 @@ Foam::turbulentDispersionModels::Gosman::~Gosman()
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-Foam::tmp<Foam::volVectorField>
-Foam::turbulentDispersionModels::Gosman::F() const
+Foam::tmp<Foam::volScalarField>
+Foam::turbulentDispersionModels::Gosman::Fprime() const
 {
     const fvMesh& mesh(pair_.phase1().mesh());
     const dragModel&
@@ -92,8 +91,7 @@ Foam::turbulentDispersionModels::Gosman::F() const
             sigma_
            *sqr(pair_.dispersed().d())
         )
-       *pair_.continuous().rho()
-       *fvc::grad(pair_.dispersed());
+       *pair_.continuous().rho();
 }
 
 
diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/turbulentDispersionModels/Gosman/Gosman.H b/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/turbulentDispersionModels/Gosman/Gosman.H
index be62e78629cb56c3d37d86aa8cb5b620de233f0e..7de9411f56c53a5799ecc27662203e877f47c133 100644
--- a/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/turbulentDispersionModels/Gosman/Gosman.H
+++ b/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/turbulentDispersionModels/Gosman/Gosman.H
@@ -92,8 +92,9 @@ public:
 
     // Member Functions
 
-        //- Turbulent dispersion force
-        virtual tmp<volVectorField> F() const;
+        //- Turbulent dispersion force coefficient
+        //  multiplying the gradient of the phase-fraction
+        virtual tmp<volScalarField> Fprime() const;
 };
 
 
diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/turbulentDispersionModels/LopezDeBertodano/LopezDeBertodano.C b/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/turbulentDispersionModels/LopezDeBertodano/LopezDeBertodano.C
index 6debe6931e1f31b6c5084ea1a9edb32875d08bb2..419c25ca0ded5e200b65c4a492cd0a18ae67f392 100644
--- a/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/turbulentDispersionModels/LopezDeBertodano/LopezDeBertodano.C
+++ b/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/turbulentDispersionModels/LopezDeBertodano/LopezDeBertodano.C
@@ -25,7 +25,6 @@ License
 
 #include "LopezDeBertodano.H"
 #include "phasePair.H"
-#include "fvc.H"
 #include "PhaseCompressibleTurbulenceModel.H"
 #include "addToRunTimeSelectionTable.H"
 
@@ -67,14 +66,13 @@ Foam::turbulentDispersionModels::LopezDeBertodano::~LopezDeBertodano()
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-Foam::tmp<Foam::volVectorField>
-Foam::turbulentDispersionModels::LopezDeBertodano::F() const
+Foam::tmp<Foam::volScalarField>
+Foam::turbulentDispersionModels::LopezDeBertodano::Fprime() const
 {
     return
         Ctd_
        *pair_.continuous().rho()
-       *pair_.continuous().turbulence().k()
-       *fvc::grad(pair_.dispersed());
+       *pair_.continuous().turbulence().k();
 }
 
 // ************************************************************************* //
diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/turbulentDispersionModels/LopezDeBertodano/LopezDeBertodano.H b/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/turbulentDispersionModels/LopezDeBertodano/LopezDeBertodano.H
index 31af43db93f6a4ef4707283ad46ca6b7397b2b97..6f83bf0389d6778c665f9f27fed25bd0bcd587b5 100644
--- a/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/turbulentDispersionModels/LopezDeBertodano/LopezDeBertodano.H
+++ b/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/turbulentDispersionModels/LopezDeBertodano/LopezDeBertodano.H
@@ -98,8 +98,9 @@ public:
 
     // Member Functions
 
-        //- Turbulent dispersion force
-        virtual tmp<volVectorField> F() const;
+        //- Turbulent dispersion force coefficient
+        //  multiplying the gradient of the phase-fraction
+        virtual tmp<volScalarField> Fprime() const;
 };
 
 
diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/turbulentDispersionModels/constantTurbulentDispersionCoefficient/constantTurbulentDispersionCoefficient.C b/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/turbulentDispersionModels/constantTurbulentDispersionCoefficient/constantTurbulentDispersionCoefficient.C
index ba202d84fef5b90aedf09a9447a0cb8461dcf105..e1754a72ce084207da00c7179b8cec181e03fb95 100644
--- a/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/turbulentDispersionModels/constantTurbulentDispersionCoefficient/constantTurbulentDispersionCoefficient.C
+++ b/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/turbulentDispersionModels/constantTurbulentDispersionCoefficient/constantTurbulentDispersionCoefficient.C
@@ -25,7 +25,6 @@ License
 
 #include "constantTurbulentDispersionCoefficient.H"
 #include "phasePair.H"
-#include "fvc.H"
 #include "PhaseCompressibleTurbulenceModel.H"
 #include "addToRunTimeSelectionTable.H"
 
@@ -69,16 +68,15 @@ Foam::turbulentDispersionModels::constantTurbulentDispersionCoefficient::
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-Foam::tmp<Foam::volVectorField>
+Foam::tmp<Foam::volScalarField>
 Foam::turbulentDispersionModels::constantTurbulentDispersionCoefficient::
-F() const
+Fprime() const
 {
     return
         Ctd_
        *pair_.dispersed()
        *pair_.continuous().rho()
-       *pair_.continuous().turbulence().k()
-       *fvc::grad(pair_.dispersed());
+       *pair_.continuous().turbulence().k();
 }
 
 
diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/turbulentDispersionModels/constantTurbulentDispersionCoefficient/constantTurbulentDispersionCoefficient.H b/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/turbulentDispersionModels/constantTurbulentDispersionCoefficient/constantTurbulentDispersionCoefficient.H
index fba46ea740b08f15840c83f864665014e1d6ef0e..0a4146fc15ffca1b441cb7b5bf0400df309e9ac3 100644
--- a/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/turbulentDispersionModels/constantTurbulentDispersionCoefficient/constantTurbulentDispersionCoefficient.H
+++ b/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/turbulentDispersionModels/constantTurbulentDispersionCoefficient/constantTurbulentDispersionCoefficient.H
@@ -83,8 +83,9 @@ public:
 
     // Member Functions
 
-        //- Turbulent dispersion force
-        virtual tmp<volVectorField> F() const;
+        //- Turbulent dispersion force coefficient
+        //  multiplying the gradient of the phase-fraction
+        virtual tmp<volScalarField> Fprime() const;
 };
 
 
diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/turbulentDispersionModels/noTurbulentDispersion/noTurbulentDispersion.C b/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/turbulentDispersionModels/noTurbulentDispersion/noTurbulentDispersion.C
index 99c8e8581ce8044a1e37a1ccc30ba6227eadc692..44a4c3504d015142b02edb3fc652d3f6a0b1a5e4 100644
--- a/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/turbulentDispersionModels/noTurbulentDispersion/noTurbulentDispersion.C
+++ b/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/turbulentDispersionModels/noTurbulentDispersion/noTurbulentDispersion.C
@@ -93,4 +93,32 @@ Foam::turbulentDispersionModels::noTurbulentDispersion::F() const
 }
 
 
+Foam::tmp<Foam::volScalarField>
+Foam::turbulentDispersionModels::noTurbulentDispersion::Fprime() const
+{
+    const fvMesh& mesh(this->pair_.phase1().mesh());
+
+    return
+        tmp<volScalarField>
+        (
+            new volScalarField
+            (
+                IOobject
+                (
+                    "zero",
+                    mesh.time().timeName(),
+                    mesh
+                ),
+                mesh,
+                dimensionedScalar
+                (
+                    "zero",
+                    dimensionSet(1, -2, 1, 0, 0),
+                    0
+                )
+            )
+        );
+}
+
+
 // ************************************************************************* //
diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/turbulentDispersionModels/noTurbulentDispersion/noTurbulentDispersion.H b/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/turbulentDispersionModels/noTurbulentDispersion/noTurbulentDispersion.H
index 19c12a237ce2c802b3e718573c33275544933679..2153ff6866f2f2de203efedb5ebd2c79e56d38c9 100644
--- a/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/turbulentDispersionModels/noTurbulentDispersion/noTurbulentDispersion.H
+++ b/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/turbulentDispersionModels/noTurbulentDispersion/noTurbulentDispersion.H
@@ -78,6 +78,10 @@ public:
 
         //- Turbulent dispersion force
         virtual tmp<volVectorField> F() const;
+
+        //- Turbulent dispersion force coefficient
+        //  multiplying the gradient of the phase-fraction
+        virtual tmp<volScalarField> Fprime() const;
 };
 
 
diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/turbulentDispersionModels/turbulentDispersionModel/turbulentDispersionModel.C b/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/turbulentDispersionModels/turbulentDispersionModel/turbulentDispersionModel.C
index 7d11722babecdf90c4fecd78b87f8841352815d6..d4199b440f79603f6520b4dd4d07eb3dcbe48ebc 100644
--- a/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/turbulentDispersionModels/turbulentDispersionModel/turbulentDispersionModel.C
+++ b/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/turbulentDispersionModels/turbulentDispersionModel/turbulentDispersionModel.C
@@ -25,6 +25,7 @@ License
 
 #include "turbulentDispersionModel.H"
 #include "phasePair.H"
+#include "fvcGrad.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -55,4 +56,14 @@ Foam::turbulentDispersionModel::~turbulentDispersionModel()
 {}
 
 
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+Foam::tmp<Foam::volVectorField>
+Foam::turbulentDispersionModel::F() const
+{
+    return Fprime()*fvc::grad(pair_.dispersed());
+
+}
+
+
 // ************************************************************************* //
diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/turbulentDispersionModels/turbulentDispersionModel/turbulentDispersionModel.H b/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/turbulentDispersionModels/turbulentDispersionModel/turbulentDispersionModel.H
index b68903347690f42b1b1c93e52747f20ecb3ad4c7..5862452544a1c6d9bf5539860d77b5a3d36689af 100644
--- a/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/turbulentDispersionModels/turbulentDispersionModel/turbulentDispersionModel.H
+++ b/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/turbulentDispersionModels/turbulentDispersionModel/turbulentDispersionModel.H
@@ -82,7 +82,7 @@ public:
 
     // Static data members
 
-        //- Force dimensions 
+        //- Force dimensions
         static const dimensionSet dimF;
 
 
@@ -112,7 +112,11 @@ public:
     // Member Functions
 
         //- Turbulent dispersion force
-        virtual tmp<volVectorField> F() const = 0;
+        virtual tmp<volVectorField> F() const;
+
+        //- Turbulent dispersion force coefficient
+        //  multiplying the gradient of the phase-fraction
+        virtual tmp<volScalarField> Fprime() const = 0;
 };