diff --git a/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/pU/pEqn.H b/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/pU/pEqn.H
index 9a4174e8c13b7636e6b07ead46e30c2bee86c18d..338b0d88b67ffc8a2352d12ae98ea1601f52d684 100644
--- a/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/pU/pEqn.H
+++ b/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/pU/pEqn.H
@@ -338,14 +338,14 @@ while (pimple.correct())
                 phi2 = phi - alphaf1*phir;
             }
 
-            // Compressibility correction for phase-fraction equations
+            // Set the phase dilatation rates
             if (phase1.compressible())
             {
-                phase1.D(pEqnComp1 & p_rgh);
+                phase1.divU(-pEqnComp1 & p_rgh);
             }
             if (phase2.compressible())
             {
-                phase2.D(pEqnComp2 & p_rgh);
+                phase2.divU(-pEqnComp2 & p_rgh);
             }
 
             // Optionally relax pressure for velocity correction
diff --git a/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/pUf/pEqn.H b/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/pUf/pEqn.H
index 05ed9e7d087616929a066c304209c0434d4791a1..8736d0d5a3fba16bc0605cb7b58837878c10d38c 100644
--- a/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/pUf/pEqn.H
+++ b/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/pUf/pEqn.H
@@ -323,14 +323,14 @@ while (pimple.correct())
             U2.correctBoundaryConditions();
             fvOptions.correct(U2);
 
-            // Compressibility correction for phase-fraction equations
+            // Set the phase dilatation rates
             if (phase1.compressible())
             {
-                phase1.D(pEqnComp1 & p_rgh);
+                phase1.divU(-pEqnComp1 & p_rgh);
             }
             if (phase2.compressible())
             {
-                phase2.D(pEqnComp2 & p_rgh);
+                phase2.divU(-pEqnComp2 & p_rgh);
             }
         }
     }
diff --git a/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/phaseModel/AnisothermalPhaseModel/AnisothermalPhaseModel.C b/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/phaseModel/AnisothermalPhaseModel/AnisothermalPhaseModel.C
index 99e16af8a0537ee0cf127498d939fd0519fd0285..01ddeb772f82c2dcaea61b8907ad5808edba4be5 100644
--- a/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/phaseModel/AnisothermalPhaseModel/AnisothermalPhaseModel.C
+++ b/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/phaseModel/AnisothermalPhaseModel/AnisothermalPhaseModel.C
@@ -36,16 +36,16 @@ Foam::AnisothermalPhaseModel<BasePhaseModel>::AnisothermalPhaseModel
 )
 :
     BasePhaseModel(fluid, phaseName),
-    D_
+    divU_
     (
         IOobject
         (
-            IOobject::groupName("D", this->name()),
+            IOobject::groupName("divU", this->name()),
             fluid.mesh().time().timeName(),
             fluid.mesh()
         ),
         fluid.mesh(),
-        dimensionedScalar("D", dimless/dimTime, 0)
+        dimensionedScalar("divU", dimless/dimTime, 0)
     ),
     K_
     (
@@ -140,17 +140,17 @@ bool Foam::AnisothermalPhaseModel<BasePhaseModel>::compressible() const
 
 template<class BasePhaseModel>
 Foam::tmp<Foam::volScalarField>
-Foam::AnisothermalPhaseModel<BasePhaseModel>::D() const
+Foam::AnisothermalPhaseModel<BasePhaseModel>::divU() const
 {
-    return D_;
+    return divU_;
 }
 
 
 template<class BasePhaseModel>
 void
-Foam::AnisothermalPhaseModel<BasePhaseModel>::D(const volScalarField& D)
+Foam::AnisothermalPhaseModel<BasePhaseModel>::divU(const volScalarField& divU)
 {
-    D_ = D;
+    divU_ = divU;
 }
 
 
diff --git a/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/phaseModel/AnisothermalPhaseModel/AnisothermalPhaseModel.H b/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/phaseModel/AnisothermalPhaseModel/AnisothermalPhaseModel.H
index c69156b356d17b0ac7863dd4be43865d7fd6dca4..1ff8b79d8ba32c2d47beffee6892305f1a324b52 100644
--- a/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/phaseModel/AnisothermalPhaseModel/AnisothermalPhaseModel.H
+++ b/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/phaseModel/AnisothermalPhaseModel/AnisothermalPhaseModel.H
@@ -54,8 +54,8 @@ class AnisothermalPhaseModel
 {
     // Private data
 
-        //- Dilatation
-        volScalarField D_;
+        //- Dilatation rate
+        volScalarField divU_;
 
         //- Kinetic energy
         volScalarField K_;
@@ -89,11 +89,11 @@ public:
             //- Return true if the phase is compressible otherwise false
             virtual bool compressible() const;
 
-            //- Phase dilatation rate ((alpha/rho)*Drho/Dt)
-            virtual tmp<volScalarField> D() const;
+            //- Phase dilatation rate (d(alpha)/dt + div(alpha*phi))
+            virtual tmp<volScalarField> divU() const;
 
-            //- Set phase dilatation rate ((alpha/rho)*Drho/Dt)
-            virtual void D(const volScalarField& D);
+            //- Set phase dilatation rate (d(alpha)/dt + div(alpha*phi))
+            virtual void divU(const volScalarField& divU);
 };
 
 
diff --git a/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/phaseModel/phaseModel/phaseModel.C b/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/phaseModel/phaseModel/phaseModel.C
index d1cc94e6c193519f818da71daa2b04b6b14b6d76..5f290080d4c6022b5ad8d1c746a48d610ca1ab6e 100644
--- a/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/phaseModel/phaseModel/phaseModel.C
+++ b/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/phaseModel/phaseModel/phaseModel.C
@@ -140,15 +140,15 @@ bool Foam::phaseModel::compressible() const
 }
 
 
-Foam::tmp<Foam::volScalarField> Foam::phaseModel::D() const
+Foam::tmp<Foam::volScalarField> Foam::phaseModel::divU() const
 {
     return tmp<volScalarField>();
 }
 
 
-void Foam::phaseModel::D(const volScalarField& D)
+void Foam::phaseModel::divU(const volScalarField& divU)
 {
-    WarningIn("phaseModel::D(const volScalarField& D)")
+    WarningIn("phaseModel::divU(const volScalarField& divU)")
         << "Attempt to set the dilatation rate of an incompressible phase"
         << endl;
 }
diff --git a/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/phaseModel/phaseModel/phaseModel.H b/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/phaseModel/phaseModel/phaseModel.H
index e42a3a70b441cddf1713904a5e7f3f44bb22942e..4fe4b9d1aa0283705dc7d32004f4b79bf43efcbd 100644
--- a/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/phaseModel/phaseModel/phaseModel.H
+++ b/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/phaseModel/phaseModel/phaseModel.H
@@ -169,11 +169,11 @@ public:
             //- Return true if the phase is compressible otherwise false
             virtual bool compressible() const;
 
-            //- Phase dilatation rate ((alpha/rho)*Drho/Dt)
-            virtual tmp<volScalarField> D() const;
+            //- Phase dilatation rate (d(alpha)/dt + div(alpha*phi))
+            virtual tmp<volScalarField> divU() const;
 
-            //- Set phase dilatation rate ((alpha/rho)*Drho/Dt)
-            virtual void D(const volScalarField& D);
+            //- Set phase dilatation rate (d(alpha)/dt + div(alpha*phi))
+            virtual void divU(const volScalarField& divU);
 
 
         // Thermo
diff --git a/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/twoPhaseSystem/twoPhaseSystem.C b/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/twoPhaseSystem/twoPhaseSystem.C
index 00ebfcffb98faf5e534c8c93cf73c3a3942c94dc..825089539fd9eee0134a974658dc5a5f774cecdb 100644
--- a/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/twoPhaseSystem/twoPhaseSystem.C
+++ b/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/twoPhaseSystem/twoPhaseSystem.C
@@ -189,26 +189,26 @@ void Foam::twoPhaseSystem::solve()
     {
         tdgdt =
         (
-            alpha1.dimensionedInternalField()
-           *phase2_.D()().dimensionedInternalField()
-          - alpha2.dimensionedInternalField()
-           *phase1_.D()().dimensionedInternalField()
+            alpha2.dimensionedInternalField()
+           *phase1_.divU()().dimensionedInternalField()
+          - alpha1.dimensionedInternalField()
+           *phase2_.divU()().dimensionedInternalField()
         );
     }
     else if (phase1_.compressible())
     {
         tdgdt =
         (
-          - alpha2.dimensionedInternalField()
-           *phase1_.D()().dimensionedInternalField()
+            alpha2.dimensionedInternalField()
+           *phase1_.divU()().dimensionedInternalField()
         );
     }
     else if (phase2_.compressible())
     {
         tdgdt =
         (
-            alpha1.dimensionedInternalField()
-           *phase2_.D()().dimensionedInternalField()
+          - alpha1.dimensionedInternalField()
+           *phase2_.divU()().dimensionedInternalField()
         );
     }