diff --git a/applications/solvers/combustion/XiFoam/UEqn.H b/applications/solvers/combustion/XiFoam/UEqn.H
index 2a7753e14e1f5efcc31edcf3f2e2e721d5dde045..9697c6e1ed2d0753f9c0803081cc6929050ef446 100644
--- a/applications/solvers/combustion/XiFoam/UEqn.H
+++ b/applications/solvers/combustion/XiFoam/UEqn.H
@@ -7,6 +7,8 @@
         rho*g
     );
 
+    UEqn.relax();
+
     if (momentumPredictor)
     {
         solve(UEqn == -fvc::grad(p));
diff --git a/applications/solvers/combustion/XiFoam/XiFoam.C b/applications/solvers/combustion/XiFoam/XiFoam.C
index 4c544b433e3b0aba41261d237fc198bc616412a4..6e804be977b1eec581016550cfb0c21cb44113d2 100644
--- a/applications/solvers/combustion/XiFoam/XiFoam.C
+++ b/applications/solvers/combustion/XiFoam/XiFoam.C
@@ -87,12 +87,12 @@ int main(int argc, char *argv[])
         runTime++;
         Info<< "Time = " << runTime.timeName() << nl << endl;
 
-        #include "rhoEqn.H"
-        #include "UEqn.H"
-
-        // --- PISO loop
-        for (int corr=1; corr<=nCorr; corr++)
+        // --- Pressure-velocity PIMPLE corrector loop
+        for (int oCorr=0; oCorr<nOuterCorr; oCorr++)
         {
+            #include "rhoEqn.H"
+            #include "UEqn.H"
+
             #include "ftEqn.H"
             #include "bEqn.H"
             #include "huEqn.H"
@@ -103,10 +103,14 @@ int main(int argc, char *argv[])
                 hu == h;
             }
 
-            #include "pEqn.H"
-        }
+            // --- PISO loop
+            for (int corr=1; corr<=nCorr; corr++)
+            {
+                #include "pEqn.H"
+            }
 
-        turbulence->correct();
+            turbulence->correct();
+        }
 
         rho = thermo.rho();
 
diff --git a/applications/solvers/combustion/XiFoam/bEqn.H b/applications/solvers/combustion/XiFoam/bEqn.H
index 33ef24bfe655b2d587b7e403fbd8ee1f5d25b736..d06ec2e5f0a1252c1262bcdc702c65ddb759b485 100644
--- a/applications/solvers/combustion/XiFoam/bEqn.H
+++ b/applications/solvers/combustion/XiFoam/bEqn.H
@@ -52,7 +52,7 @@ if (ign.ignited())
       + mvConvection->fvmDiv(phi, b)
       + fvm::div(phiSt, b, "div(phiSt,b)")
       - fvm::Sp(fvc::div(phiSt), b)
-      - fvm::laplacian(turbulence->muEff(), b)
+      - fvm::laplacian(turbulence->alphaEff(), b)
     );
 
 
@@ -90,7 +90,7 @@ if (ign.ignited())
     // ~~~~~~~~~~~~~~~~~
     surfaceScalarField phiXi =
         phiSt
-      - fvc::interpolate(fvc::laplacian(turbulence->muEff(), b)/mgb)*nf
+      - fvc::interpolate(fvc::laplacian(turbulence->alphaEff(), b)/mgb)*nf
       + fvc::interpolate(rho)*fvc::interpolate(Su*(1.0/Xi - Xi))*nf;
 
 
diff --git a/applications/solvers/combustion/XiFoam/ftEqn.H b/applications/solvers/combustion/XiFoam/ftEqn.H
index 519cbd7cbe4bfa08623d7a62678ebf5188bfee1b..46d7aeae8228a9d47d3cc28eac64c8b5a65cfb3e 100644
--- a/applications/solvers/combustion/XiFoam/ftEqn.H
+++ b/applications/solvers/combustion/XiFoam/ftEqn.H
@@ -1,7 +1,7 @@
 tmp<fv::convectionScheme<scalar> > mvConvection
 (
     fv::convectionScheme<scalar>::New
-    (    
+    (
         mesh,
         fields,
         phi,
@@ -17,6 +17,6 @@ if (composition.contains("ft"))
     (
         fvm::ddt(rho, ft)
       + mvConvection->fvmDiv(phi, ft)
-      - fvm::laplacian(turbulence->muEff(), ft)
+      - fvm::laplacian(turbulence->alphaEff(), ft)
     );
 }
diff --git a/applications/solvers/combustion/XiFoam/hEqn.H b/applications/solvers/combustion/XiFoam/hEqn.H
index ebce30e24e50dee92824de875abfa138f0abf3b1..513ae604419f28b96bbf455e8ed0bc1e75d258bb 100644
--- a/applications/solvers/combustion/XiFoam/hEqn.H
+++ b/applications/solvers/combustion/XiFoam/hEqn.H
@@ -1,5 +1,5 @@
 {
-    solve
+    fvScalarMatrix hEqn
     (
         fvm::ddt(rho, h)
       + mvConvection->fvmDiv(phi, h)
@@ -8,5 +8,8 @@
         DpDt
     );
 
+    hEqn.relax();
+    hEqn.solve();
+
     thermo.correct();
 }