diff --git a/applications/solvers/DNS/dnsFoam/dnsFoam.C b/applications/solvers/DNS/dnsFoam/dnsFoam.C
index eb1e69c41b..34f212bc4d 100644
--- a/applications/solvers/DNS/dnsFoam/dnsFoam.C
+++ b/applications/solvers/DNS/dnsFoam/dnsFoam.C
@@ -85,15 +85,15 @@ int main(int argc, char *argv[])
 
         for (int corr=1; corr<=1; corr++)
         {
-            volScalarField rUA = 1.0/UEqn.A();
+            volScalarField rAU = 1.0/UEqn.A();
 
-            U = rUA*UEqn.H();
+            U = rAU*UEqn.H();
             phi = (fvc::interpolate(U) & mesh.Sf())
-                + fvc::ddtPhiCorr(rUA, U, phi);
+                + fvc::ddtPhiCorr(rAU, U, phi);
 
             fvScalarMatrix pEqn
             (
-                fvm::laplacian(rUA, p) == fvc::div(phi)
+                fvm::laplacian(rAU, p) == fvc::div(phi)
             );
 
             pEqn.solve();
@@ -102,7 +102,7 @@ int main(int argc, char *argv[])
 
             #include "continuityErrs.H"
 
-            U -= rUA*fvc::grad(p);
+            U -= rAU*fvc::grad(p);
             U.correctBoundaryConditions();
         }
 
diff --git a/applications/solvers/combustion/PDRFoam/pEqn.H b/applications/solvers/combustion/PDRFoam/pEqn.H
index 524c8eac74..9042837f54 100644
--- a/applications/solvers/combustion/PDRFoam/pEqn.H
+++ b/applications/solvers/combustion/PDRFoam/pEqn.H
@@ -1,6 +1,6 @@
 rho = thermo.rho();
 
-volScalarField rUA = 1.0/UEqn.A();
+volScalarField rAU = 1.0/UEqn.A();
 U = invA & UEqn.H();
 
 if (transonic)
@@ -11,7 +11,7 @@ if (transonic)
         fvc::interpolate(psi)
        *(
             (fvc::interpolate(U) & mesh.Sf())
-          + fvc::ddtPhiCorr(rUA, rho, U, phi)
+          + fvc::ddtPhiCorr(rAU, rho, U, phi)
         )
     );
 
@@ -38,7 +38,7 @@ else
         fvc::interpolate(rho)*
         (
             (fvc::interpolate(U) & mesh.Sf())
-          + fvc::ddtPhiCorr(rUA, rho, U, phi)
+          + fvc::ddtPhiCorr(rAU, rho, U, phi)
         );
 
     for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
diff --git a/applications/solvers/combustion/XiFoam/pEqn.H b/applications/solvers/combustion/XiFoam/pEqn.H
index 9443f909a3..3d414f4af6 100644
--- a/applications/solvers/combustion/XiFoam/pEqn.H
+++ b/applications/solvers/combustion/XiFoam/pEqn.H
@@ -1,7 +1,7 @@
 rho = thermo.rho();
 
-volScalarField rUA = 1.0/UEqn.A();
-U = rUA*UEqn.H();
+volScalarField rAU = 1.0/UEqn.A();
+U = rAU*UEqn.H();
 
 if (transonic)
 {
@@ -11,7 +11,7 @@ if (transonic)
         fvc::interpolate(psi)
        *(
             (fvc::interpolate(U) & mesh.Sf())
-          + fvc::ddtPhiCorr(rUA, rho, U, phi)
+          + fvc::ddtPhiCorr(rAU, rho, U, phi)
         )
     );
 
@@ -21,7 +21,7 @@ if (transonic)
         (
             fvm::ddt(psi, p)
           + fvm::div(phid, p)
-          - fvm::laplacian(rho*rUA, p)
+          - fvm::laplacian(rho*rAU, p)
         );
 
         pEqn.solve();
@@ -38,7 +38,7 @@ else
         fvc::interpolate(rho)
        *(
             (fvc::interpolate(U) & mesh.Sf())
-          + fvc::ddtPhiCorr(rUA, rho, U, phi)
+          + fvc::ddtPhiCorr(rAU, rho, U, phi)
         );
 
     for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
@@ -47,7 +47,7 @@ else
         (
             fvm::ddt(psi, p)
           + fvc::div(phi)
-          - fvm::laplacian(rho*rUA, p)
+          - fvm::laplacian(rho*rAU, p)
         );
 
         pEqn.solve();
@@ -62,7 +62,7 @@ else
 #include "rhoEqn.H"
 #include "compressibleContinuityErrs.H"
 
-U -= rUA*fvc::grad(p);
+U -= rAU*fvc::grad(p);
 U.correctBoundaryConditions();
 
 DpDt = fvc::DDt(surfaceScalarField("phiU", phi/fvc::interpolate(rho)), p);
diff --git a/applications/solvers/combustion/dieselFoam/pEqn.H b/applications/solvers/combustion/dieselFoam/pEqn.H
index d74947fe53..f1b66877a0 100644
--- a/applications/solvers/combustion/dieselFoam/pEqn.H
+++ b/applications/solvers/combustion/dieselFoam/pEqn.H
@@ -1,7 +1,7 @@
 rho = thermo.rho();
 
-volScalarField rUA = 1.0/UEqn.A();
-U = rUA*UEqn.H();
+volScalarField rAU = 1.0/UEqn.A();
+U = rAU*UEqn.H();
 
 if (transonic)
 {
@@ -11,7 +11,7 @@ if (transonic)
         fvc::interpolate(psi)
        *(
             (fvc::interpolate(U) & mesh.Sf())
-          + fvc::ddtPhiCorr(rUA, rho, U, phi)
+          + fvc::ddtPhiCorr(rAU, rho, U, phi)
         )
     );
 
@@ -21,7 +21,7 @@ if (transonic)
         (
             fvm::ddt(psi, p)
           + fvm::div(phid, p)
-          - fvm::laplacian(rho*rUA, p)
+          - fvm::laplacian(rho*rAU, p)
          ==
             Sevap
         );
@@ -40,7 +40,7 @@ else
         fvc::interpolate(rho)
        *(
             (fvc::interpolate(U) & mesh.Sf())
-          + fvc::ddtPhiCorr(rUA, rho, U, phi)
+          + fvc::ddtPhiCorr(rAU, rho, U, phi)
         );
 
     for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
@@ -49,7 +49,7 @@ else
         (
             fvm::ddt(psi, p)
           + fvc::div(phi)
-          - fvm::laplacian(rho*rUA, p)
+          - fvm::laplacian(rho*rAU, p)
          ==
             Sevap
         );
@@ -66,7 +66,7 @@ else
 #include "rhoEqn.H"
 #include "compressibleContinuityErrs.H"
 
-U -= rUA*fvc::grad(p);
+U -= rAU*fvc::grad(p);
 U.correctBoundaryConditions();
 
 DpDt = fvc::DDt(surfaceScalarField("phiU", phi/fvc::interpolate(rho)), p);
diff --git a/applications/solvers/combustion/engineFoam/pEqn.H b/applications/solvers/combustion/engineFoam/pEqn.H
index 39b4967312..913180f794 100644
--- a/applications/solvers/combustion/engineFoam/pEqn.H
+++ b/applications/solvers/combustion/engineFoam/pEqn.H
@@ -1,7 +1,7 @@
 rho = thermo.rho();
 
-volScalarField rUA = 1.0/UEqn.A();
-U = rUA*UEqn.H();
+volScalarField rAU = 1.0/UEqn.A();
+U = rAU*UEqn.H();
 
 if (transonic)
 {
@@ -18,7 +18,7 @@ if (transonic)
         (
             fvm::ddt(psi, p)
           + fvm::div(phid, p, "div(phid,p)")
-          - fvm::laplacian(rho*rUA, p)
+          - fvm::laplacian(rho*rAU, p)
         );
 
         pEqn.solve();
@@ -40,7 +40,7 @@ else
         (
             fvm::ddt(psi, p)
           + fvc::div(phi)
-          - fvm::laplacian(rho*rUA, p)
+          - fvm::laplacian(rho*rAU, p)
         );
 
         pEqn.solve();
@@ -55,7 +55,7 @@ else
 #include "rhoEqn.H"
 #include "compressibleContinuityErrs.H"
 
-U -= rUA*fvc::grad(p);
+U -= rAU*fvc::grad(p);
 U.correctBoundaryConditions();
 
 DpDt = fvc::DDt(surfaceScalarField("phiU", phi/fvc::interpolate(rho)), p);
diff --git a/applications/solvers/combustion/fireFoam/pEqn.H b/applications/solvers/combustion/fireFoam/pEqn.H
index b31a7c6636..135b55d879 100644
--- a/applications/solvers/combustion/fireFoam/pEqn.H
+++ b/applications/solvers/combustion/fireFoam/pEqn.H
@@ -1,29 +1,29 @@
 rho = thermo.rho();
 
-volScalarField rUA = 1.0/UEqn.A();
-surfaceScalarField rhorUAf("(rho*(1|A(U)))", fvc::interpolate(rho*rUA));
-U = rUA*UEqn.H();
+volScalarField rAU = 1.0/UEqn.A();
+surfaceScalarField rhorAUf("(rho*(1|A(U)))", fvc::interpolate(rho*rAU));
+U = rAU*UEqn.H();
 
 surfaceScalarField phiU
 (
     fvc::interpolate(rho)
    *(
         (fvc::interpolate(U) & mesh.Sf())
-      + fvc::ddtPhiCorr(rUA, rho, U, phi)
+      + fvc::ddtPhiCorr(rAU, rho, U, phi)
     )
 );
 
-phi = phiU - rhorUAf*ghf*fvc::snGrad(rho)*mesh.magSf();
+phi = phiU - rhorAUf*ghf*fvc::snGrad(rho)*mesh.magSf();
 
 for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
 {
-    surfaceScalarField rhorUAf = fvc::interpolate(rho*rUA);
+    surfaceScalarField rhorAUf = fvc::interpolate(rho*rAU);
 
     fvScalarMatrix p_rghEqn
     (
         fvm::ddt(psi, p_rgh) + fvc::ddt(psi, rho)*gh
       + fvc::div(phi)
-      - fvm::laplacian(rhorUAf, p_rgh)
+      - fvm::laplacian(rhorAUf, p_rgh)
     );
 
     p_rghEqn.solve
@@ -52,7 +52,7 @@ p = p_rgh + rho*gh;
 #include "rhoEqn.H"
 #include "compressibleContinuityErrs.H"
 
-U += rUA*fvc::reconstruct((phi - phiU)/rhorUAf);
+U += rAU*fvc::reconstruct((phi - phiU)/rhorAUf);
 U.correctBoundaryConditions();
 
 DpDt = fvc::DDt(surfaceScalarField("phiU", phi/fvc::interpolate(rho)), p);
diff --git a/applications/solvers/combustion/reactingFoam/pEqn.H b/applications/solvers/combustion/reactingFoam/pEqn.H
index 9443f909a3..3d414f4af6 100644
--- a/applications/solvers/combustion/reactingFoam/pEqn.H
+++ b/applications/solvers/combustion/reactingFoam/pEqn.H
@@ -1,7 +1,7 @@
 rho = thermo.rho();
 
-volScalarField rUA = 1.0/UEqn.A();
-U = rUA*UEqn.H();
+volScalarField rAU = 1.0/UEqn.A();
+U = rAU*UEqn.H();
 
 if (transonic)
 {
@@ -11,7 +11,7 @@ if (transonic)
         fvc::interpolate(psi)
        *(
             (fvc::interpolate(U) & mesh.Sf())
-          + fvc::ddtPhiCorr(rUA, rho, U, phi)
+          + fvc::ddtPhiCorr(rAU, rho, U, phi)
         )
     );
 
@@ -21,7 +21,7 @@ if (transonic)
         (
             fvm::ddt(psi, p)
           + fvm::div(phid, p)
-          - fvm::laplacian(rho*rUA, p)
+          - fvm::laplacian(rho*rAU, p)
         );
 
         pEqn.solve();
@@ -38,7 +38,7 @@ else
         fvc::interpolate(rho)
        *(
             (fvc::interpolate(U) & mesh.Sf())
-          + fvc::ddtPhiCorr(rUA, rho, U, phi)
+          + fvc::ddtPhiCorr(rAU, rho, U, phi)
         );
 
     for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
@@ -47,7 +47,7 @@ else
         (
             fvm::ddt(psi, p)
           + fvc::div(phi)
-          - fvm::laplacian(rho*rUA, p)
+          - fvm::laplacian(rho*rAU, p)
         );
 
         pEqn.solve();
@@ -62,7 +62,7 @@ else
 #include "rhoEqn.H"
 #include "compressibleContinuityErrs.H"
 
-U -= rUA*fvc::grad(p);
+U -= rAU*fvc::grad(p);
 U.correctBoundaryConditions();
 
 DpDt = fvc::DDt(surfaceScalarField("phiU", phi/fvc::interpolate(rho)), p);
diff --git a/applications/solvers/combustion/rhoReactingFoam/pEqn.H b/applications/solvers/combustion/rhoReactingFoam/pEqn.H
index c76ab5110b..45ae3c1e5f 100644
--- a/applications/solvers/combustion/rhoReactingFoam/pEqn.H
+++ b/applications/solvers/combustion/rhoReactingFoam/pEqn.H
@@ -5,14 +5,14 @@
     // pressure solution - done in 2 parts. Part 1:
     thermo.rho() -= psi*p;
 
-    volScalarField rUA = 1.0/UEqn.A();
-    U = rUA*UEqn.H();
+    volScalarField rAU = 1.0/UEqn.A();
+    U = rAU*UEqn.H();
 
     if (transonic)
     {
         surfaceScalarField phiv =
             (fvc::interpolate(U) & mesh.Sf())
-          + fvc::ddtPhiCorr(rUA, rho, U, phi);
+          + fvc::ddtPhiCorr(rAU, rho, U, phi);
 
         phi = fvc::interpolate(rho)*phiv;
 
@@ -28,7 +28,7 @@
             (
                 fvc::ddt(rho) + fvc::div(phi)
               + correction(fvm::ddt(psi, p) + fvm::div(phid, p))
-              - fvm::laplacian(rho*rUA, p)
+              - fvm::laplacian(rho*rAU, p)
             );
 
             pEqn.solve
@@ -58,7 +58,7 @@
             fvc::interpolate(rho)
            *(
                 (fvc::interpolate(U) & mesh.Sf())
-              + fvc::ddtPhiCorr(rUA, rho, U, phi)
+              + fvc::ddtPhiCorr(rAU, rho, U, phi)
             );
 
         for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
@@ -67,7 +67,7 @@
             (
                 fvc::ddt(rho) + psi*correction(fvm::ddt(p))
               + fvc::div(phi)
-              - fvm::laplacian(rho*rUA, p)
+              - fvm::laplacian(rho*rAU, p)
             );
 
             pEqn.solve
@@ -98,7 +98,7 @@
     #include "rhoEqn.H"
     #include "compressibleContinuityErrs.H"
 
-    U -= rUA*fvc::grad(p);
+    U -= rAU*fvc::grad(p);
     U.correctBoundaryConditions();
 
     DpDt = fvc::DDt(surfaceScalarField("phiU", phi/fvc::interpolate(rho)), p);
diff --git a/applications/solvers/compressible/rhoPimpleFoam/UEqn.H b/applications/solvers/compressible/rhoPimpleFoam/UEqn.H
index 954907d203..45aa8f7790 100644
--- a/applications/solvers/compressible/rhoPimpleFoam/UEqn.H
+++ b/applications/solvers/compressible/rhoPimpleFoam/UEqn.H
@@ -9,7 +9,7 @@ tmp<fvVectorMatrix> UEqn
 
 UEqn().relax();
 
-volScalarField rUA = 1.0/UEqn().A();
+volScalarField rAU = 1.0/UEqn().A();
 
 if (momentumPredictor)
 {
@@ -17,6 +17,6 @@ if (momentumPredictor)
 }
 else
 {
-    U = rUA*(UEqn().H() - fvc::grad(p));
+    U = rAU*(UEqn().H() - fvc::grad(p));
     U.correctBoundaryConditions();
 }
diff --git a/applications/solvers/compressible/rhoPimpleFoam/pEqn.H b/applications/solvers/compressible/rhoPimpleFoam/pEqn.H
index c0dba194e7..84e6451890 100644
--- a/applications/solvers/compressible/rhoPimpleFoam/pEqn.H
+++ b/applications/solvers/compressible/rhoPimpleFoam/pEqn.H
@@ -1,6 +1,6 @@
 rho = thermo.rho();
 
-U = rUA*UEqn().H();
+U = rAU*UEqn().H();
 
 if (nCorr <= 1)
 {
@@ -15,7 +15,7 @@ if (transonic)
         fvc::interpolate(psi)
        *(
             (fvc::interpolate(U) & mesh.Sf())
-          + fvc::ddtPhiCorr(rUA, rho, U, phi)
+          + fvc::ddtPhiCorr(rAU, rho, U, phi)
         )
     );
 
@@ -25,7 +25,7 @@ if (transonic)
         (
             fvm::ddt(psi, p)
           + fvm::div(phid, p)
-          - fvm::laplacian(rho*rUA, p)
+          - fvm::laplacian(rho*rAU, p)
         );
 
         pEqn.solve
@@ -55,7 +55,7 @@ else
         fvc::interpolate(rho)*
         (
             (fvc::interpolate(U) & mesh.Sf())
-          + fvc::ddtPhiCorr(rUA, rho, U, phi)
+          + fvc::ddtPhiCorr(rAU, rho, U, phi)
         );
 
     for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
@@ -65,7 +65,7 @@ else
         (
             fvm::ddt(psi, p)
           + fvc::div(phi)
-          - fvm::laplacian(rho*rUA, p)
+          - fvm::laplacian(rho*rAU, p)
         );
 
         pEqn.solve
@@ -101,7 +101,7 @@ rho = thermo.rho();
 Info<< "rho max/min : " << max(rho).value()
     << " " << min(rho).value() << endl;
 
-U -= rUA*fvc::grad(p);
+U -= rAU*fvc::grad(p);
 U.correctBoundaryConditions();
 
 DpDt = fvc::DDt(surfaceScalarField("phiU", phi/fvc::interpolate(rho)), p);
diff --git a/applications/solvers/compressible/rhoPorousMRFPimpleFoam/UEqn.H b/applications/solvers/compressible/rhoPorousMRFPimpleFoam/UEqn.H
index a331d7dbb8..0b0364b706 100644
--- a/applications/solvers/compressible/rhoPorousMRFPimpleFoam/UEqn.H
+++ b/applications/solvers/compressible/rhoPorousMRFPimpleFoam/UEqn.H
@@ -12,7 +12,7 @@ UEqn().relax();
 mrfZones.addCoriolis(rho, UEqn());
 pZones.addResistance(UEqn());
 
-volScalarField rUA = 1.0/UEqn().A();
+volScalarField rAU = 1.0/UEqn().A();
 
 if (momentumPredictor)
 {
@@ -20,6 +20,6 @@ if (momentumPredictor)
 }
 else
 {
-    U = rUA*(UEqn().H() - fvc::grad(p));
+    U = rAU*(UEqn().H() - fvc::grad(p));
     U.correctBoundaryConditions();
 }
diff --git a/applications/solvers/compressible/rhoPorousMRFPimpleFoam/pEqn.H b/applications/solvers/compressible/rhoPorousMRFPimpleFoam/pEqn.H
index 6be6584202..590a3c1950 100644
--- a/applications/solvers/compressible/rhoPorousMRFPimpleFoam/pEqn.H
+++ b/applications/solvers/compressible/rhoPorousMRFPimpleFoam/pEqn.H
@@ -1,7 +1,7 @@
 rho = thermo.rho();
 
-volScalarField rUA = 1.0/UEqn().A();
-U = rUA*UEqn().H();
+volScalarField rAU = 1.0/UEqn().A();
+U = rAU*UEqn().H();
 
 if (nCorr <= 1)
 {
@@ -16,7 +16,7 @@ if (transonic)
         fvc::interpolate(psi)
        *(
             (fvc::interpolate(U) & mesh.Sf())
-          + fvc::ddtPhiCorr(rUA, rho, U, phi)
+          + fvc::ddtPhiCorr(rAU, rho, U, phi)
         )
     );
     mrfZones.relativeFlux(fvc::interpolate(psi), phid);
@@ -27,7 +27,7 @@ if (transonic)
         (
             fvm::ddt(psi, p)
           + fvm::div(phid, p)
-          - fvm::laplacian(rho*rUA, p)
+          - fvm::laplacian(rho*rAU, p)
         );
 
         pEqn.solve
@@ -57,7 +57,7 @@ else
         fvc::interpolate(rho)*
         (
             (fvc::interpolate(U) & mesh.Sf())
-        //+ fvc::ddtPhiCorr(rUA, rho, U, phi)
+        //+ fvc::ddtPhiCorr(rAU, rho, U, phi)
         );
     mrfZones.relativeFlux(fvc::interpolate(rho), phi);
 
@@ -68,7 +68,7 @@ else
         (
             fvm::ddt(psi, p)
           + fvc::div(phi)
-          - fvm::laplacian(rho*rUA, p)
+          - fvm::laplacian(rho*rAU, p)
         );
 
         pEqn.solve
@@ -109,7 +109,7 @@ else
         << " " << min(rho).value() << endl;
 }
 
-U -= rUA*fvc::grad(p);
+U -= rAU*fvc::grad(p);
 U.correctBoundaryConditions();
 
 DpDt = fvc::DDt(surfaceScalarField("phiU", phi/fvc::interpolate(rho)), p);
diff --git a/applications/solvers/compressible/rhoSimpleFoam/pEqn.H b/applications/solvers/compressible/rhoSimpleFoam/pEqn.H
index 96d8280559..a082cac880 100644
--- a/applications/solvers/compressible/rhoSimpleFoam/pEqn.H
+++ b/applications/solvers/compressible/rhoSimpleFoam/pEqn.H
@@ -3,8 +3,8 @@ rho = max(rho, rhoMin);
 rho = min(rho, rhoMax);
 rho.relax();
 
-volScalarField rUA = 1.0/UEqn().A();
-U = rUA*UEqn().H();
+volScalarField rAU = 1.0/UEqn().A();
+U = rAU*UEqn().H();
 UEqn.clear();
 
 bool closedVolume = false;
@@ -22,7 +22,7 @@ if (transonic)
         fvScalarMatrix pEqn
         (
             fvm::div(phid, p)
-          - fvm::laplacian(rho*rUA, p)
+          - fvm::laplacian(rho*rAU, p)
         );
 
         // Relax the pressure equation to ensure diagonal-dominance
@@ -47,7 +47,7 @@ else
     {
         fvScalarMatrix pEqn
         (
-            fvm::laplacian(rho*rUA, p) == fvc::div(phi)
+            fvm::laplacian(rho*rAU, p) == fvc::div(phi)
         );
 
         pEqn.setReference(pRefCell, pRefValue);
@@ -67,7 +67,7 @@ else
 // Explicitly relax pressure for momentum corrector
 p.relax();
 
-U -= rUA*fvc::grad(p);
+U -= rAU*fvc::grad(p);
 U.correctBoundaryConditions();
 
 // For closed-volume cases adjust the pressure and density levels
diff --git a/applications/solvers/compressible/sonicFoam/pEqn.H b/applications/solvers/compressible/sonicFoam/pEqn.H
index 96a500d4c2..8d6040b530 100644
--- a/applications/solvers/compressible/sonicFoam/pEqn.H
+++ b/applications/solvers/compressible/sonicFoam/pEqn.H
@@ -1,7 +1,7 @@
 rho = thermo.rho();
 
-volScalarField rUA = 1.0/UEqn.A();
-U = rUA*UEqn.H();
+volScalarField rAU = 1.0/UEqn.A();
+U = rAU*UEqn.H();
 
 surfaceScalarField phid
 (
@@ -9,7 +9,7 @@ surfaceScalarField phid
     fvc::interpolate(psi)
    *(
         (fvc::interpolate(U) & mesh.Sf())
-      + fvc::ddtPhiCorr(rUA, rho, U, phi)
+      + fvc::ddtPhiCorr(rAU, rho, U, phi)
     )
 );
 
@@ -19,7 +19,7 @@ for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
     (
         fvm::ddt(psi, p)
       + fvm::div(phid, p)
-      - fvm::laplacian(rho*rUA, p)
+      - fvm::laplacian(rho*rAU, p)
     );
 
     pEqn.solve();
@@ -33,5 +33,5 @@ for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
 #include "rhoEqn.H"
 #include "compressibleContinuityErrs.H"
 
-U -= rUA*fvc::grad(p);
+U -= rAU*fvc::grad(p);
 U.correctBoundaryConditions();
diff --git a/applications/solvers/compressible/sonicLiquidFoam/sonicLiquidFoam.C b/applications/solvers/compressible/sonicLiquidFoam/sonicLiquidFoam.C
index c785343af7..c92f3d82f9 100644
--- a/applications/solvers/compressible/sonicLiquidFoam/sonicLiquidFoam.C
+++ b/applications/solvers/compressible/sonicLiquidFoam/sonicLiquidFoam.C
@@ -71,8 +71,8 @@ int main(int argc, char *argv[])
 
         for (int corr=0; corr<nCorr; corr++)
         {
-            volScalarField rUA = 1.0/UEqn.A();
-            U = rUA*UEqn.H();
+            volScalarField rAU = 1.0/UEqn.A();
+            U = rAU*UEqn.H();
 
             surfaceScalarField phid
             (
@@ -80,7 +80,7 @@ int main(int argc, char *argv[])
                 psi
                *(
                     (fvc::interpolate(U) & mesh.Sf())
-                  + fvc::ddtPhiCorr(rUA, rho, U, phi)
+                  + fvc::ddtPhiCorr(rAU, rho, U, phi)
                 )
             );
 
@@ -91,7 +91,7 @@ int main(int argc, char *argv[])
                 fvm::ddt(psi, p)
               + fvc::div(phi)
               + fvm::div(phid, p)
-              - fvm::laplacian(rho*rUA, p)
+              - fvm::laplacian(rho*rAU, p)
             );
 
             pEqn.solve();
@@ -100,7 +100,7 @@ int main(int argc, char *argv[])
 
             #include "compressibleContinuityErrs.H"
 
-            U -= rUA*fvc::grad(p);
+            U -= rAU*fvc::grad(p);
             U.correctBoundaryConditions();
         }
 
diff --git a/applications/solvers/electromagnetics/mhdFoam/mhdFoam.C b/applications/solvers/electromagnetics/mhdFoam/mhdFoam.C
index 98512923c7..7ebefc77f9 100644
--- a/applications/solvers/electromagnetics/mhdFoam/mhdFoam.C
+++ b/applications/solvers/electromagnetics/mhdFoam/mhdFoam.C
@@ -92,18 +92,18 @@ int main(int argc, char *argv[])
 
             for (int corr=0; corr<nCorr; corr++)
             {
-                volScalarField rUA = 1.0/UEqn.A();
+                volScalarField rAU = 1.0/UEqn.A();
 
-                U = rUA*UEqn.H();
+                U = rAU*UEqn.H();
 
                 phi = (fvc::interpolate(U) & mesh.Sf())
-                    + fvc::ddtPhiCorr(rUA, U, phi);
+                    + fvc::ddtPhiCorr(rAU, U, phi);
 
                 for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
                 {
                     fvScalarMatrix pEqn
                     (
-                        fvm::laplacian(rUA, p) == fvc::div(phi)
+                        fvm::laplacian(rAU, p) == fvc::div(phi)
                     );
 
                     pEqn.setReference(pRefCell, pRefValue);
@@ -117,7 +117,7 @@ int main(int argc, char *argv[])
 
                 #include "continuityErrs.H"
 
-                U -= rUA*fvc::grad(p);
+                U -= rAU*fvc::grad(p);
                 U.correctBoundaryConditions();
             }
         }
diff --git a/applications/solvers/heatTransfer/buoyantBoussinesqPimpleFoam/pEqn.H b/applications/solvers/heatTransfer/buoyantBoussinesqPimpleFoam/pEqn.H
index 60828e18dc..1e972269ff 100644
--- a/applications/solvers/heatTransfer/buoyantBoussinesqPimpleFoam/pEqn.H
+++ b/applications/solvers/heatTransfer/buoyantBoussinesqPimpleFoam/pEqn.H
@@ -1,20 +1,20 @@
 {
-    volScalarField rUA("rUA", 1.0/UEqn.A());
-    surfaceScalarField rUAf("(1|A(U))", fvc::interpolate(rUA));
+    volScalarField rAU("rAU", 1.0/UEqn.A());
+    surfaceScalarField rAUf("(1|A(U))", fvc::interpolate(rAU));
 
-    U = rUA*UEqn.H();
+    U = rAU*UEqn.H();
 
     phi = (fvc::interpolate(U) & mesh.Sf())
-        + fvc::ddtPhiCorr(rUA, U, phi);
+        + fvc::ddtPhiCorr(rAU, U, phi);
 
-    surfaceScalarField buoyancyPhi = rUAf*ghf*fvc::snGrad(rhok)*mesh.magSf();
+    surfaceScalarField buoyancyPhi = rAUf*ghf*fvc::snGrad(rhok)*mesh.magSf();
     phi -= buoyancyPhi;
 
     for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
     {
         fvScalarMatrix p_rghEqn
         (
-            fvm::laplacian(rUAf, p_rgh) == fvc::div(phi)
+            fvm::laplacian(rAUf, p_rgh) == fvc::div(phi)
         );
 
         p_rghEqn.setReference(pRefCell, getRefCellValue(p_rgh, pRefCell));
@@ -44,7 +44,7 @@
 
             // Correct the momentum source with the pressure gradient flux
             // calculated from the relaxed pressure
-            U -= rUA*fvc::reconstruct((buoyancyPhi + p_rghEqn.flux())/rUAf);
+            U -= rAU*fvc::reconstruct((buoyancyPhi + p_rghEqn.flux())/rAUf);
             U.correctBoundaryConditions();
         }
     }
diff --git a/applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/pEqn.H b/applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/pEqn.H
index 50149ed360..daa7867fd5 100644
--- a/applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/pEqn.H
+++ b/applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/pEqn.H
@@ -1,21 +1,21 @@
 {
-    volScalarField rUA("rUA", 1.0/UEqn().A());
-    surfaceScalarField rUAf("(1|A(U))", fvc::interpolate(rUA));
+    volScalarField rAU("rAU", 1.0/UEqn().A());
+    surfaceScalarField rAUf("(1|A(U))", fvc::interpolate(rAU));
 
-    U = rUA*UEqn().H();
+    U = rAU*UEqn().H();
     UEqn.clear();
 
     phi = fvc::interpolate(U) & mesh.Sf();
     adjustPhi(phi, U, p_rgh);
 
-    surfaceScalarField buoyancyPhi = rUAf*ghf*fvc::snGrad(rhok)*mesh.magSf();
+    surfaceScalarField buoyancyPhi = rAUf*ghf*fvc::snGrad(rhok)*mesh.magSf();
     phi -= buoyancyPhi;
 
     for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
     {
         fvScalarMatrix p_rghEqn
         (
-            fvm::laplacian(rUAf, p_rgh) == fvc::div(phi)
+            fvm::laplacian(rAUf, p_rgh) == fvc::div(phi)
         );
 
         p_rghEqn.setReference(pRefCell, getRefCellValue(p_rgh, pRefCell));
@@ -32,7 +32,7 @@
 
             // Correct the momentum source with the pressure gradient flux
             // calculated from the relaxed pressure
-            U -= rUA*fvc::reconstruct((buoyancyPhi + p_rghEqn.flux())/rUAf);
+            U -= rAU*fvc::reconstruct((buoyancyPhi + p_rghEqn.flux())/rAUf);
             U.correctBoundaryConditions();
         }
     }
diff --git a/applications/solvers/heatTransfer/buoyantPimpleFoam/pEqn.H b/applications/solvers/heatTransfer/buoyantPimpleFoam/pEqn.H
index e625f122a3..05863b264c 100644
--- a/applications/solvers/heatTransfer/buoyantPimpleFoam/pEqn.H
+++ b/applications/solvers/heatTransfer/buoyantPimpleFoam/pEqn.H
@@ -5,18 +5,18 @@
     // pressure solution - done in 2 parts. Part 1:
     thermo.rho() -= psi*p_rgh;
 
-    volScalarField rUA = 1.0/UEqn.A();
-    surfaceScalarField rhorUAf("(rho*(1|A(U)))", fvc::interpolate(rho*rUA));
+    volScalarField rAU = 1.0/UEqn.A();
+    surfaceScalarField rhorAUf("(rho*(1|A(U)))", fvc::interpolate(rho*rAU));
 
-    U = rUA*UEqn.H();
+    U = rAU*UEqn.H();
 
     phi = fvc::interpolate(rho)*
     (
         (fvc::interpolate(U) & mesh.Sf())
-      + fvc::ddtPhiCorr(rUA, rho, U, phi)
+      + fvc::ddtPhiCorr(rAU, rho, U, phi)
     );
 
-    surfaceScalarField buoyancyPhi = -rhorUAf*ghf*fvc::snGrad(rho)*mesh.magSf();
+    surfaceScalarField buoyancyPhi = -rhorAUf*ghf*fvc::snGrad(rho)*mesh.magSf();
     phi += buoyancyPhi;
 
     for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
@@ -25,7 +25,7 @@
         (
             fvc::ddt(rho) + psi*correction(fvm::ddt(p_rgh))
           + fvc::div(phi)
-          - fvm::laplacian(rhorUAf, p_rgh)
+          - fvm::laplacian(rhorAUf, p_rgh)
         );
 
         p_rghEqn.solve
@@ -53,7 +53,7 @@
 
             // Correct the momentum source with the pressure gradient flux
             // calculated from the relaxed pressure
-            U += rUA*fvc::reconstruct((buoyancyPhi + p_rghEqn.flux())/rhorUAf);
+            U += rAU*fvc::reconstruct((buoyancyPhi + p_rghEqn.flux())/rhorAUf);
             U.correctBoundaryConditions();
         }
     }
diff --git a/applications/solvers/heatTransfer/buoyantSimpleFoam/pEqn.H b/applications/solvers/heatTransfer/buoyantSimpleFoam/pEqn.H
index f1a62d4770..197005e0f7 100644
--- a/applications/solvers/heatTransfer/buoyantSimpleFoam/pEqn.H
+++ b/applications/solvers/heatTransfer/buoyantSimpleFoam/pEqn.H
@@ -2,23 +2,23 @@
     rho = thermo.rho();
     rho.relax();
 
-    volScalarField rUA = 1.0/UEqn().A();
-    surfaceScalarField rhorUAf("(rho*(1|A(U)))", fvc::interpolate(rho*rUA));
+    volScalarField rAU = 1.0/UEqn().A();
+    surfaceScalarField rhorAUf("(rho*(1|A(U)))", fvc::interpolate(rho*rAU));
 
-    U = rUA*UEqn().H();
+    U = rAU*UEqn().H();
     UEqn.clear();
 
     phi = fvc::interpolate(rho)*(fvc::interpolate(U) & mesh.Sf());
     bool closedVolume = adjustPhi(phi, U, p_rgh);
 
-    surfaceScalarField buoyancyPhi = rhorUAf*ghf*fvc::snGrad(rho)*mesh.magSf();
+    surfaceScalarField buoyancyPhi = rhorAUf*ghf*fvc::snGrad(rho)*mesh.magSf();
     phi -= buoyancyPhi;
 
     for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
     {
         fvScalarMatrix p_rghEqn
         (
-            fvm::laplacian(rhorUAf, p_rgh) == fvc::div(phi)
+            fvm::laplacian(rhorAUf, p_rgh) == fvc::div(phi)
         );
 
         p_rghEqn.setReference(pRefCell, getRefCellValue(p_rgh, pRefCell));
@@ -34,7 +34,7 @@
 
             // Correct the momentum source with the pressure gradient flux
             // calculated from the relaxed pressure
-            U -= rUA*fvc::reconstruct((buoyancyPhi + p_rghEqn.flux())/rhorUAf);
+            U -= rAU*fvc::reconstruct((buoyancyPhi + p_rghEqn.flux())/rhorAUf);
             U.correctBoundaryConditions();
         }
     }
diff --git a/applications/solvers/heatTransfer/buoyantSimpleFoam_old/pEqn.H b/applications/solvers/heatTransfer/buoyantSimpleFoam_old/pEqn.H
index 43625e053e..34536155ef 100644
--- a/applications/solvers/heatTransfer/buoyantSimpleFoam_old/pEqn.H
+++ b/applications/solvers/heatTransfer/buoyantSimpleFoam_old/pEqn.H
@@ -1,24 +1,24 @@
 {
     rho = thermo.rho();
 
-    volScalarField rUA = 1.0/UEqn().A();
-    surfaceScalarField rhorUAf("(rho*(1|A(U)))", fvc::interpolate(rho*rUA));
+    volScalarField rAU = 1.0/UEqn().A();
+    surfaceScalarField rhorAUf("(rho*(1|A(U)))", fvc::interpolate(rho*rAU));
 
-    U = rUA*UEqn().H();
+    U = rAU*UEqn().H();
     UEqn.clear();
 
     phi = fvc::interpolate(rho)*(fvc::interpolate(U) & mesh.Sf());
     bool closedVolume = adjustPhi(phi, U, p);
 
     surfaceScalarField buoyancyPhi =
-        rhorUAf*fvc::interpolate(rho)*(g & mesh.Sf());
+        rhorAUf*fvc::interpolate(rho)*(g & mesh.Sf());
     phi += buoyancyPhi;
 
     for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
     {
         fvScalarMatrix pEqn
         (
-            fvm::laplacian(rhorUAf, p) == fvc::div(phi)
+            fvm::laplacian(rhorAUf, p) == fvc::div(phi)
         );
 
         pEqn.setReference(pRefCell, pRefValue);
@@ -42,8 +42,8 @@
 
             // Correct the momentum source with the pressure gradient flux
             // calculated from the relaxed pressure
-            U += rUA*(rho*g - fvc::grad(p));
-            //U += rUA*fvc::reconstruct((buoyancyPhi - pEqn.flux())/rhorUAf);
+            U += rAU*(rho*g - fvc::grad(p));
+            //U += rAU*fvc::reconstruct((buoyancyPhi - pEqn.flux())/rhorAUf);
             U.correctBoundaryConditions();
         }
     }
diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/pEqn.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/pEqn.H
index 317d377f1b..b0cbbf80bf 100644
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/pEqn.H
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/pEqn.H
@@ -5,16 +5,16 @@
     rho = min(rho, rhoMax[i]);
     rho.relax();
 
-    volScalarField rUA = 1.0/UEqn().A();
-    surfaceScalarField rhorUAf("(rho*(1|A(U)))", fvc::interpolate(rho*rUA));
+    volScalarField rAU = 1.0/UEqn().A();
+    surfaceScalarField rhorAUf("(rho*(1|A(U)))", fvc::interpolate(rho*rAU));
 
-    U = rUA*UEqn().H();
+    U = rAU*UEqn().H();
     UEqn.clear();
 
     phi = fvc::interpolate(rho)*(fvc::interpolate(U) & mesh.Sf());
     bool closedVolume = adjustPhi(phi, U, p_rgh);
 
-    surfaceScalarField buoyancyPhi = rhorUAf*ghf*fvc::snGrad(rho)*mesh.magSf();
+    surfaceScalarField buoyancyPhi = rhorAUf*ghf*fvc::snGrad(rho)*mesh.magSf();
     phi -= buoyancyPhi;
 
     // Solve pressure
@@ -22,7 +22,7 @@
     {
         fvScalarMatrix p_rghEqn
         (
-            fvm::laplacian(rhorUAf, p_rgh) == fvc::div(phi)
+            fvm::laplacian(rhorAUf, p_rgh) == fvc::div(phi)
         );
 
         p_rghEqn.setReference(pRefCell, getRefCellValue(p_rgh, pRefCell));
@@ -39,7 +39,7 @@
 
             // Correct the momentum source with the pressure gradient flux
             // calculated from the relaxed pressure
-            U -= rUA*fvc::reconstruct((buoyancyPhi + p_rghEqn.flux())/rhorUAf);
+            U -= rAU*fvc::reconstruct((buoyancyPhi + p_rghEqn.flux())/rhorAUf);
             U.correctBoundaryConditions();
         }
     }
diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/pEqn.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/pEqn.H
index 9bf3c4acc5..686a8edcee 100644
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/pEqn.H
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/pEqn.H
@@ -3,21 +3,21 @@
 
     rho = thermo.rho();
 
-    volScalarField rUA = 1.0/UEqn().A();
-    surfaceScalarField rhorUAf("(rho*(1|A(U)))", fvc::interpolate(rho*rUA));
+    volScalarField rAU = 1.0/UEqn().A();
+    surfaceScalarField rhorAUf("(rho*(1|A(U)))", fvc::interpolate(rho*rAU));
 
-    U = rUA*UEqn().H();
+    U = rAU*UEqn().H();
 
     surfaceScalarField phiU
     (
         fvc::interpolate(rho)
        *(
             (fvc::interpolate(U) & mesh.Sf())
-          + fvc::ddtPhiCorr(rUA, rho, U, phi)
+          + fvc::ddtPhiCorr(rAU, rho, U, phi)
         )
     );
 
-    phi = phiU - rhorUAf*ghf*fvc::snGrad(rho)*mesh.magSf();
+    phi = phiU - rhorAUf*ghf*fvc::snGrad(rho)*mesh.magSf();
 
     for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
     {
@@ -25,7 +25,7 @@
         (
             fvm::ddt(psi, p_rgh) + fvc::ddt(psi, rho)*gh
           + fvc::div(phi)
-          - fvm::laplacian(rhorUAf, p_rgh)
+          - fvm::laplacian(rhorAUf, p_rgh)
         );
 
         p_rghEqn.solve
@@ -50,7 +50,7 @@
     }
 
     // Correct velocity field
-    U += rUA*fvc::reconstruct((phi - phiU)/rhorUAf);
+    U += rAU*fvc::reconstruct((phi - phiU)/rhorAUf);
     U.correctBoundaryConditions();
 
     p = p_rgh + rho*gh;
diff --git a/applications/solvers/incompressible/channelFoam/channelFoam.C b/applications/solvers/incompressible/channelFoam/channelFoam.C
index db1e0d4457..8918a4d8eb 100644
--- a/applications/solvers/incompressible/channelFoam/channelFoam.C
+++ b/applications/solvers/incompressible/channelFoam/channelFoam.C
@@ -77,13 +77,13 @@ int main(int argc, char *argv[])
 
         // --- PISO loop
 
-        volScalarField rUA = 1.0/UEqn.A();
+        volScalarField rAU = 1.0/UEqn.A();
 
         for (int corr=0; corr<nCorr; corr++)
         {
-            U = rUA*UEqn.H();
+            U = rAU*UEqn.H();
             phi = (fvc::interpolate(U) & mesh.Sf())
-                + fvc::ddtPhiCorr(rUA, U, phi);
+                + fvc::ddtPhiCorr(rAU, U, phi);
 
             adjustPhi(phi, U, p);
 
@@ -91,7 +91,7 @@ int main(int argc, char *argv[])
             {
                 fvScalarMatrix pEqn
                 (
-                    fvm::laplacian(rUA, p) == fvc::div(phi)
+                    fvm::laplacian(rAU, p) == fvc::div(phi)
                 );
 
                 pEqn.setReference(pRefCell, pRefValue);
@@ -113,7 +113,7 @@ int main(int argc, char *argv[])
 
             #include "continuityErrs.H"
 
-            U -= rUA*fvc::grad(p);
+            U -= rAU*fvc::grad(p);
             U.correctBoundaryConditions();
         }
 
@@ -127,9 +127,9 @@ int main(int argc, char *argv[])
         // Calculate the pressure gradient increment needed to
         // adjust the average flow-rate to the correct value
         dimensionedScalar gragPplus =
-            (magUbar - magUbarStar)/rUA.weightedAverage(mesh.V());
+            (magUbar - magUbarStar)/rAU.weightedAverage(mesh.V());
 
-        U += flowDirection*rUA*gragPplus;
+        U += flowDirection*rAU*gragPplus;
 
         gradP += gragPplus;
 
diff --git a/applications/solvers/incompressible/icoFoam/icoFoam.C b/applications/solvers/incompressible/icoFoam/icoFoam.C
index 8a4b61761d..98a4c4324b 100644
--- a/applications/solvers/incompressible/icoFoam/icoFoam.C
+++ b/applications/solvers/incompressible/icoFoam/icoFoam.C
@@ -66,11 +66,11 @@ int main(int argc, char *argv[])
 
         for (int corr=0; corr<nCorr; corr++)
         {
-            volScalarField rUA = 1.0/UEqn.A();
+            volScalarField rAU = 1.0/UEqn.A();
 
-            U = rUA*UEqn.H();
+            U = rAU*UEqn.H();
             phi = (fvc::interpolate(U) & mesh.Sf())
-                + fvc::ddtPhiCorr(rUA, U, phi);
+                + fvc::ddtPhiCorr(rAU, U, phi);
 
             adjustPhi(phi, U, p);
 
@@ -78,7 +78,7 @@ int main(int argc, char *argv[])
             {
                 fvScalarMatrix pEqn
                 (
-                    fvm::laplacian(rUA, p) == fvc::div(phi)
+                    fvm::laplacian(rAU, p) == fvc::div(phi)
                 );
 
                 pEqn.setReference(pRefCell, pRefValue);
@@ -92,7 +92,7 @@ int main(int argc, char *argv[])
 
             #include "continuityErrs.H"
 
-            U -= rUA*fvc::grad(p);
+            U -= rAU*fvc::grad(p);
             U.correctBoundaryConditions();
         }
 
diff --git a/applications/solvers/incompressible/nonNewtonianIcoFoam/nonNewtonianIcoFoam.C b/applications/solvers/incompressible/nonNewtonianIcoFoam/nonNewtonianIcoFoam.C
index 022796b064..dc62df435c 100644
--- a/applications/solvers/incompressible/nonNewtonianIcoFoam/nonNewtonianIcoFoam.C
+++ b/applications/solvers/incompressible/nonNewtonianIcoFoam/nonNewtonianIcoFoam.C
@@ -69,11 +69,11 @@ int main(int argc, char *argv[])
 
         for (int corr=0; corr<nCorr; corr++)
         {
-            volScalarField rUA = 1.0/UEqn.A();
+            volScalarField rAU = 1.0/UEqn.A();
 
-            U = rUA*UEqn.H();
+            U = rAU*UEqn.H();
             phi = (fvc::interpolate(U) & mesh.Sf())
-                + fvc::ddtPhiCorr(rUA, U, phi);
+                + fvc::ddtPhiCorr(rAU, U, phi);
 
             adjustPhi(phi, U, p);
 
@@ -81,7 +81,7 @@ int main(int argc, char *argv[])
             {
                 fvScalarMatrix pEqn
                 (
-                    fvm::laplacian(rUA, p) == fvc::div(phi)
+                    fvm::laplacian(rAU, p) == fvc::div(phi)
                 );
 
                 pEqn.setReference(pRefCell, pRefValue);
@@ -95,7 +95,7 @@ int main(int argc, char *argv[])
 
             #include "continuityErrs.H"
 
-            U -= rUA*fvc::grad(p);
+            U -= rAU*fvc::grad(p);
             U.correctBoundaryConditions();
         }
 
diff --git a/applications/solvers/incompressible/pimpleFoam/pimpleDyMFoam/UEqn.H b/applications/solvers/incompressible/pimpleFoam/pimpleDyMFoam/UEqn.H
new file mode 100644
index 0000000000..652773fea7
--- /dev/null
+++ b/applications/solvers/incompressible/pimpleFoam/pimpleDyMFoam/UEqn.H
@@ -0,0 +1,22 @@
+// Solve the Momentum equation
+
+tmp<fvVectorMatrix> UEqn
+(
+    fvm::ddt(U)
+  + fvm::div(phi, U)
+  + turbulence->divDevReff(U)
+);
+
+UEqn().relax();
+
+rAU = 1.0/UEqn().A();
+
+if (momentumPredictor)
+{
+    solve(UEqn() == -fvc::grad(p));
+}
+else
+{
+    U = rAU*(UEqn().H() - fvc::grad(p));
+    U.correctBoundaryConditions();
+}
diff --git a/applications/solvers/incompressible/pisoFoam/pisoFoam.C b/applications/solvers/incompressible/pisoFoam/pisoFoam.C
index 5f7e6dd8e0..47414bf173 100644
--- a/applications/solvers/incompressible/pisoFoam/pisoFoam.C
+++ b/applications/solvers/incompressible/pisoFoam/pisoFoam.C
@@ -79,11 +79,11 @@ int main(int argc, char *argv[])
 
             for (int corr=0; corr<nCorr; corr++)
             {
-                volScalarField rUA = 1.0/UEqn.A();
+                volScalarField rAU = 1.0/UEqn.A();
 
-                U = rUA*UEqn.H();
+                U = rAU*UEqn.H();
                 phi = (fvc::interpolate(U) & mesh.Sf())
-                    + fvc::ddtPhiCorr(rUA, U, phi);
+                    + fvc::ddtPhiCorr(rAU, U, phi);
 
                 adjustPhi(phi, U, p);
 
@@ -94,7 +94,7 @@ int main(int argc, char *argv[])
 
                     fvScalarMatrix pEqn
                     (
-                        fvm::laplacian(rUA, p) == fvc::div(phi)
+                        fvm::laplacian(rAU, p) == fvc::div(phi)
                     );
 
                     pEqn.setReference(pRefCell, pRefValue);
@@ -120,7 +120,7 @@ int main(int argc, char *argv[])
 
                 #include "continuityErrs.H"
 
-                U -= rUA*fvc::grad(p);
+                U -= rAU*fvc::grad(p);
                 U.correctBoundaryConditions();
             }
         }
diff --git a/applications/solvers/incompressible/shallowWaterFoam/shallowWaterFoam.C b/applications/solvers/incompressible/shallowWaterFoam/shallowWaterFoam.C
index a7b5602c57..274c8d1a63 100644
--- a/applications/solvers/incompressible/shallowWaterFoam/shallowWaterFoam.C
+++ b/applications/solvers/incompressible/shallowWaterFoam/shallowWaterFoam.C
@@ -89,22 +89,22 @@ int main(int argc, char *argv[])
             // --- PISO loop
             for (int corr=0; corr<nCorr; corr++)
             {
-                volScalarField rUA = 1.0/hUEqn.A();
-                surfaceScalarField ghrUAf = magg*fvc::interpolate(h*rUA);
+                volScalarField rAU = 1.0/hUEqn.A();
+                surfaceScalarField ghrAUf = magg*fvc::interpolate(h*rAU);
 
-                surfaceScalarField phih0 = ghrUAf*mesh.magSf()*fvc::snGrad(h0);
+                surfaceScalarField phih0 = ghrAUf*mesh.magSf()*fvc::snGrad(h0);
 
                 if (rotating)
                 {
-                    hU = rUA*(hUEqn.H() - (F ^ hU));
+                    hU = rAU*(hUEqn.H() - (F ^ hU));
                 }
                 else
                 {
-                    hU = rUA*hUEqn.H();
+                    hU = rAU*hUEqn.H();
                 }
 
                 phi = (fvc::interpolate(hU) & mesh.Sf())
-                    + fvc::ddtPhiCorr(rUA, h, hU, phi)
+                    + fvc::ddtPhiCorr(rAU, h, hU, phi)
                     - phih0;
 
                 for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
@@ -113,7 +113,7 @@ int main(int argc, char *argv[])
                     (
                         fvm::ddt(h)
                       + fvc::div(phi)
-                      - fvm::laplacian(ghrUAf, h)
+                      - fvm::laplacian(ghrAUf, h)
                     );
 
                     if (ucorr < nOuterCorr-1 || corr < nCorr-1)
@@ -131,7 +131,7 @@ int main(int argc, char *argv[])
                     }
                 }
 
-                hU -= rUA*h*magg*fvc::grad(h + h0);
+                hU -= rAU*h*magg*fvc::grad(h + h0);
 
                 // Constrain the momentum to be in the geometry if 3D geometry
                 if (mesh.nGeometricD() == 3)
diff --git a/applications/solvers/lagrangian/coalChemistryFoam/pEqn.H b/applications/solvers/lagrangian/coalChemistryFoam/pEqn.H
index d3efde0b34..3f007e247d 100644
--- a/applications/solvers/lagrangian/coalChemistryFoam/pEqn.H
+++ b/applications/solvers/lagrangian/coalChemistryFoam/pEqn.H
@@ -1,7 +1,7 @@
 rho = thermo.rho();
 
-volScalarField rUA = 1.0/UEqn.A();
-U = rUA*UEqn.H();
+volScalarField rAU = 1.0/UEqn.A();
+U = rAU*UEqn.H();
 
 if (transonic)
 {
@@ -11,7 +11,7 @@ if (transonic)
         fvc::interpolate(psi)
        *(
             (fvc::interpolate(U) & mesh.Sf())
-          + fvc::ddtPhiCorr(rUA, rho, U, phi)
+          + fvc::ddtPhiCorr(rAU, rho, U, phi)
         )
     );
 
@@ -21,7 +21,7 @@ if (transonic)
         (
             fvm::ddt(psi, p)
           + fvm::div(phid, p)
-          - fvm::laplacian(rho*rUA, p)
+          - fvm::laplacian(rho*rAU, p)
          ==
             coalParcels.Srho()
         );
@@ -53,7 +53,7 @@ else
         fvc::interpolate(rho)*
         (
             (fvc::interpolate(U) & mesh.Sf())
-          + fvc::ddtPhiCorr(rUA, rho, U, phi)
+          + fvc::ddtPhiCorr(rAU, rho, U, phi)
         );
 
     for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
@@ -62,7 +62,7 @@ else
         (
             fvm::ddt(psi, p)
           + fvc::div(phi)
-          - fvm::laplacian(rho*rUA, p)
+          - fvm::laplacian(rho*rAU, p)
          ==
             coalParcels.Srho()
         );
@@ -92,7 +92,7 @@ else
 #include "rhoEqn.H"
 #include "compressibleContinuityErrs.H"
 
-U -= rUA*fvc::grad(p);
+U -= rAU*fvc::grad(p);
 U.correctBoundaryConditions();
 
 DpDt = fvc::DDt(surfaceScalarField("phiU", phi/fvc::interpolate(rho)), p);
diff --git a/applications/solvers/lagrangian/reactingParcelFilmFoam/pEqn.H b/applications/solvers/lagrangian/reactingParcelFilmFoam/pEqn.H
index 93764f87de..f957ab017e 100644
--- a/applications/solvers/lagrangian/reactingParcelFilmFoam/pEqn.H
+++ b/applications/solvers/lagrangian/reactingParcelFilmFoam/pEqn.H
@@ -1,7 +1,7 @@
 rho = thermo.rho();
 
-volScalarField rUA = 1.0/UEqn.A();
-U = rUA*UEqn.H();
+volScalarField rAU = 1.0/UEqn.A();
+U = rAU*UEqn.H();
 
 if (transonic)
 {
@@ -11,7 +11,7 @@ if (transonic)
         fvc::interpolate(psi)
        *(
             (fvc::interpolate(U) & mesh.Sf())
-          + fvc::ddtPhiCorr(rUA, rho, U, phi)
+          + fvc::ddtPhiCorr(rAU, rho, U, phi)
         )
     );
 
@@ -21,7 +21,7 @@ if (transonic)
         (
             fvm::ddt(psi, p)
           + fvm::div(phid, p)
-          - fvm::laplacian(rho*rUA, p)
+          - fvm::laplacian(rho*rAU, p)
          ==
             parcels.Srho()
           + surfaceFilm.Srho()
@@ -41,7 +41,7 @@ else
         fvc::interpolate(rho)
        *(
             (fvc::interpolate(U) & mesh.Sf())
-          + fvc::ddtPhiCorr(rUA, rho, U, phi)
+          + fvc::ddtPhiCorr(rAU, rho, U, phi)
         );
 
     for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
@@ -50,7 +50,7 @@ else
         (
             fvm::ddt(psi, p)
           + fvc::div(phi)
-          - fvm::laplacian(rho*rUA, p)
+          - fvm::laplacian(rho*rAU, p)
          ==
             parcels.Srho()
           + surfaceFilm.Srho()
@@ -68,7 +68,7 @@ else
 #include "rhoEqn.H"
 #include "compressibleContinuityErrs.H"
 
-U -= rUA*fvc::grad(p);
+U -= rAU*fvc::grad(p);
 U.correctBoundaryConditions();
 
 DpDt = fvc::DDt(surfaceScalarField("phiU", phi/fvc::interpolate(rho)), p);
diff --git a/applications/solvers/lagrangian/reactingParcelFoam/pEqn.H b/applications/solvers/lagrangian/reactingParcelFoam/pEqn.H
index d27afcce8a..d5de608bd2 100644
--- a/applications/solvers/lagrangian/reactingParcelFoam/pEqn.H
+++ b/applications/solvers/lagrangian/reactingParcelFoam/pEqn.H
@@ -1,7 +1,7 @@
 rho = thermo.rho();
 
-volScalarField rUA = 1.0/UEqn.A();
-U = rUA*UEqn.H();
+volScalarField rAU = 1.0/UEqn.A();
+U = rAU*UEqn.H();
 
 if (transonic)
 {
@@ -11,7 +11,7 @@ if (transonic)
         fvc::interpolate(psi)
        *(
             (fvc::interpolate(U) & mesh.Sf())
-          + fvc::ddtPhiCorr(rUA, rho, U, phi)
+          + fvc::ddtPhiCorr(rAU, rho, U, phi)
         )
     );
 
@@ -21,7 +21,7 @@ if (transonic)
         (
             fvm::ddt(psi, p)
           + fvm::div(phid, p)
-          - fvm::laplacian(rho*rUA, p)
+          - fvm::laplacian(rho*rAU, p)
          ==
             parcels.Srho()
         );
@@ -40,7 +40,7 @@ else
         fvc::interpolate(rho)
        *(
             (fvc::interpolate(U) & mesh.Sf())
-          + fvc::ddtPhiCorr(rUA, rho, U, phi)
+          + fvc::ddtPhiCorr(rAU, rho, U, phi)
         );
 
     for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
@@ -49,7 +49,7 @@ else
         (
             fvm::ddt(psi, p)
           + fvc::div(phi)
-          - fvm::laplacian(rho*rUA, p)
+          - fvm::laplacian(rho*rAU, p)
          ==
             parcels.Srho()
         );
@@ -66,7 +66,7 @@ else
 #include "rhoEqn.H"
 #include "compressibleContinuityErrs.H"
 
-U -= rUA*fvc::grad(p);
+U -= rAU*fvc::grad(p);
 U.correctBoundaryConditions();
 
 DpDt = fvc::DDt(surfaceScalarField("phiU", phi/fvc::interpolate(rho)), p);
diff --git a/applications/solvers/multiphase/cavitatingFoam/pEqn.H b/applications/solvers/multiphase/cavitatingFoam/pEqn.H
index bcaedb7a1c..e6606b7b62 100644
--- a/applications/solvers/multiphase/cavitatingFoam/pEqn.H
+++ b/applications/solvers/multiphase/cavitatingFoam/pEqn.H
@@ -11,16 +11,16 @@
 
     surfaceScalarField rhof = fvc::interpolate(rho, "rhof");
 
-    volScalarField rUA = 1.0/UEqn.A();
-    surfaceScalarField rUAf("rUAf", rhof*fvc::interpolate(rUA));
-    volVectorField HbyA = rUA*UEqn.H();
+    volScalarField rAU = 1.0/UEqn.A();
+    surfaceScalarField rAUf("rAUf", rhof*fvc::interpolate(rAU));
+    volVectorField HbyA = rAU*UEqn.H();
 
     phiv = (fvc::interpolate(HbyA) & mesh.Sf())
-         + fvc::ddtPhiCorr(rUA, rho, U, phiv);
+         + fvc::ddtPhiCorr(rAU, rho, U, phiv);
 
     p.boundaryField().updateCoeffs();
 
-    surfaceScalarField phiGradp = rUAf*mesh.magSf()*fvc::snGrad(p);
+    surfaceScalarField phiGradp = rAUf*mesh.magSf()*fvc::snGrad(p);
 
     phiv -= phiGradp/rhof;
 
@@ -34,7 +34,7 @@
           - (rhol0 + (psil - psiv)*pSat)*fvc::ddt(gamma) - pSat*fvc::ddt(psi)
           + fvc::div(phiv, rho)
           + fvc::div(phiGradp)
-          - fvm::laplacian(rUAf, p)
+          - fvm::laplacian(rAUf, p)
         );
 
         if (corr == nCorr-1 && nonOrth == nNonOrthCorr)
@@ -79,7 +79,7 @@
 
     // Correct velocity
 
-    U = HbyA - rUA*fvc::grad(p);
+    U = HbyA - rAU*fvc::grad(p);
 
     // Remove the swirl component of velocity for "wedge" cases
     if (piso.found("removeSwirl"))
diff --git a/applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/pEqn.H b/applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/pEqn.H
index 0df0847c07..ace51b5d84 100644
--- a/applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/pEqn.H
+++ b/applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/pEqn.H
@@ -1,6 +1,6 @@
 {
-    volScalarField rUA = 1.0/UEqn.A();
-    surfaceScalarField rUAf = fvc::interpolate(rUA);
+    volScalarField rAU = 1.0/UEqn.A();
+    surfaceScalarField rAUf = fvc::interpolate(rAU);
 
     tmp<fvScalarMatrix> p_rghEqnComp;
 
@@ -24,27 +24,27 @@
     }
 
 
-    U = rUA*UEqn.H();
+    U = rAU*UEqn.H();
 
     surfaceScalarField phiU
     (
         "phiU",
         (fvc::interpolate(U) & mesh.Sf())
-      + fvc::ddtPhiCorr(rUA, rho, U, phi)
+      + fvc::ddtPhiCorr(rAU, rho, U, phi)
     );
 
     phi = phiU +
         (
             fvc::interpolate(interface.sigmaK())*fvc::snGrad(alpha1)
           - ghf*fvc::snGrad(rho)
-        )*rUAf*mesh.magSf();
+        )*rAUf*mesh.magSf();
 
     for(int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
     {
         fvScalarMatrix p_rghEqnIncomp
         (
             fvc::div(phi)
-          - fvm::laplacian(rUAf, p_rgh)
+          - fvm::laplacian(rAUf, p_rgh)
         );
 
         solve
@@ -75,7 +75,7 @@
         }
     }
 
-    U += rUA*fvc::reconstruct((phi - phiU)/rUAf);
+    U += rAU*fvc::reconstruct((phi - phiU)/rAUf);
     U.correctBoundaryConditions();
 
     p = max
diff --git a/applications/solvers/multiphase/compressibleInterFoam/pEqn.H b/applications/solvers/multiphase/compressibleInterFoam/pEqn.H
index eefe6bec08..7054fe1d44 100644
--- a/applications/solvers/multiphase/compressibleInterFoam/pEqn.H
+++ b/applications/solvers/multiphase/compressibleInterFoam/pEqn.H
@@ -1,6 +1,6 @@
 {
-    volScalarField rUA = 1.0/UEqn.A();
-    surfaceScalarField rUAf = fvc::interpolate(rUA);
+    volScalarField rAU = 1.0/UEqn.A();
+    surfaceScalarField rAUf = fvc::interpolate(rAU);
 
     tmp<fvScalarMatrix> p_rghEqnComp;
 
@@ -24,27 +24,27 @@
     }
 
 
-    U = rUA*UEqn.H();
+    U = rAU*UEqn.H();
 
     surfaceScalarField phiU
     (
         "phiU",
         (fvc::interpolate(U) & mesh.Sf())
-      + fvc::ddtPhiCorr(rUA, rho, U, phi)
+      + fvc::ddtPhiCorr(rAU, rho, U, phi)
     );
 
     phi = phiU +
         (
             fvc::interpolate(interface.sigmaK())*fvc::snGrad(alpha1)
           - ghf*fvc::snGrad(rho)
-        )*rUAf*mesh.magSf();
+        )*rAUf*mesh.magSf();
 
     for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
     {
         fvScalarMatrix p_rghEqnIncomp
         (
             fvc::div(phi)
-          - fvm::laplacian(rUAf, p_rgh)
+          - fvm::laplacian(rAUf, p_rgh)
         );
 
         solve
@@ -75,7 +75,7 @@
         }
     }
 
-    U += rUA*fvc::reconstruct((phi - phiU)/rUAf);
+    U += rAU*fvc::reconstruct((phi - phiU)/rAUf);
     U.correctBoundaryConditions();
 
     p = max
diff --git a/applications/solvers/multiphase/interPhaseChangeFoam/pEqn.H b/applications/solvers/multiphase/interPhaseChangeFoam/pEqn.H
index cb42a4474a..4f4e68cb59 100644
--- a/applications/solvers/multiphase/interPhaseChangeFoam/pEqn.H
+++ b/applications/solvers/multiphase/interPhaseChangeFoam/pEqn.H
@@ -1,14 +1,14 @@
 {
-    volScalarField rUA = 1.0/UEqn.A();
-    surfaceScalarField rUAf = fvc::interpolate(rUA);
+    volScalarField rAU = 1.0/UEqn.A();
+    surfaceScalarField rAUf = fvc::interpolate(rAU);
 
-    U = rUA*UEqn.H();
+    U = rAU*UEqn.H();
 
     surfaceScalarField phiU
     (
         "phiU",
         (fvc::interpolate(U) & mesh.Sf())
-      + fvc::ddtPhiCorr(rUA, rho, U, phi)
+      + fvc::ddtPhiCorr(rAU, rho, U, phi)
     );
 
     adjustPhi(phiU, U, p_rgh);
@@ -17,7 +17,7 @@
         (
             fvc::interpolate(interface.sigmaK())*fvc::snGrad(alpha1)
           - ghf*fvc::snGrad(rho)
-        )*rUAf*mesh.magSf();
+        )*rAUf*mesh.magSf();
 
     Pair<tmp<volScalarField> > vDotP = twoPhaseProperties->vDotP();
     const volScalarField& vDotcP = vDotP[0]();
@@ -27,7 +27,7 @@
     {
         fvScalarMatrix p_rghEqn
         (
-            fvc::div(phi) - fvm::laplacian(rUAf, p_rgh)
+            fvc::div(phi) - fvm::laplacian(rAUf, p_rgh)
           - (vDotvP - vDotcP)*(pSat - rho*gh) + fvm::Sp(vDotvP - vDotcP, p_rgh)
         );
 
@@ -52,7 +52,7 @@
         }
     }
 
-    U += rUA*fvc::reconstruct((phi - phiU)/rUAf);
+    U += rAU*fvc::reconstruct((phi - phiU)/rAUf);
     U.correctBoundaryConditions();
 
     #include "continuityErrs.H"
diff --git a/applications/solvers/multiphase/settlingFoam/pEqn.H b/applications/solvers/multiphase/settlingFoam/pEqn.H
index 3d1e6816d7..806909ffe7 100644
--- a/applications/solvers/multiphase/settlingFoam/pEqn.H
+++ b/applications/solvers/multiphase/settlingFoam/pEqn.H
@@ -1,27 +1,27 @@
-volScalarField rUA = 1.0/UEqn.A();
+volScalarField rAU = 1.0/UEqn.A();
 
-surfaceScalarField rUAf
+surfaceScalarField rAUf
 (
     "(rho*(1|A(U)))",
-    fvc::interpolate(rho)*fvc::interpolate(rUA)
+    fvc::interpolate(rho)*fvc::interpolate(rAU)
 );
 
-U = rUA*UEqn.H();
+U = rAU*UEqn.H();
 phi =
     fvc::interpolate(rho)
    *(
        (fvc::interpolate(U) & mesh.Sf())
-     + fvc::ddtPhiCorr(rUA, rho, U, phi)
+     + fvc::ddtPhiCorr(rAU, rho, U, phi)
     );
 
 surfaceScalarField phiU("phiU", phi);
-phi -= ghf*fvc::snGrad(rho)*rUAf*mesh.magSf();
+phi -= ghf*fvc::snGrad(rho)*rAUf*mesh.magSf();
 
 for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
 {
     fvScalarMatrix p_rghEqn
     (
-        fvm::laplacian(rUAf, p_rgh) == fvc::ddt(rho) + fvc::div(phi)
+        fvm::laplacian(rAUf, p_rgh) == fvc::ddt(rho) + fvc::div(phi)
     );
 
     p_rghEqn.setReference(pRefCell, getRefCellValue(p_rgh, pRefCell));
@@ -49,5 +49,5 @@ if (p_rgh.needReference())
 #include "rhoEqn.H"
 #include "compressibleContinuityErrs.H"
 
-U += rUA*fvc::reconstruct((phi - phiU)/rUAf);
+U += rAU*fvc::reconstruct((phi - phiU)/rAUf);
 U.correctBoundaryConditions();
diff --git a/applications/test/LduMatrix/LduMatrixTest3.C b/applications/test/LduMatrix/LduMatrixTest3.C
index 741e07d990..35e8a4ab85 100644
--- a/applications/test/LduMatrix/LduMatrixTest3.C
+++ b/applications/test/LduMatrix/LduMatrixTest3.C
@@ -102,11 +102,11 @@ int main(int argc, char *argv[])
 
         for (int corr=0; corr<nCorr; corr++)
         {
-            volScalarField rUA = 1.0/UEqn.A();
+            volScalarField rAU = 1.0/UEqn.A();
 
-            U = rUA*UEqn.H();
+            U = rAU*UEqn.H();
             phi = (fvc::interpolate(U) & mesh.Sf()) 
-                + fvc::ddtPhiCorr(rUA, U, phi);
+                + fvc::ddtPhiCorr(rAU, U, phi);
 
             adjustPhi(phi, U, p);
 
@@ -114,7 +114,7 @@ int main(int argc, char *argv[])
             {
                 fvScalarMatrix pEqn
                 (
-                    fvm::laplacian(rUA, p) == fvc::div(phi)
+                    fvm::laplacian(rAU, p) == fvc::div(phi)
                 );
 
                 pEqn.setReference(pRefCell, pRefValue);
@@ -128,7 +128,7 @@ int main(int argc, char *argv[])
 
 #           include "continuityErrs.H"
 
-            U -= rUA*fvc::grad(p);
+            U -= rAU*fvc::grad(p);
             U.correctBoundaryConditions();
         }
 
diff --git a/src/finiteVolume/cfdTools/general/fieldSources/pressureGradientExplicitSource/pressureGradientExplicitSource.C b/src/finiteVolume/cfdTools/general/fieldSources/pressureGradientExplicitSource/pressureGradientExplicitSource.C
index c08ded14a2..fd0ed25b2c 100644
--- a/src/finiteVolume/cfdTools/general/fieldSources/pressureGradientExplicitSource/pressureGradientExplicitSource.C
+++ b/src/finiteVolume/cfdTools/general/fieldSources/pressureGradientExplicitSource/pressureGradientExplicitSource.C
@@ -161,13 +161,13 @@ Foam::pressureGradientExplicitSource::Su() const
 
 void Foam::pressureGradientExplicitSource::update()
 {
-    const volScalarField& rUA =
+    const volScalarField& rAU =
         mesh_.lookupObject<volScalarField>("(1|A(" + U_.name() + "))");
 
     // Integrate flow variables over cell set
     scalar volTot = 0.0;
     scalar magUbarAve = 0.0;
-    scalar rUAave = 0.0;
+    scalar rAUave = 0.0;
     forAllConstIter(cellSet, selectedCellSet_, iter)
     {
         label cellI = iter.key();
@@ -176,27 +176,27 @@ void Foam::pressureGradientExplicitSource::update()
         volTot += volCell;
 
         magUbarAve += (flowDir_ & U_[cellI])*volCell;
-        rUAave += rUA[cellI]*volCell;
+        rAUave += rAU[cellI]*volCell;
     }
 
     // Collect across all processors
     reduce(volTot, sumOp<scalar>());
     reduce(magUbarAve, sumOp<scalar>());
-    reduce(rUAave, sumOp<scalar>());
+    reduce(rAUave, sumOp<scalar>());
 
     // Volume averages
     magUbarAve /= volTot;
-    rUAave /= volTot;
+    rAUave /= volTot;
 
     // Calculate the pressure gradient increment needed to adjust the average
     // flow-rate to the desired value
-    scalar gradPplus = (mag(Ubar_) - magUbarAve)/rUAave;
+    scalar gradPplus = (mag(Ubar_) - magUbarAve)/rAUave;
 
     // Apply correction to velocity field
     forAllConstIter(cellSet, selectedCellSet_, iter)
     {
         label cellI = iter.key();
-        U_[cellI] += flowDir_*rUA[cellI]*gradPplus;
+        U_[cellI] += flowDir_*rAU[cellI]*gradPplus;
     }
 
     // Update pressure gradient
diff --git a/tutorials/lagrangian/rhoPisoTwinParcelFoam/rhoPisoTwinParcelFoam/pEqn.H b/tutorials/lagrangian/rhoPisoTwinParcelFoam/rhoPisoTwinParcelFoam/pEqn.H
index 9443f909a3..3d414f4af6 100644
--- a/tutorials/lagrangian/rhoPisoTwinParcelFoam/rhoPisoTwinParcelFoam/pEqn.H
+++ b/tutorials/lagrangian/rhoPisoTwinParcelFoam/rhoPisoTwinParcelFoam/pEqn.H
@@ -1,7 +1,7 @@
 rho = thermo.rho();
 
-volScalarField rUA = 1.0/UEqn.A();
-U = rUA*UEqn.H();
+volScalarField rAU = 1.0/UEqn.A();
+U = rAU*UEqn.H();
 
 if (transonic)
 {
@@ -11,7 +11,7 @@ if (transonic)
         fvc::interpolate(psi)
        *(
             (fvc::interpolate(U) & mesh.Sf())
-          + fvc::ddtPhiCorr(rUA, rho, U, phi)
+          + fvc::ddtPhiCorr(rAU, rho, U, phi)
         )
     );
 
@@ -21,7 +21,7 @@ if (transonic)
         (
             fvm::ddt(psi, p)
           + fvm::div(phid, p)
-          - fvm::laplacian(rho*rUA, p)
+          - fvm::laplacian(rho*rAU, p)
         );
 
         pEqn.solve();
@@ -38,7 +38,7 @@ else
         fvc::interpolate(rho)
        *(
             (fvc::interpolate(U) & mesh.Sf())
-          + fvc::ddtPhiCorr(rUA, rho, U, phi)
+          + fvc::ddtPhiCorr(rAU, rho, U, phi)
         );
 
     for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
@@ -47,7 +47,7 @@ else
         (
             fvm::ddt(psi, p)
           + fvc::div(phi)
-          - fvm::laplacian(rho*rUA, p)
+          - fvm::laplacian(rho*rAU, p)
         );
 
         pEqn.solve();
@@ -62,7 +62,7 @@ else
 #include "rhoEqn.H"
 #include "compressibleContinuityErrs.H"
 
-U -= rUA*fvc::grad(p);
+U -= rAU*fvc::grad(p);
 U.correctBoundaryConditions();
 
 DpDt = fvc::DDt(surfaceScalarField("phiU", phi/fvc::interpolate(rho)), p);
diff --git a/tutorials/multiphase/cavitatingFoam/les/throttle/system/fvSchemes b/tutorials/multiphase/cavitatingFoam/les/throttle/system/fvSchemes
index ea089af844..6de4135a48 100644
--- a/tutorials/multiphase/cavitatingFoam/les/throttle/system/fvSchemes
+++ b/tutorials/multiphase/cavitatingFoam/les/throttle/system/fvSchemes
@@ -44,7 +44,7 @@ laplacianSchemes
     laplacian(nuf,rhoU) Gauss linear corrected;
     laplacian(muEff,U) Gauss linear corrected;
     laplacian(rrhoUAf,p) Gauss linear corrected;
-    laplacian(rUAf,p) Gauss linear corrected;
+    laplacian(rAUf,p) Gauss linear corrected;
     laplacian(DkEff,k) Gauss linear corrected;
     laplacian(1,p)  Gauss linear corrected;
 }
diff --git a/tutorials/multiphase/cavitatingFoam/les/throttle3D/system/fvSchemes b/tutorials/multiphase/cavitatingFoam/les/throttle3D/system/fvSchemes
index ea089af844..6de4135a48 100644
--- a/tutorials/multiphase/cavitatingFoam/les/throttle3D/system/fvSchemes
+++ b/tutorials/multiphase/cavitatingFoam/les/throttle3D/system/fvSchemes
@@ -44,7 +44,7 @@ laplacianSchemes
     laplacian(nuf,rhoU) Gauss linear corrected;
     laplacian(muEff,U) Gauss linear corrected;
     laplacian(rrhoUAf,p) Gauss linear corrected;
-    laplacian(rUAf,p) Gauss linear corrected;
+    laplacian(rAUf,p) Gauss linear corrected;
     laplacian(DkEff,k) Gauss linear corrected;
     laplacian(1,p)  Gauss linear corrected;
 }
diff --git a/tutorials/multiphase/cavitatingFoam/ras/throttle/system/fvSchemes b/tutorials/multiphase/cavitatingFoam/ras/throttle/system/fvSchemes
index 48be074223..f3c3cc3d66 100644
--- a/tutorials/multiphase/cavitatingFoam/ras/throttle/system/fvSchemes
+++ b/tutorials/multiphase/cavitatingFoam/ras/throttle/system/fvSchemes
@@ -45,7 +45,7 @@ laplacianSchemes
     laplacian(nuf,rhoU) Gauss linear corrected;
     laplacian(muEff,U) Gauss linear corrected;
     laplacian(rrhoUAf,p) Gauss linear corrected;
-    laplacian(rUAf,p) Gauss linear corrected;
+    laplacian(rAUf,p) Gauss linear corrected;
     laplacian(DomegaEff,omega) Gauss linear corrected;
     laplacian(DkEff,k) Gauss linear corrected;
     laplacian(1,p)  Gauss linear corrected;