diff --git a/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/CourantNo.H b/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/CourantNo.H
index 93c97ac4b..d7d1cc66d 100644
--- a/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/CourantNo.H
+++ b/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/CourantNo.H
@@ -39,12 +39,12 @@ if (mesh.nInternalFaces())
         fvc::surfaceSum(mag(phi))().internalField()
     );
 
-    forAll(fluid.phases(), phasei)
+    forAll(phases, phasei)
     {
         sumPhi = max
         (
             sumPhi,
-            fvc::surfaceSum(mag(fluid.phases()[phasei].phi()))().internalField()
+            fvc::surfaceSum(mag(phases[phasei].phi()))().internalField()
         );
     }
 
diff --git a/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/EEqns.H b/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/EEqns.H
index 15a8f6216..386331e0e 100644
--- a/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/EEqns.H
+++ b/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/EEqns.H
@@ -6,9 +6,9 @@
 
     phaseSystem::heatTransferTable& heatTransfer = heatTransferPtr();
 
-    forAll(fluid.phases(), phasei)
+    forAll(phases, phasei)
     {
-        phaseModel& phase = fluid.phases()[phasei];
+        phaseModel& phase = phases[phasei];
 
         const volScalarField& alpha = phase;
         const volScalarField& rho = phase.rho();
@@ -36,9 +36,9 @@
 
 fluid.correctThermo();
 
-forAll(fluid.phases(), phasei)
+forAll(phases, phasei)
 {
-    phaseModel& phase = fluid.phases()[phasei];
+    phaseModel& phase = phases[phasei];
 
     Info<< phase.name() << " min/max T "
         << min(phase.thermo().T()).value()
diff --git a/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/YEqns.H b/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/YEqns.H
index 5a44011fc..4094e185a 100644
--- a/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/YEqns.H
+++ b/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/YEqns.H
@@ -5,9 +5,9 @@
     phaseSystem::massTransferTable&
         massTransfer(massTransferPtr());
 
-    forAll(fluid.phases(), phasei)
+    forAll(phases, phasei)
     {
-        phaseModel& phase = fluid.phases()[phasei];
+        phaseModel& phase = phases[phasei];
 
         PtrList<volScalarField>& Y = phase.Y();
         const volScalarField& alpha = phase;
diff --git a/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/createFields.H b/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/createFields.H
index 7f89ea5f4..1375db8c4 100644
--- a/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/createFields.H
+++ b/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/createFields.H
@@ -8,6 +8,7 @@ autoPtr<multiphaseSystem> fluidPtr
     multiphaseSystem::New(mesh)
 );
 multiphaseSystem& fluid = fluidPtr();
+multiphaseSystem::phaseModelList& phases = fluid.phases();
 
 surfaceScalarField& phi = fluid.phi();
 
@@ -20,7 +21,7 @@ dimensionedScalar pMin
 
 #include "gh.H"
 
-volScalarField& p = fluid.phases().first().thermo().p();
+volScalarField& p = phases[0].thermo().p();
 
 Info<< "Reading field p_rgh\n" << endl;
 volScalarField p_rgh
diff --git a/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/pU/UEqns.H b/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/pU/UEqns.H
index 3b6ee9551..32cedb0fe 100644
--- a/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/pU/UEqns.H
+++ b/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/pU/UEqns.H
@@ -1,6 +1,6 @@
 Info<< "Constructing momentum equations" << endl;
 
-PtrList<fvVectorMatrix> UEqns(fluid.phases().size());
+PtrList<fvVectorMatrix> UEqns(phases.size());
 
 {
     autoPtr<phaseSystem::momentumTransferTable>
@@ -9,9 +9,9 @@ PtrList<fvVectorMatrix> UEqns(fluid.phases().size());
     phaseSystem::momentumTransferTable&
         momentumTransfer(momentumTransferPtr());
 
-    forAll(fluid.phases(), phasei)
+    forAll(phases, phasei)
     {
-        phaseModel& phase = fluid.phases()[phasei];
+        phaseModel& phase = phases[phasei];
 
         const volScalarField& alpha = phase;
         const volScalarField& rho = phase.rho();
diff --git a/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/pU/pEqn.H b/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/pU/pEqn.H
index ff759f42c..9fd728ec6 100644
--- a/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/pU/pEqn.H
+++ b/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/pU/pEqn.H
@@ -1,10 +1,10 @@
-PtrList<surfaceScalarField> alphafs(fluid.phases().size());
-PtrList<volScalarField> rAUs(fluid.phases().size());
-PtrList<surfaceScalarField> alpharAUfs(fluid.phases().size());
+PtrList<surfaceScalarField> alphafs(phases.size());
+PtrList<volScalarField> rAUs(phases.size());
+PtrList<surfaceScalarField> alpharAUfs(phases.size());
 
-forAll(fluid.phases(), phasei)
+forAll(phases, phasei)
 {
-    phaseModel& phase = fluid.phases()[phasei];
+    phaseModel& phase = phases[phasei];
     const volScalarField& alpha = phase;
 
     alphafs.set(phasei, fvc::interpolate(alpha).ptr());
@@ -35,13 +35,13 @@ forAll(fluid.phases(), phasei)
 }
 
 // Turbulent diffusion, particle-pressure, lift and wall-lubrication fluxes
-PtrList<surfaceScalarField> phiFs(fluid.phases().size());
+PtrList<surfaceScalarField> phiFs(phases.size());
 {
     autoPtr<PtrList<volVectorField> > Fs = fluid.Fs();
 
-    forAll(fluid.phases(), phasei)
+    forAll(phases, phasei)
     {
-        phaseModel& phase = fluid.phases()[phasei];
+        phaseModel& phase = phases[phasei];
 
         phiFs.set
         (
@@ -61,11 +61,11 @@ while (pimple.correct())
     // Update continuity errors due to temperature changes
     fluid.correct();
 
-    PtrList<volVectorField> HbyAs(fluid.phases().size());
+    PtrList<volVectorField> HbyAs(phases.size());
 
-    forAll(fluid.phases(), phasei)
+    forAll(phases, phasei)
     {
-        phaseModel& phase = fluid.phases()[phasei];
+        phaseModel& phase = phases[phasei];
         const volScalarField& alpha = phase;
 
         // Correct fixed-flux BCs to be consistent with the velocity BCs
@@ -99,10 +99,10 @@ while (pimple.correct())
         ghf*fvc::snGrad(rho)*mesh.magSf()
     );
 
-    PtrList<surfaceScalarField> phigs(fluid.phases().size());
-    forAll(fluid.phases(), phasei)
+    PtrList<surfaceScalarField> phigs(phases.size());
+    forAll(phases, phasei)
     {
-        phaseModel& phase = fluid.phases()[phasei];
+        phaseModel& phase = phases[phasei];
 
         phigs.set
         (
@@ -118,7 +118,7 @@ while (pimple.correct())
         );
     }
 
-    PtrList<surfaceScalarField> phiHbyAs(fluid.phases().size());
+    PtrList<surfaceScalarField> phiHbyAs(phases.size());
 
     surfaceScalarField phiHbyA
     (
@@ -134,9 +134,9 @@ while (pimple.correct())
         dimensionedScalar("phiHbyA", dimArea*dimVelocity, 0)
     );
 
-    forAll(fluid.phases(), phasei)
+    forAll(phases, phasei)
     {
-        phaseModel& phase = fluid.phases()[phasei];
+        phaseModel& phase = phases[phasei];
         const volScalarField& alpha = phase;
 
         // ddtPhiCorr filter -- only apply in pure(ish) phases
@@ -227,7 +227,7 @@ while (pimple.correct())
         dimensionedScalar("rAUf", dimensionSet(-1, 3, 1, 0, 0), 0)
     );
 
-    forAll(fluid.phases(), phasei)
+    forAll(phases, phasei)
     {
         rAUf += alphafs[phasei]*alpharAUfs[phasei];
     }
@@ -238,9 +238,9 @@ while (pimple.correct())
     {
         surfaceScalarField::GeometricBoundaryField phib(phi.boundaryField());
         phib = 0;
-        forAll(fluid.phases(), phasei)
+        forAll(phases, phasei)
         {
-            phaseModel& phase = fluid.phases()[phasei];
+            phaseModel& phase = phases[phasei];
             phib += alphafs[phasei].boundaryField()*phase.phi().boundaryField();
         }
 
@@ -253,10 +253,10 @@ while (pimple.correct())
         );
     }
 
-    PtrList<fvScalarMatrix> pEqnComps(fluid.phases().size());
-    forAll(fluid.phases(), phasei)
+    PtrList<fvScalarMatrix> pEqnComps(phases.size());
+    forAll(phases, phasei)
     {
-        phaseModel& phase = fluid.phases()[phasei];
+        phaseModel& phase = phases[phasei];
 
         if (phase.compressible())
         {
@@ -336,9 +336,9 @@ while (pimple.correct())
         {
             fvScalarMatrix pEqn(pEqnIncomp);
 
-            forAll(fluid.phases(), phasei)
+            forAll(phases, phasei)
             {
-                phaseModel& phase = fluid.phases()[phasei];
+                phaseModel& phase = phases[phasei];
 
                 if (phase.compressible())
                 {
@@ -360,9 +360,9 @@ while (pimple.correct())
 
             surfaceScalarField mSfGradp("mSfGradp", pEqnIncomp.flux()/rAUf);
 
-            forAll(fluid.phases(), phasei)
+            forAll(phases, phasei)
             {
-                phaseModel& phase = fluid.phases()[phasei];
+                phaseModel& phase = phases[phasei];
 
                 phase.phi() = phiHbyAs[phasei] + alpharAUfs[phasei]*mSfGradp;
 
@@ -378,9 +378,9 @@ while (pimple.correct())
 
             mSfGradp = pEqnIncomp.flux()/rAUf;
 
-            forAll(fluid.phases(), phasei)
+            forAll(phases, phasei)
             {
-                phaseModel& phase = fluid.phases()[phasei];
+                phaseModel& phase = phases[phasei];
 
                 phase.U() =
                     HbyAs[phasei]
@@ -403,9 +403,9 @@ while (pimple.correct())
     p_rgh = p - rho*gh;
 
     // Update densities from change in p_rgh
-    forAll(fluid.phases(), phasei)
+    forAll(phases, phasei)
     {
-        phaseModel& phase = fluid.phases()[phasei];
+        phaseModel& phase = phases[phasei];
         phase.rho()() += phase.thermo().psi()*(p_rgh - p_rgh_0);
     }
 
diff --git a/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/setRDeltaT.H b/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/setRDeltaT.H
index aae1d9072..9e9ca5423 100644
--- a/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/setRDeltaT.H
+++ b/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/setRDeltaT.H
@@ -12,9 +12,9 @@
 
     surfaceScalarField maxPhi("maxPhi", phi);
 
-    forAll(fluid.phases(), phasei)
+    forAll(phases, phasei)
     {
-        maxPhi = max(maxPhi, mag(fluid.phases()[phasei].phi()));
+        maxPhi = max(maxPhi, mag(phases[phasei].phi()));
     }
 
     // Set the reciprocal time-step from the local Courant number