diff --git a/applications/solvers/multiphase/compressibleInterFoam/alphaEqns.H b/applications/solvers/multiphase/compressibleInterFoam/alphaEqns.H
index 41ed1cf9c0..51dad81e77 100644
--- a/applications/solvers/multiphase/compressibleInterFoam/alphaEqns.H
+++ b/applications/solvers/multiphase/compressibleInterFoam/alphaEqns.H
@@ -45,7 +45,7 @@
         }
 
 
-        surfaceScalarField phiAlpha1
+        surfaceScalarField alphaPhi1
         (
             fvc::flux
             (
@@ -66,7 +66,7 @@
             geometricOneField(),
             alpha1,
             phi,
-            phiAlpha1,
+            alphaPhi1,
             Sp,
             Su,
             1,
@@ -75,7 +75,7 @@
 
         surfaceScalarField rho1f(fvc::interpolate(rho1));
         surfaceScalarField rho2f(fvc::interpolate(rho2));
-        rhoPhi = phiAlpha1*(rho1f - rho2f) + phi*rho2f;
+        rhoPhi = alphaPhi1*(rho1f - rho2f) + phi*rho2f;
 
         alpha2 = scalar(1) - alpha1;
     }
diff --git a/applications/solvers/multiphase/compressibleMultiphaseInterFoam/multiphaseMixtureThermo/multiphaseMixtureThermo.C b/applications/solvers/multiphase/compressibleMultiphaseInterFoam/multiphaseMixtureThermo/multiphaseMixtureThermo.C
index 6eda25f9ff..7733cf1364 100644
--- a/applications/solvers/multiphase/compressibleMultiphaseInterFoam/multiphaseMixtureThermo/multiphaseMixtureThermo.C
+++ b/applications/solvers/multiphase/compressibleMultiphaseInterFoam/multiphaseMixtureThermo/multiphaseMixtureThermo.C
@@ -942,14 +942,14 @@ void Foam::multiphaseMixtureThermo::solveAlphas
     surfaceScalarField phic(mag(phi_/mesh_.magSf()));
     phic = min(cAlpha*phic, max(phic));
 
-    PtrList<surfaceScalarField> phiAlphaCorrs(phases_.size());
+    PtrList<surfaceScalarField> alphaPhiCorrs(phases_.size());
     int phasei = 0;
 
     forAllIter(PtrDictionary<phaseModel>, phases_, phase)
     {
         phaseModel& alpha = phase();
 
-        phiAlphaCorrs.set
+        alphaPhiCorrs.set
         (
             phasei,
             new surfaceScalarField
@@ -964,7 +964,7 @@ void Foam::multiphaseMixtureThermo::solveAlphas
             )
         );
 
-        surfaceScalarField& phiAlphaCorr = phiAlphaCorrs[phasei];
+        surfaceScalarField& alphaPhiCorr = alphaPhiCorrs[phasei];
 
         forAllIter(PtrDictionary<phaseModel>, phases_, phase2)
         {
@@ -974,7 +974,7 @@ void Foam::multiphaseMixtureThermo::solveAlphas
 
             surfaceScalarField phir(phic*nHatf(alpha, alpha2));
 
-            phiAlphaCorr += fvc::flux
+            alphaPhiCorr += fvc::flux
             (
                 -fvc::flux(-phir, alpha2, alpharScheme),
                 alpha,
@@ -988,7 +988,7 @@ void Foam::multiphaseMixtureThermo::solveAlphas
             geometricOneField(),
             alpha,
             phi_,
-            phiAlphaCorr,
+            alphaPhiCorr,
             zeroField(),
             zeroField(),
             1,
@@ -999,7 +999,7 @@ void Foam::multiphaseMixtureThermo::solveAlphas
         phasei++;
     }
 
-    MULES::limitSum(phiAlphaCorrs);
+    MULES::limitSum(alphaPhiCorrs);
 
     rhoPhi_ = dimensionedScalar("0", dimensionSet(1, 0, -1, 0, 0), 0);
 
@@ -1025,8 +1025,8 @@ void Foam::multiphaseMixtureThermo::solveAlphas
     {
         phaseModel& alpha = phase();
 
-        surfaceScalarField& phiAlpha = phiAlphaCorrs[phasei];
-        phiAlpha += upwind<scalar>(mesh_, phi_).flux(alpha);
+        surfaceScalarField& alphaPhi = alphaPhiCorrs[phasei];
+        alphaPhi += upwind<scalar>(mesh_, phi_).flux(alpha);
 
         volScalarField::DimensionedInternalField Sp
         (
@@ -1096,12 +1096,12 @@ void Foam::multiphaseMixtureThermo::solveAlphas
         (
             geometricOneField(),
             alpha,
-            phiAlpha,
+            alphaPhi,
             Sp,
             Su
         );
 
-        rhoPhi_ += fvc::interpolate(alpha.thermo().rho())*phiAlpha;
+        rhoPhi_ += fvc::interpolate(alpha.thermo().rho())*alphaPhi;
 
         Info<< alpha.name() << " volume fraction, min, max = "
             << alpha.weightedAverage(mesh_.V()).value()
diff --git a/applications/solvers/multiphase/driftFluxFoam/alphaEqn.H b/applications/solvers/multiphase/driftFluxFoam/alphaEqn.H
index 954b2a8b3a..68f027cd24 100644
--- a/applications/solvers/multiphase/driftFluxFoam/alphaEqn.H
+++ b/applications/solvers/multiphase/driftFluxFoam/alphaEqn.H
@@ -23,32 +23,32 @@
             << "  Max(" << alpha1.name() << ") = " << max(alpha1).value()
             << endl;
 
-        tmp<surfaceScalarField> tphiAlphaUD(alpha1Eqn.flux());
-        phiAlpha = tphiAlphaUD();
+        tmp<surfaceScalarField> talphaPhiUD(alpha1Eqn.flux());
+        alphaPhi = talphaPhiUD();
 
-        if (alphaApplyPrevCorr && tphiAlphaCorr0.valid())
+        if (alphaApplyPrevCorr && talphaPhiCorr0.valid())
         {
             Info<< "Applying the previous iteration correction flux" << endl;
 
             MULES::correct
             (
                 alpha1,
-                phiAlpha,
-                tphiAlphaCorr0(),
+                alphaPhi,
+                talphaPhiCorr0(),
                 mixture.alphaMax(),
                 0
             );
 
-            phiAlpha += tphiAlphaCorr0();
+            alphaPhi += talphaPhiCorr0();
         }
 
         // Cache the upwind-flux
-        tphiAlphaCorr0 = tphiAlphaUD;
+        talphaPhiCorr0 = talphaPhiUD;
     }
 
     for (int aCorr=0; aCorr<nAlphaCorr; aCorr++)
     {
-        tmp<surfaceScalarField> tphiAlphaUn
+        tmp<surfaceScalarField> talphaPhiUn
         (
             fvc::flux
             (
@@ -66,14 +66,14 @@
 
         if (MULESCorr)
         {
-            tmp<surfaceScalarField> tphiAlphaCorr(tphiAlphaUn() - phiAlpha);
+            tmp<surfaceScalarField> talphaPhiCorr(talphaPhiUn() - alphaPhi);
             volScalarField alpha10("alpha10", alpha1);
 
             MULES::correct
             (
                 alpha1,
-                tphiAlphaUn(),
-                tphiAlphaCorr(),
+                talphaPhiUn(),
+                talphaPhiCorr(),
                 mixture.alphaMax(),
                 0
             );
@@ -81,23 +81,23 @@
             // Under-relax the correction for all but the 1st corrector
             if (aCorr == 0)
             {
-                phiAlpha += tphiAlphaCorr();
+                alphaPhi += talphaPhiCorr();
             }
             else
             {
                 alpha1 = 0.5*alpha1 + 0.5*alpha10;
-                phiAlpha += 0.5*tphiAlphaCorr();
+                alphaPhi += 0.5*talphaPhiCorr();
             }
         }
         else
         {
-            phiAlpha = tphiAlphaUn;
+            alphaPhi = talphaPhiUn;
 
             MULES::explicitSolve
             (
                 alpha1,
                 phi,
-                phiAlpha,
+                alphaPhi,
                 mixture.alphaMax(),
                 0
             );
@@ -106,7 +106,7 @@
 
     if (alphaApplyPrevCorr && MULESCorr)
     {
-        tphiAlphaCorr0 = phiAlpha - tphiAlphaCorr0;
+        talphaPhiCorr0 = alphaPhi - talphaPhiCorr0;
     }
 
     alpha2 = 1.0 - alpha1;
diff --git a/applications/solvers/multiphase/driftFluxFoam/alphaEqnSubCycle.H b/applications/solvers/multiphase/driftFluxFoam/alphaEqnSubCycle.H
index 956bacf48d..b586e77fad 100644
--- a/applications/solvers/multiphase/driftFluxFoam/alphaEqnSubCycle.H
+++ b/applications/solvers/multiphase/driftFluxFoam/alphaEqnSubCycle.H
@@ -1,9 +1,9 @@
 {
-    surfaceScalarField phiAlpha
+    surfaceScalarField alphaPhi
     (
         IOobject
         (
-            "phiAlpha",
+            "alphaPhi",
             runTime.timeName(),
             mesh
         ),
@@ -19,11 +19,11 @@
     if (nAlphaSubCycles > 1)
     {
         dimensionedScalar totalDeltaT = runTime.deltaT();
-        surfaceScalarField phiAlphaSum
+        surfaceScalarField alphaPhiSum
         (
             IOobject
             (
-                "phiAlphaSum",
+                "alphaPhiSum",
                 runTime.timeName(),
                 mesh
             ),
@@ -38,10 +38,10 @@
         )
         {
             #include "alphaEqn.H"
-            phiAlphaSum += (runTime.deltaT()/totalDeltaT)*phiAlpha;
+            alphaPhiSum += (runTime.deltaT()/totalDeltaT)*alphaPhi;
         }
 
-        phiAlpha = phiAlphaSum;
+        alphaPhi = alphaPhiSum;
     }
     else
     {
@@ -59,7 +59,7 @@
 
         alpha1Eqn.solve(mesh.solver("alpha1Diffusion"));
 
-        phiAlpha += alpha1Eqn.flux();
+        alphaPhi += alpha1Eqn.flux();
         alpha2 = 1.0 - alpha1;
 
         Info<< "Phase-1 volume fraction = "
@@ -69,6 +69,6 @@
             << endl;
     }
 
-    rhoPhi = phiAlpha*(rho1 - rho2) + phi*rho2;
+    rhoPhi = alphaPhi*(rho1 - rho2) + phi*rho2;
     rho = mixture.rho();
 }
diff --git a/applications/solvers/multiphase/driftFluxFoam/createFields.H b/applications/solvers/multiphase/driftFluxFoam/createFields.H
index d2e8d1d93e..65d3d38eae 100644
--- a/applications/solvers/multiphase/driftFluxFoam/createFields.H
+++ b/applications/solvers/multiphase/driftFluxFoam/createFields.H
@@ -135,4 +135,4 @@ mesh.setFluxRequired(p_rgh.name());
 mesh.setFluxRequired(alpha1.name());
 
 // MULES Correction
-tmp<surfaceScalarField> tphiAlphaCorr0;
+tmp<surfaceScalarField> talphaPhiCorr0;
diff --git a/applications/solvers/multiphase/interFoam/alphaEqn.H b/applications/solvers/multiphase/interFoam/alphaEqn.H
index 4405c58015..a2e9097e30 100644
--- a/applications/solvers/multiphase/interFoam/alphaEqn.H
+++ b/applications/solvers/multiphase/interFoam/alphaEqn.H
@@ -98,19 +98,19 @@
             << "  Max(" << alpha1.name() << ") = " << max(alpha1).value()
             << endl;
 
-        tmp<surfaceScalarField> tphiAlphaUD(alpha1Eqn.flux());
-        phiAlpha = tphiAlphaUD();
+        tmp<surfaceScalarField> talphaPhiUD(alpha1Eqn.flux());
+        alphaPhi = talphaPhiUD();
 
-        if (alphaApplyPrevCorr && tphiAlphaCorr0.valid())
+        if (alphaApplyPrevCorr && talphaPhiCorr0.valid())
         {
             Info<< "Applying the previous iteration compression flux" << endl;
-            MULES::correct(alpha1, phiAlpha, tphiAlphaCorr0(), 1, 0);
+            MULES::correct(alpha1, alphaPhi, talphaPhiCorr0(), 1, 0);
 
-            phiAlpha += tphiAlphaCorr0();
+            alphaPhi += talphaPhiCorr0();
         }
 
         // Cache the upwind-flux
-        tphiAlphaCorr0 = tphiAlphaUD;
+        talphaPhiCorr0 = talphaPhiUD;
 
         alpha2 = 1.0 - alpha1;
 
@@ -122,7 +122,7 @@
     {
         surfaceScalarField phir(phic*mixture.nHatf());
 
-        tmp<surfaceScalarField> tphiAlphaUn
+        tmp<surfaceScalarField> talphaPhiUn
         (
             fvc::flux
             (
@@ -141,33 +141,33 @@
         // Calculate the Crank-Nicolson off-centred alpha flux
         if (ocCoeff > 0)
         {
-            tphiAlphaUn =
-                cnCoeff*tphiAlphaUn + (1.0 - cnCoeff)*phiAlpha.oldTime();
+            talphaPhiUn =
+                cnCoeff*talphaPhiUn + (1.0 - cnCoeff)*alphaPhi.oldTime();
         }
 
         if (MULESCorr)
         {
-            tmp<surfaceScalarField> tphiAlphaCorr(tphiAlphaUn() - phiAlpha);
+            tmp<surfaceScalarField> talphaPhiCorr(talphaPhiUn() - alphaPhi);
             volScalarField alpha10("alpha10", alpha1);
 
-            MULES::correct(alpha1, tphiAlphaUn(), tphiAlphaCorr(), 1, 0);
+            MULES::correct(alpha1, talphaPhiUn(), talphaPhiCorr(), 1, 0);
 
             // Under-relax the correction for all but the 1st corrector
             if (aCorr == 0)
             {
-                phiAlpha += tphiAlphaCorr();
+                alphaPhi += talphaPhiCorr();
             }
             else
             {
                 alpha1 = 0.5*alpha1 + 0.5*alpha10;
-                phiAlpha += 0.5*tphiAlphaCorr();
+                alphaPhi += 0.5*talphaPhiCorr();
             }
         }
         else
         {
-            phiAlpha = tphiAlphaUn;
+            alphaPhi = talphaPhiUn;
 
-            MULES::explicitSolve(alpha1, phiCN, phiAlpha, 1, 0);
+            MULES::explicitSolve(alpha1, phiCN, alphaPhi, 1, 0);
         }
 
         alpha2 = 1.0 - alpha1;
@@ -177,7 +177,7 @@
 
     if (alphaApplyPrevCorr && MULESCorr)
     {
-        tphiAlphaCorr0 = phiAlpha - tphiAlphaCorr0;
+        talphaPhiCorr0 = alphaPhi - talphaPhiCorr0;
     }
 
     if
@@ -186,18 +186,18 @@
      == fv::EulerDdtScheme<vector>::typeName
     )
     {
-        rhoPhi = phiAlpha*(rho1 - rho2) + phiCN*rho2;
+        rhoPhi = alphaPhi*(rho1 - rho2) + phiCN*rho2;
     }
     else
     {
         if (ocCoeff > 0)
         {
             // Calculate the end-of-time-step alpha flux
-            phiAlpha = (phiAlpha - (1.0 - cnCoeff)*phiAlpha.oldTime())/cnCoeff;
+            alphaPhi = (alphaPhi - (1.0 - cnCoeff)*alphaPhi.oldTime())/cnCoeff;
         }
 
         // Calculate the end-of-time-step mass flux
-        rhoPhi = phiAlpha*(rho1 - rho2) + phi*rho2;
+        rhoPhi = alphaPhi*(rho1 - rho2) + phi*rho2;
     }
 
     Info<< "Phase-1 volume fraction = "
diff --git a/applications/solvers/multiphase/interFoam/createFields.H b/applications/solvers/multiphase/interFoam/createFields.H
index fc682190ce..c3e65f24e3 100644
--- a/applications/solvers/multiphase/interFoam/createFields.H
+++ b/applications/solvers/multiphase/interFoam/createFields.H
@@ -121,11 +121,11 @@ mesh.setFluxRequired(p_rgh.name());
 mesh.setFluxRequired(alpha1.name());
 
 // MULES flux from previous time-step
-surfaceScalarField phiAlpha
+surfaceScalarField alphaPhi
 (
     IOobject
     (
-        "phiAlpha",
+        "alphaPhi",
         runTime.timeName(),
         mesh,
         IOobject::READ_IF_PRESENT,
@@ -135,4 +135,4 @@ surfaceScalarField phiAlpha
 );
 
 // MULES Correction
-tmp<surfaceScalarField> tphiAlphaCorr0;
+tmp<surfaceScalarField> talphaPhiCorr0;
diff --git a/applications/solvers/multiphase/interFoam/interMixingFoam/alphaEqns.H b/applications/solvers/multiphase/interFoam/interMixingFoam/alphaEqns.H
index 9585905c08..d720aa97f6 100644
--- a/applications/solvers/multiphase/interFoam/interMixingFoam/alphaEqns.H
+++ b/applications/solvers/multiphase/interFoam/interMixingFoam/alphaEqns.H
@@ -40,7 +40,7 @@
 
 
         // Create the complete convection flux for alpha1
-        surfaceScalarField phiAlpha1
+        surfaceScalarField alphaPhi1
         (
             fvc::flux
             (
@@ -63,13 +63,13 @@
         );
 
         // Create the bounded (upwind) flux for alpha1
-        surfaceScalarField phiAlpha1BD
+        surfaceScalarField alphaPhi1BD
         (
             upwind<scalar>(mesh, phi).flux(alpha1)
         );
 
         // Calculate the flux correction for alpha1
-        phiAlpha1 -= phiAlpha1BD;
+        alphaPhi1 -= alphaPhi1BD;
 
         // Calculate the limiter for alpha1
         if (LTS)
@@ -83,8 +83,8 @@
                 rDeltaT,
                 geometricOneField(),
                 alpha1,
-                phiAlpha1BD,
-                phiAlpha1,
+                alphaPhi1BD,
+                alphaPhi1,
                 zeroField(),
                 zeroField(),
                 1,
@@ -99,8 +99,8 @@
                 1.0/runTime.deltaT().value(),
                 geometricOneField(),
                 alpha1,
-                phiAlpha1BD,
-                phiAlpha1,
+                alphaPhi1BD,
+                alphaPhi1,
                 zeroField(),
                 zeroField(),
                 1,
@@ -109,7 +109,7 @@
         }
 
         // Create the complete flux for alpha2
-        surfaceScalarField phiAlpha2
+        surfaceScalarField alphaPhi2
         (
             fvc::flux
             (
@@ -126,13 +126,13 @@
         );
 
         // Create the bounded (upwind) flux for alpha2
-        surfaceScalarField phiAlpha2BD
+        surfaceScalarField alphaPhi2BD
         (
             upwind<scalar>(mesh, phi).flux(alpha2)
         );
 
         // Calculate the flux correction for alpha2
-        phiAlpha2 -= phiAlpha2BD;
+        alphaPhi2 -= alphaPhi2BD;
 
         // Further limit the limiter for alpha2
         if (LTS)
@@ -146,8 +146,8 @@
                 rDeltaT,
                 geometricOneField(),
                 alpha2,
-                phiAlpha2BD,
-                phiAlpha2,
+                alphaPhi2BD,
+                alphaPhi2,
                 zeroField(),
                 zeroField(),
                 1,
@@ -162,8 +162,8 @@
                 1.0/runTime.deltaT().value(),
                 geometricOneField(),
                 alpha2,
-                phiAlpha2BD,
-                phiAlpha2,
+                alphaPhi2BD,
+                alphaPhi2,
                 zeroField(),
                 zeroField(),
                 1,
@@ -172,32 +172,32 @@
         }
 
         // Construct the limited fluxes
-        phiAlpha1 = phiAlpha1BD + lambda*phiAlpha1;
-        phiAlpha2 = phiAlpha2BD + lambda*phiAlpha2;
+        alphaPhi1 = alphaPhi1BD + lambda*alphaPhi1;
+        alphaPhi2 = alphaPhi2BD + lambda*alphaPhi2;
 
         // Solve for alpha1
-        solve(fvm::ddt(alpha1) + fvc::div(phiAlpha1));
+        solve(fvm::ddt(alpha1) + fvc::div(alphaPhi1));
 
         // Create the diffusion coefficients for alpha2<->alpha3
         volScalarField Dc23(D23*max(alpha3, scalar(0))*pos(alpha2));
         volScalarField Dc32(D23*max(alpha2, scalar(0))*pos(alpha3));
 
         // Add the diffusive flux for alpha3->alpha2
-        phiAlpha2 -= fvc::interpolate(Dc32)*mesh.magSf()*fvc::snGrad(alpha1);
+        alphaPhi2 -= fvc::interpolate(Dc32)*mesh.magSf()*fvc::snGrad(alpha1);
 
         // Solve for alpha2
         fvScalarMatrix alpha2Eqn
         (
             fvm::ddt(alpha2)
-          + fvc::div(phiAlpha2)
+          + fvc::div(alphaPhi2)
           - fvm::laplacian(Dc23 + Dc32, alpha2)
         );
         alpha2Eqn.solve();
 
         // Construct the complete mass flux
         rhoPhi =
-              phiAlpha1*(rho1 - rho3)
-            + (phiAlpha2 + alpha2Eqn.flux())*(rho2 - rho3)
+              alphaPhi1*(rho1 - rho3)
+            + (alphaPhi2 + alpha2Eqn.flux())*(rho2 - rho3)
             + phi*rho3;
 
         alpha3 = 1.0 - alpha1 - alpha2;
diff --git a/applications/solvers/multiphase/interPhaseChangeFoam/alphaEqn.H b/applications/solvers/multiphase/interPhaseChangeFoam/alphaEqn.H
index bf506c2a70..172676b97f 100644
--- a/applications/solvers/multiphase/interPhaseChangeFoam/alphaEqn.H
+++ b/applications/solvers/multiphase/interPhaseChangeFoam/alphaEqn.H
@@ -10,7 +10,7 @@
     const volScalarField& vDotvAlphal = vDotAlphal[1]();
     const volScalarField vDotvmcAlphal(vDotvAlphal - vDotcAlphal);
 
-    tmp<surfaceScalarField> tphiAlpha;
+    tmp<surfaceScalarField> talphaPhi;
 
     if (MULESCorr)
     {
@@ -37,14 +37,14 @@
             << "  Max(" << alpha1.name() << ") = " << max(alpha1).value()
             << endl;
 
-        tphiAlpha = alpha1Eqn.flux();
+        talphaPhi = alpha1Eqn.flux();
     }
 
     volScalarField alpha10("alpha10", alpha1);
 
     for (int aCorr=0; aCorr<nAlphaCorr; aCorr++)
     {
-        tmp<surfaceScalarField> tphiAlphaCorr
+        tmp<surfaceScalarField> talphaPhiCorr
         (
             fvc::flux
             (
@@ -62,7 +62,7 @@
 
         if (MULESCorr)
         {
-            tphiAlphaCorr() -= tphiAlpha();
+            talphaPhiCorr() -= talphaPhi();
 
             volScalarField alpha100("alpha100", alpha10);
             alpha10 = alpha1;
@@ -71,8 +71,8 @@
             (
                 geometricOneField(),
                 alpha1,
-                tphiAlpha(),
-                tphiAlphaCorr(),
+                talphaPhi(),
+                talphaPhiCorr(),
                 vDotvmcAlphal,
                 (
                     divU*(alpha10 - alpha100)
@@ -85,12 +85,12 @@
             // Under-relax the correction for all but the 1st corrector
             if (aCorr == 0)
             {
-                tphiAlpha() += tphiAlphaCorr();
+                talphaPhi() += talphaPhiCorr();
             }
             else
             {
                 alpha1 = 0.5*alpha1 + 0.5*alpha10;
-                tphiAlpha() += 0.5*tphiAlphaCorr();
+                talphaPhi() += 0.5*talphaPhiCorr();
             }
         }
         else
@@ -100,20 +100,20 @@
                 geometricOneField(),
                 alpha1,
                 phi,
-                tphiAlphaCorr(),
+                talphaPhiCorr(),
                 vDotvmcAlphal,
                 (divU*alpha1 + vDotcAlphal)(),
                 1,
                 0
             );
 
-            tphiAlpha = tphiAlphaCorr;
+            talphaPhi = talphaPhiCorr;
         }
 
         alpha2 = 1.0 - alpha1;
     }
 
-    rhoPhi = tphiAlpha()*(rho1 - rho2) + phi*rho2;
+    rhoPhi = talphaPhi()*(rho1 - rho2) + phi*rho2;
 
     Info<< "Liquid phase volume fraction = "
         << alpha1.weightedAverage(mesh.V()).value()
diff --git a/applications/solvers/multiphase/multiphaseInterFoam/multiphaseMixture/multiphaseMixture.C b/applications/solvers/multiphase/multiphaseInterFoam/multiphaseMixture/multiphaseMixture.C
index 16d1b5f809..f1305897ba 100644
--- a/applications/solvers/multiphase/multiphaseInterFoam/multiphaseMixture/multiphaseMixture.C
+++ b/applications/solvers/multiphase/multiphaseInterFoam/multiphaseMixture/multiphaseMixture.C
@@ -573,14 +573,14 @@ void Foam::multiphaseMixture::solveAlphas
     surfaceScalarField phic(mag(phi_/mesh_.magSf()));
     phic = min(cAlpha*phic, max(phic));
 
-    PtrList<surfaceScalarField> phiAlphaCorrs(phases_.size());
+    PtrList<surfaceScalarField> alphaPhiCorrs(phases_.size());
     int phasei = 0;
 
     forAllIter(PtrDictionary<phase>, phases_, iter)
     {
         phase& alpha = iter();
 
-        phiAlphaCorrs.set
+        alphaPhiCorrs.set
         (
             phasei,
             new surfaceScalarField
@@ -595,7 +595,7 @@ void Foam::multiphaseMixture::solveAlphas
             )
         );
 
-        surfaceScalarField& phiAlphaCorr = phiAlphaCorrs[phasei];
+        surfaceScalarField& alphaPhiCorr = alphaPhiCorrs[phasei];
 
         forAllIter(PtrDictionary<phase>, phases_, iter2)
         {
@@ -605,7 +605,7 @@ void Foam::multiphaseMixture::solveAlphas
 
             surfaceScalarField phir(phic*nHatf(alpha, alpha2));
 
-            phiAlphaCorr += fvc::flux
+            alphaPhiCorr += fvc::flux
             (
                 -fvc::flux(-phir, alpha2, alpharScheme),
                 alpha,
@@ -619,7 +619,7 @@ void Foam::multiphaseMixture::solveAlphas
             geometricOneField(),
             alpha,
             phi_,
-            phiAlphaCorr,
+            alphaPhiCorr,
             zeroField(),
             zeroField(),
             1,
@@ -630,7 +630,7 @@ void Foam::multiphaseMixture::solveAlphas
         phasei++;
     }
 
-    MULES::limitSum(phiAlphaCorrs);
+    MULES::limitSum(alphaPhiCorrs);
 
     rhoPhi_ = dimensionedScalar("0", dimensionSet(1, 0, -1, 0, 0), 0);
 
@@ -652,19 +652,19 @@ void Foam::multiphaseMixture::solveAlphas
     {
         phase& alpha = iter();
 
-        surfaceScalarField& phiAlpha = phiAlphaCorrs[phasei];
-        phiAlpha += upwind<scalar>(mesh_, phi_).flux(alpha);
+        surfaceScalarField& alphaPhi = alphaPhiCorrs[phasei];
+        alphaPhi += upwind<scalar>(mesh_, phi_).flux(alpha);
 
         MULES::explicitSolve
         (
             geometricOneField(),
             alpha,
-            phiAlpha,
+            alphaPhi,
             zeroField(),
             zeroField()
         );
 
-        rhoPhi_ += phiAlpha*alpha.rho();
+        rhoPhi_ += alphaPhi*alpha.rho();
 
         Info<< alpha.name() << " volume fraction, min, max = "
             << alpha.weightedAverage(mesh_.V()).value()
diff --git a/applications/solvers/multiphase/twoLiquidMixingFoam/alphaEqn.H b/applications/solvers/multiphase/twoLiquidMixingFoam/alphaEqn.H
index bd10646871..5344fb75f2 100644
--- a/applications/solvers/multiphase/twoLiquidMixingFoam/alphaEqn.H
+++ b/applications/solvers/multiphase/twoLiquidMixingFoam/alphaEqn.H
@@ -1,7 +1,7 @@
 {
     word alphaScheme("div(phi,alpha)");
 
-    surfaceScalarField phiAlpha
+    surfaceScalarField alphaPhi
     (
         phi.name() + alpha1.name(),
         fvc::flux
@@ -12,9 +12,9 @@
         )
     );
 
-    MULES::explicitSolve(alpha1, phi, phiAlpha, 1, 0);
+    MULES::explicitSolve(alpha1, phi, alphaPhi, 1, 0);
 
-    rhoPhi = phiAlpha*(rho1 - rho2) + phi*rho2;
+    rhoPhi = alphaPhi*(rho1 - rho2) + phi*rho2;
 
     Info<< "Phase-1 volume fraction = "
         << alpha1.weightedAverage(mesh.Vsc()).value()