diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/turbulentDispersionModels/Burns/Burns.C b/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/turbulentDispersionModels/Burns/Burns.C
index 7fa406498..1de7c87de 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 edd33dc59..0d414bd5f 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 7f7cc5652..005ef7c1a 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 be62e7862..7de9411f5 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 6debe6931..419c25ca0 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 31af43db9..6f83bf038 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 ba202d84f..e1754a72c 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 fba46ea74..0a4146fc1 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 99c8e8581..44a4c3504 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 19c12a237..2153ff686 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 7d11722ba..d4199b440 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 b68903347..586245254 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;
 };