diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/EEqns.H b/applications/solvers/multiphase/twoPhaseEulerFoam/EEqns.H
index 3127450f0090dbe406d994f192ceb674be5f5936..f53973f0cc514263c30fcc3418eebd63da8ef55c 100644
--- a/applications/solvers/multiphase/twoPhaseEulerFoam/EEqns.H
+++ b/applications/solvers/multiphase/twoPhaseEulerFoam/EEqns.H
@@ -23,7 +23,7 @@
         )/rho1
 
       //***HGW- fvm::laplacian(alpha1*turbulence1->alphaEff(), he1)
-      - fvm::laplacian(alpha1*turbulence1->nuEff(), he1)
+      - fvm::laplacian(alpha1*phase1.turbulence().nuEff(), he1)
      ==
         heatTransferCoeff*(thermo2.T() - thermo1.T())/rho1
       + heatTransferCoeff*he1/Cpv1/rho1
@@ -46,7 +46,7 @@
         )/rho2
 
       //***HGW- fvm::laplacian(alpha2*turbulence2->alphaEff(), he2)
-      - fvm::laplacian(alpha2*turbulence2->nuEff(), he2)
+      - fvm::laplacian(alpha2*phase2.turbulence().nuEff(), he2)
      ==
         heatTransferCoeff*(thermo1.T() - thermo2.T())/rho2
       + heatTransferCoeff*he2/Cpv2/rho2
diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/UEqns.H b/applications/solvers/multiphase/twoPhaseEulerFoam/UEqns.H
index dea48be335460a4e698b13fa867467eedb0af76c..fad73cfeb8dd8ec1a517f164d2cfdf3f20ba0f68 100644
--- a/applications/solvers/multiphase/twoPhaseEulerFoam/UEqns.H
+++ b/applications/solvers/multiphase/twoPhaseEulerFoam/UEqns.H
@@ -26,7 +26,7 @@ volScalarField dragCoeff(fluid.dragCoeff());
               - fvm::Sp(fvc::div(phi1), U1)
             )
 
-          + turbulence1->divDevReff(U1)
+          + phase1.turbulence().divDevReff(U1)
          ==
           - fvm::Sp(dragCoeff/rho1, U1)
           - alpha1*alpha2/rho1*(liftForce - fluid.Cvm()*rho2*DDtU2)
@@ -50,7 +50,7 @@ volScalarField dragCoeff(fluid.dragCoeff());
               + fvm::div(phi2, U2)
               - fvm::Sp(fvc::div(phi2), U2)
             )
-          + turbulence2->divDevReff(U2)
+          + phase2.turbulence().divDevReff(U2)
          ==
           - fvm::Sp(dragCoeff/rho2, U2)
           + alpha1*alpha2/rho2*(liftForce + fluid.Cvm()*rho2*DDtU1)
diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/alphaEqn.H b/applications/solvers/multiphase/twoPhaseEulerFoam/alphaEqn.H
deleted file mode 100644
index be420a36e1f9e01b96cb07097e275f919e632012..0000000000000000000000000000000000000000
--- a/applications/solvers/multiphase/twoPhaseEulerFoam/alphaEqn.H
+++ /dev/null
@@ -1,165 +0,0 @@
-{
-    word alphaScheme("div(phi," + alpha1.name() + ')');
-    word alpharScheme("div(phir," + alpha1.name() + ')');
-
-    alpha1.correctBoundaryConditions();
-
-    surfaceScalarField phic("phic", phi);
-    surfaceScalarField phir("phir", phi1 - phi2);
-
-    surfaceScalarField alpha1f(fvc::interpolate(max(alpha1, scalar(0))));
-
-    tmp<surfaceScalarField> pPrimeByA;
-
-    if (implicitPhasePressure)
-    {
-        pPrimeByA =
-            fvc::interpolate((1.0/rho1)*rAU1*turbulence1().pPrime())
-          + fvc::interpolate((1.0/rho2)*rAU2*turbulence2().pPrime());
-
-        surfaceScalarField phiP
-        (
-            pPrimeByA()*fvc::snGrad(alpha1, "bounded")*mesh.magSf()
-        );
-
-        phic += alpha1f*phiP;
-        phir += phiP;
-    }
-
-    for (int acorr=0; acorr<nAlphaCorr; acorr++)
-    {
-        volScalarField::DimensionedInternalField Sp
-        (
-            IOobject
-            (
-                "Sp",
-                runTime.timeName(),
-                mesh
-            ),
-            mesh,
-            dimensionedScalar("Sp", dgdt.dimensions(), 0.0)
-        );
-
-        volScalarField::DimensionedInternalField Su
-        (
-            IOobject
-            (
-                "Su",
-                runTime.timeName(),
-                mesh
-            ),
-            // Divergence term is handled explicitly to be
-            // consistent with the explicit transport solution
-            fvc::div(phi)*min(alpha1, scalar(1))
-        );
-
-        forAll(dgdt, celli)
-        {
-            if (dgdt[celli] > 0.0 && alpha1[celli] > 0.0)
-            {
-                Sp[celli] -= dgdt[celli]*alpha1[celli];
-                Su[celli] += dgdt[celli]*alpha1[celli];
-            }
-            else if (dgdt[celli] < 0.0 && alpha1[celli] < 1.0)
-            {
-                Sp[celli] += dgdt[celli]*(1.0 - alpha1[celli]);
-            }
-        }
-
-        dimensionedScalar totalDeltaT = runTime.deltaT();
-        if (nAlphaSubCycles > 1)
-        {
-            alphaPhi1 = dimensionedScalar("0", alphaPhi1.dimensions(), 0);
-        }
-
-        for
-        (
-            subCycle<volScalarField> alphaSubCycle(alpha1, nAlphaSubCycles);
-            !(++alphaSubCycle).end();
-        )
-        {
-            surfaceScalarField alphaPhic1
-            (
-                fvc::flux
-                (
-                    phic,
-                    alpha1,
-                    alphaScheme
-                )
-              + fvc::flux
-                (
-                    -fvc::flux(-phir, scalar(1) - alpha1, alpharScheme),
-                    alpha1,
-                    alpharScheme
-                )
-            );
-
-            // Ensure that the flux at inflow BCs is preserved
-            forAll(alphaPhic1.boundaryField(), patchi)
-            {
-                fvsPatchScalarField& alphaPhic1p =
-                    alphaPhic1.boundaryField()[patchi];
-
-                if (!alphaPhic1p.coupled())
-                {
-                    const scalarField& phi1p = phi1.boundaryField()[patchi];
-                    const scalarField& alpha1p = alpha1.boundaryField()[patchi];
-
-                    forAll(alphaPhic1p, facei)
-                    {
-                        if (phi1p[facei] < 0)
-                        {
-                            alphaPhic1p[facei] = alpha1p[facei]*phi1p[facei];
-                        }
-                    }
-                }
-            }
-
-            MULES::explicitSolve
-            (
-                geometricOneField(),
-                alpha1,
-                phi,
-                alphaPhic1,
-                Sp,
-                Su,
-                1,
-                0
-            );
-
-            if (nAlphaSubCycles > 1)
-            {
-                alphaPhi1 += (runTime.deltaT()/totalDeltaT)*alphaPhic1;
-            }
-            else
-            {
-                alphaPhi1 = alphaPhic1;
-            }
-        }
-
-        if (implicitPhasePressure)
-        {
-            fvScalarMatrix alpha1Eqn
-            (
-                fvm::ddt(alpha1) - fvc::ddt(alpha1)
-              - fvm::laplacian(alpha1f*pPrimeByA, alpha1, "bounded")
-            );
-
-            alpha1Eqn.relax();
-            alpha1Eqn.solve();
-
-            alphaPhi1 += alpha1Eqn.flux();
-        }
-
-        alphaPhi2 = phi - alphaPhi1;
-        alpha2 = scalar(1) - alpha1;
-
-        Info<< "Dispersed phase volume fraction = "
-            << alpha1.weightedAverage(mesh.V()).value()
-            << "  Min(alpha1) = " << min(alpha1).value()
-            << "  Max(alpha1) = " << max(alpha1).value()
-            << endl;
-    }
-}
-
-rho = fluid.rho();
diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/createFields.H b/applications/solvers/multiphase/twoPhaseEulerFoam/createFields.H
index f279445cf639bc7992037f31a6f4b92ed487c714..c3da8d972a299a103d802aacf2d74a964a9711d0 100644
--- a/applications/solvers/multiphase/twoPhaseEulerFoam/createFields.H
+++ b/applications/solvers/multiphase/twoPhaseEulerFoam/createFields.H
@@ -10,19 +10,13 @@
 
     volVectorField& U1 = phase1.U();
     surfaceScalarField& phi1 = phase1.phi();
-    surfaceScalarField alphaPhi1
-    (
-        IOobject::groupName("alphaPhi", phase1.name()),
-        fvc::interpolate(alpha1)*phi1
-    );
+    surfaceScalarField& alphaPhi1 = phase1.phiAlpha();
 
     volVectorField& U2 = phase2.U();
     surfaceScalarField& phi2 = phase2.phi();
-    surfaceScalarField alphaPhi2
-    (
-        IOobject::groupName("alphaPhi", phase2.name()),
-        fvc::interpolate(alpha2)*phi2
-    );
+    surfaceScalarField& alphaPhi2 = phase2.phiAlpha();
+
+    surfaceScalarField& phi = fluid.phi();
 
     dimensionedScalar pMin
     (
@@ -55,19 +49,6 @@
         fluid.U()
     );
 
-    surfaceScalarField phi
-    (
-        IOobject
-        (
-            "phi",
-            runTime.timeName(),
-            mesh,
-            IOobject::NO_READ,
-            IOobject::AUTO_WRITE
-        ),
-        fluid.phi()
-    );
-
     volScalarField rho
     (
         IOobject
@@ -133,13 +114,6 @@
     scalar pRefValue = 0.0;
     setRefCell(p, mesh.solutionDict().subDict("PIMPLE"), pRefCell, pRefValue);
 
-
-    volScalarField dgdt
-    (
-        pos(alpha2)*fvc::div(phi)/max(alpha2, scalar(0.0001))
-    );
-
-
     Info<< "Creating field dpdt\n" << endl;
     volScalarField dpdt
     (
@@ -157,29 +131,3 @@
     Info<< "Creating field kinetic energy K\n" << endl;
     volScalarField K1(IOobject::groupName("K", phase1.name()), 0.5*magSqr(U1));
     volScalarField K2(IOobject::groupName("K", phase2.name()), 0.5*magSqr(U2));
-
-    autoPtr<PhaseIncompressibleTurbulenceModel<phaseModel> >
-    turbulence1
-    (
-        PhaseIncompressibleTurbulenceModel<phaseModel>::New
-        (
-            alpha1,
-            U1,
-            alphaPhi1,
-            phi1,
-            phase1
-        )
-    );
-
-    autoPtr<PhaseIncompressibleTurbulenceModel<phaseModel> >
-    turbulence2
-    (
-        PhaseIncompressibleTurbulenceModel<phaseModel>::New
-        (
-            alpha2,
-            U2,
-            alphaPhi2,
-            phi2,
-            phase2
-        )
-    );
diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/pEqn.H b/applications/solvers/multiphase/twoPhaseEulerFoam/pEqn.H
index 9e99c5cc61bfd9daec53705af3f31add2ec095ab..b097abf7451d8c2dc6c500e6dd8d91ad4d281e2e 100644
--- a/applications/solvers/multiphase/twoPhaseEulerFoam/pEqn.H
+++ b/applications/solvers/multiphase/twoPhaseEulerFoam/pEqn.H
@@ -31,7 +31,7 @@
     surfaceScalarField phiP1
     (
         "phiP1",
-        fvc::interpolate((1.0/rho1)*rAU1*turbulence1().pPrime())
+        fvc::interpolate((1.0/rho1)*rAU1*phase1.turbulence().pPrime())
        *fvc::snGrad(alpha1)*mesh.magSf()
     );
 
@@ -39,7 +39,7 @@
     surfaceScalarField phiP2
     (
         "phiP2",
-        fvc::interpolate((1.0/rho2)*rAU2*turbulence2().pPrime())
+        fvc::interpolate((1.0/rho2)*rAU2*phase2.turbulence().pPrime())
        *fvc::snGrad(alpha2)*mesh.magSf()
     );
 
@@ -199,7 +199,7 @@
 
             phi = alpha1f*phi1 + alpha2f*phi2;
 
-            dgdt =
+            fluid.dgdt() =
             (
                 pos(alpha2)*(pEqnComp2 & p)/rho2
               - pos(alpha1)*(pEqnComp1 & p)/rho1
diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/readTwoPhaseEulerFoamControls.H b/applications/solvers/multiphase/twoPhaseEulerFoam/readTwoPhaseEulerFoamControls.H
deleted file mode 100644
index 29353a8fa1deea189e0fa017c28ba5685013245c..0000000000000000000000000000000000000000
--- a/applications/solvers/multiphase/twoPhaseEulerFoam/readTwoPhaseEulerFoamControls.H
+++ /dev/null
@@ -1,7 +0,0 @@
-    #include "readTimeControls.H"
-    #include "alphaControls.H"
-
-    Switch implicitPhasePressure
-    (
-        alphaControls.lookupOrDefault<Switch>("implicitPhasePressure", false)
-    );
diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseEulerFoam.C b/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseEulerFoam.C
index 1e7cbafea75b13045ce2df010d2a4f164bb7f3ad..69682a9a3d52f6209b787b43d6e0eca11b4b1aea 100644
--- a/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseEulerFoam.C
+++ b/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseEulerFoam.C
@@ -31,13 +31,10 @@ Description
 \*---------------------------------------------------------------------------*/
 
 #include "fvCFD.H"
-#include "MULES.H"
-#include "subCycle.H"
-#include "rhoThermo.H"
 #include "twoPhaseSystem.H"
+#include "PhaseIncompressibleTurbulenceModel.H"
 #include "dragModel.H"
 #include "heatTransferModel.H"
-#include "PhaseIncompressibleTurbulenceModel.H"
 #include "pimpleControl.H"
 #include "IOMRFZoneList.H"
 #include "fixedFluxPressureFvPatchScalarField.H"
@@ -66,7 +63,7 @@ int main(int argc, char *argv[])
 
     while (runTime.run())
     {
-        #include "readTwoPhaseEulerFoamControls.H"
+        #include "readTimeControls.H"
         #include "CourantNos.H"
         #include "setDeltaT.H"
 
@@ -76,7 +73,10 @@ int main(int argc, char *argv[])
         // --- Pressure-velocity PIMPLE corrector loop
         while (pimple.loop())
         {
-            #include "alphaEqn.H"
+            fluid.solve();
+            rho = fluid.rho();
+            fluid.correct();
+
             #include "EEqns.H"
             #include "UEqns.H"
 
@@ -90,8 +90,7 @@ int main(int argc, char *argv[])
 
             if (pimple.turbCorr())
             {
-                turbulence1->correct();
-                turbulence2->correct();
+                fluid.correctTurbulence();
             }
         }
 
diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/Make/files b/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/Make/files
index fb49c3cef761ab0d94d4b702a807f0378cd0cc74..0de9e8939cb8c1efbe265dc11d65836963512eb6 100644
--- a/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/Make/files
+++ b/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/Make/files
@@ -4,6 +4,12 @@ diameterModels/diameterModel/newDiameterModel.C
 diameterModels/constantDiameter/constantDiameter.C
 diameterModels/isothermalDiameter/isothermalDiameter.C
 
+diameterModels/IATE/IATE.C
+diameterModels/IATE/IATEsources/IATEsource/IATEsource.C
+diameterModels/IATE/IATEsources/wakeEntrainmentCoalescence/wakeEntrainmentCoalescence.C
+diameterModels/IATE/IATEsources/turbulentBreakUp/turbulentBreakUp.C
+diameterModels/IATE/IATEsources/randomCoalescence/randomCoalescence.C
+
 twoPhaseSystem.C
 
 LIB = $(FOAM_LIBBIN)/libcompressibleTwoPhaseSystem
diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/Make/options b/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/Make/options
index ab3c396f5796ea50f8afd47eeb3a4c6d5c950316..eec60d23e7ea35a3951697ccf22f57070385b9d1 100644
--- a/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/Make/options
+++ b/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/Make/options
@@ -3,7 +3,10 @@ EXE_INC = \
     -I../interfacialModels/lnInclude \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
-    -I$(LIB_SRC)/transportModels/incompressible/lnInclude
+    -I$(LIB_SRC)/transportModels/incompressible/lnInclude \
+    -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \
+    -I$(LIB_SRC)/TurbulenceModels/incompressible/lnInclude \
+    -I$(LIB_SRC)/TurbulenceModels/phaseIncompressible/lnInclude
 
 LIB_LIBS = \
     -lincompressibleTransportModels \
diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/diameterModels/IATE/IATE.C b/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/diameterModels/IATE/IATE.C
new file mode 100644
index 0000000000000000000000000000000000000000..ef4bc12d0c79f51af6b87d92c32d543474d725f0
--- /dev/null
+++ b/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/diameterModels/IATE/IATE.C
@@ -0,0 +1,188 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "IATE.H"
+#include "IATEsource.H"
+#include "twoPhaseSystem.H"
+#include "fvmDdt.H"
+#include "fvmDiv.H"
+#include "fvmSup.H"
+#include "fvcDdt.H"
+#include "fvcDiv.H"
+#include "fvcAverage.H"
+#include "mathematicalConstants.H"
+#include "fundamentalConstants.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace diameterModels
+{
+    defineTypeNameAndDebug(IATE, 0);
+
+    addToRunTimeSelectionTable
+    (
+        diameterModel,
+        IATE,
+        dictionary
+    );
+}
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::diameterModels::IATE::IATE
+(
+    const dictionary& diameterProperties,
+    const phaseModel& phase
+)
+:
+    diameterModel(diameterProperties, phase),
+    kappai_
+    (
+        IOobject
+        (
+            IOobject::groupName("kappai", phase.name()),
+            phase_.U().time().timeName(),
+            phase_.U().mesh(),
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        phase_.U().mesh()
+    ),
+    dMax_("dMax", dimLength, diameterProperties_.lookup("dMax")),
+    dMin_("dMin", dimLength, diameterProperties_.lookup("dMin")),
+    d_
+    (
+        IOobject
+        (
+            IOobject::groupName("d", phase.name()),
+            phase_.U().time().timeName(),
+            phase_.U().mesh(),
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        dsm()
+    ),
+    sources_
+    (
+        diameterProperties_.lookup("sources"),
+        IATEsource::iNew(*this)
+    )
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::diameterModels::IATE::~IATE()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+Foam::tmp<Foam::volScalarField> Foam::diameterModels::IATE::dsm() const
+{
+    return max(6/max(kappai_, 6/dMax_), dMin_);
+}
+
+// Placeholder for the nucleation/condensation model
+// Foam::tmp<Foam::volScalarField> Foam::diameterModels::IATE::Rph() const
+// {
+//     const volScalarField& T = phase_.thermo().T();
+//     const volScalarField& p = phase_.thermo().p();
+//
+//     scalar A, B, C, sigma, vm, Rph;
+//
+//     volScalarField ps(1e5*pow(10, A - B/(T + C)));
+//     volScalarField Dbc
+//     (
+//         4*sigma*vm/(constant::physicoChemical::k*T*log(p/ps))
+//     );
+//
+//     return constant::mathematical::pi*sqr(Dbc)*Rph;
+// }
+
+void Foam::diameterModels::IATE::correct()
+{
+    // Initialise the accumulated source term to the dilatation effect
+    volScalarField R
+    (
+        (
+            (2.0/3.0)
+           /max
+            (
+                fvc::average(phase_ + phase_.oldTime()),
+                phase_.fluid().residualPhaseFraction()
+            )
+        )
+       *(fvc::ddt(phase_) + fvc::div(phase_.phiAlpha()))
+    );
+
+    // Accumulate the run-time selectable sources
+    forAll(sources_, j)
+    {
+        R -= sources_[j].R();
+    }
+
+    // Construct the interfacial curvature equation
+    fvScalarMatrix kappaiEqn
+    (
+        fvm::ddt(kappai_) + fvm::div(phase_.phi(), kappai_)
+      - fvm::Sp(fvc::div(phase_.phi()), kappai_)
+     ==
+      - fvm::SuSp(R, kappai_)
+    //+ Rph() // Omit the nucleation/condensation term
+    );
+
+    kappaiEqn.relax();
+    kappaiEqn.solve();
+
+    // Update the Sauter-mean diameter
+    d_ = dsm();
+}
+
+
+bool Foam::diameterModels::IATE::read(const dictionary& phaseProperties)
+{
+    diameterModel::read(phaseProperties);
+
+    diameterProperties_.lookup("dMax") >> dMax_;
+    diameterProperties_.lookup("dMin") >> dMin_;
+
+    // Re-create all the sources updating number, type and coefficients
+    PtrList<IATEsource>
+    (
+        diameterProperties_.lookup("sources"),
+        IATEsource::iNew(*this)
+    ).transfer(sources_);
+
+    return true;
+}
+
+
+// ************************************************************************* //
diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/diameterModels/IATE/IATE.H b/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/diameterModels/IATE/IATE.H
new file mode 100644
index 0000000000000000000000000000000000000000..f599fb9959f858e422c2ad11b570413dbe9b1afd
--- /dev/null
+++ b/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/diameterModels/IATE/IATE.H
@@ -0,0 +1,154 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::diameterModels::IATE
+
+Description
+    IATE (Interfacial Area Transport Equation) bubble diameter model.
+
+    Solves for the interfacial curvature per unit volume of the phase rather
+    than interfacial area per unit volume to avoid stability issues relating to
+    the consistency requirements between the phase fraction and interfacial area
+    per unit volume.  In every other respect this model is as presented in the
+    paper:
+
+    \verbatim
+        "Development of Interfacial Area Transport Equation"
+        M. Ishii,
+        S. Kim,
+        J Kelly,
+        Nuclear Engineering and Technology, Vol.37 No.6 December 2005
+    \endverbatim
+
+SourceFiles
+    IATE.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef IATE_H
+#define IATE_H
+
+#include "diameterModel.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace diameterModels
+{
+
+// Forward declaration of classes
+class IATEsource;
+
+/*---------------------------------------------------------------------------*\
+                           Class IATE Declaration
+\*---------------------------------------------------------------------------*/
+
+class IATE
+:
+    public diameterModel
+{
+    // Private data
+
+        //- Interfacial curvature (alpha*interfacial area)
+        volScalarField kappai_;
+
+        //- Maximum diameter used for stabilisation in the limit kappai->0
+        dimensionedScalar dMax_;
+
+        //- Minimum diameter used for stabilisation in the limit kappai->inf
+        dimensionedScalar dMin_;
+
+        //- The Sauter-mean diameter of the phase
+        volScalarField d_;
+
+        //- IATE sources
+        PtrList<IATEsource> sources_;
+
+
+    // Private member functions
+
+        tmp<volScalarField> dsm() const;
+
+
+public:
+
+    friend class IATEsource;
+
+    //- Runtime type information
+    TypeName("IATE");
+
+
+    // Constructors
+
+        //- Construct from components
+        IATE
+        (
+            const dictionary& diameterProperties,
+            const phaseModel& phase
+        );
+
+
+    //- Destructor
+    virtual ~IATE();
+
+
+    // Member Functions
+
+        //- Return the interfacial curvature
+        const volScalarField& kappai() const
+        {
+            return kappai_;
+        }
+
+        //- Return the interfacial area
+        tmp<volScalarField> a() const
+        {
+            return phase_*kappai_;
+        }
+
+        //- Return the Sauter-mean diameter
+        virtual tmp<volScalarField> d() const
+        {
+            return d_;
+        }
+
+        //- Correct the diameter field
+        virtual void correct();
+
+        //- Read phaseProperties dictionary
+        virtual bool read(const dictionary& phaseProperties);
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace diameterModels
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/diameterModels/IATE/IATEsources/IATEsource/IATEsource.C b/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/diameterModels/IATE/IATEsources/IATEsource/IATEsource.C
new file mode 100644
index 0000000000000000000000000000000000000000..4d8343133128eec13da6f2bc46b0dd0414f0bed5
--- /dev/null
+++ b/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/diameterModels/IATE/IATEsources/IATEsource/IATEsource.C
@@ -0,0 +1,147 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "IATEsource.H"
+#include "twoPhaseSystem.H"
+#include "fvMatrix.H"
+#include "PhaseIncompressibleTurbulenceModel.H"
+#include "uniformDimensionedFields.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace diameterModels
+{
+    defineTypeNameAndDebug(IATEsource, 0);
+    defineRunTimeSelectionTable(IATEsource, dictionary);
+}
+}
+
+
+// * * * * * * * * * * * * * * * * * Selectors * * * * * * * * * * * * * * * //
+
+Foam::autoPtr<Foam::diameterModels::IATEsource>
+Foam::diameterModels::IATEsource::New
+(
+    const word& type,
+    const IATE& iate,
+    const dictionary& dict
+)
+{
+    dictionaryConstructorTable::iterator cstrIter =
+        dictionaryConstructorTablePtr_->find(type);
+
+    if (cstrIter == dictionaryConstructorTablePtr_->end())
+    {
+        FatalErrorIn
+        (
+            "IATEsource::New"
+            "(const word& type, const IATE&, const dictionary&)"
+        )   << "Unknown IATE source type "
+            << type << nl << nl
+            << "Valid IATE source types : " << endl
+            << dictionaryConstructorTablePtr_->sortedToc()
+            << exit(FatalError);
+    }
+
+    return autoPtr<IATEsource>(cstrIter()(iate, dict));
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+Foam::tmp<Foam::volScalarField> Foam::diameterModels::IATEsource::Ur() const
+{
+    const uniformDimensionedVectorField& g =
+        phase().U().db().lookupObject<uniformDimensionedVectorField>("g");
+
+    return
+        sqrt(2.0)
+       *pow025
+        (
+            fluid().sigma()*mag(g)
+           *(otherPhase().rho() - phase().rho())
+           /sqr(otherPhase().rho())
+        )
+       *pow(max(1 - phase(), scalar(0)), 1.75);
+}
+
+Foam::tmp<Foam::volScalarField> Foam::diameterModels::IATEsource::Ut() const
+{
+    return sqrt(2*otherPhase().turbulence().k());
+}
+
+Foam::tmp<Foam::volScalarField> Foam::diameterModels::IATEsource::Re() const
+{
+    return max(Ur()*phase().d()/otherPhase().nu(), scalar(1.0e-3));
+}
+
+Foam::tmp<Foam::volScalarField> Foam::diameterModels::IATEsource::CD() const
+{
+    const volScalarField Eo(this->Eo());
+    const volScalarField Re(this->Re());
+
+    return
+        max
+        (
+            min
+            (
+                (16/Re)*(1 + 0.15*pow(Re, 0.687)),
+                48/Re
+            ),
+            8*Eo/(3*(Eo + 4))
+        );
+}
+
+Foam::tmp<Foam::volScalarField> Foam::diameterModels::IATEsource::Mo() const
+{
+    const uniformDimensionedVectorField& g =
+        phase().U().db().lookupObject<uniformDimensionedVectorField>("g");
+
+    return
+        mag(g)*pow4(otherPhase().nu())*sqr(otherPhase().rho())
+       *(otherPhase().rho() - phase().rho())
+       /pow3(fluid().sigma());
+}
+
+Foam::tmp<Foam::volScalarField> Foam::diameterModels::IATEsource::Eo() const
+{
+    const uniformDimensionedVectorField& g =
+        phase().U().db().lookupObject<uniformDimensionedVectorField>("g");
+
+    return
+        mag(g)*sqr(phase().d())
+       *(otherPhase().rho() - phase().rho())
+       /fluid().sigma();
+}
+
+Foam::tmp<Foam::volScalarField> Foam::diameterModels::IATEsource::We() const
+{
+    return otherPhase().rho()*sqr(Ur())*phase().d()/fluid().sigma();
+}
+
+
+// ************************************************************************* //
diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/diameterModels/IATE/IATEsources/IATEsource/IATEsource.H b/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/diameterModels/IATE/IATEsources/IATEsource/IATEsource.H
new file mode 100644
index 0000000000000000000000000000000000000000..f21f4f387679fadcc06523097ad84c6e8e938023
--- /dev/null
+++ b/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/diameterModels/IATE/IATEsources/IATEsource/IATEsource.H
@@ -0,0 +1,192 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::diameterModels::IATEsource
+
+Description
+    IATE (Interfacial Area Transport Equation) bubble diameter model
+    run-time selectable sources.
+
+SourceFiles
+    IATEsource.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef IATEsource_H
+#define IATEsource_H
+
+#include "IATE.H"
+#include "mathematicalConstants.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace diameterModels
+{
+
+/*---------------------------------------------------------------------------*\
+                           Class IATEsource Declaration
+\*---------------------------------------------------------------------------*/
+
+class IATEsource
+{
+
+protected:
+
+    // Protected data
+
+        //- Reference to the IATE this source applies to
+        const IATE& iate_;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("IATEsource");
+
+
+    // Declare run-time constructor selection table
+
+        declareRunTimeSelectionTable
+        (
+            autoPtr,
+            IATEsource,
+            dictionary,
+            (
+                const IATE& iate,
+                const dictionary& dict
+            ),
+            (iate, dict)
+        );
+
+
+    //- Class used for the read-construction of
+    //  PtrLists of IATE sources
+    class iNew
+    {
+        const IATE& iate_;
+
+    public:
+
+        iNew(const IATE& iate)
+        :
+            iate_(iate)
+        {}
+
+        autoPtr<IATEsource> operator()(Istream& is) const
+        {
+            word type(is);
+            dictionary dict(is);
+            return IATEsource::New(type, iate_, dict);
+        }
+    };
+
+
+    // Constructors
+
+        IATEsource(const IATE& iate)
+        :
+            iate_(iate)
+        {}
+
+        autoPtr<IATEsource> clone() const
+        {
+            notImplemented("autoPtr<IATEsource> clone() const");
+            return autoPtr<IATEsource>(NULL);
+        }
+
+
+    // Selectors
+
+        static autoPtr<IATEsource> New
+        (
+            const word& type,
+            const IATE& iate,
+            const dictionary& dict
+        );
+
+
+    //- Destructor
+    virtual ~IATEsource()
+    {}
+
+
+    // Member Functions
+
+        const phaseModel& phase() const
+        {
+            return iate_.phase();
+        }
+
+        const twoPhaseSystem& fluid() const
+        {
+            return iate_.phase().fluid();
+        }
+
+        const phaseModel& otherPhase() const
+        {
+            return phase().otherPhase();
+        }
+
+        scalar phi() const
+        {
+            return 1.0/(36*constant::mathematical::pi);
+        }
+
+        //- Return the bubble relative velocity
+        tmp<volScalarField> Ur() const;
+
+        //- Return the bubble turbulent velocity
+        tmp<volScalarField> Ut() const;
+
+        //- Return the bubble Reynolds number
+        tmp<volScalarField> Re() const;
+
+        //- Return the bubble drag coefficient
+        tmp<volScalarField> CD() const;
+
+        //- Return the bubble Morton number
+        tmp<volScalarField> Mo() const;
+
+        //- Return the bubble Eotvos number
+        tmp<volScalarField> Eo() const;
+
+        //- Return the bubble Webber number
+        tmp<volScalarField> We() const;
+
+        virtual tmp<volScalarField> R() const = 0;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace diameterModels
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/diameterModels/IATE/IATEsources/dummy/dummy.C b/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/diameterModels/IATE/IATEsources/dummy/dummy.C
new file mode 100644
index 0000000000000000000000000000000000000000..d442373b0c3193311f0e48855ff056a6f3c3be0b
--- /dev/null
+++ b/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/diameterModels/IATE/IATEsources/dummy/dummy.C
@@ -0,0 +1,66 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "dummy.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace diameterModels
+{
+namespace IATEsources
+{
+    defineTypeNameAndDebug(dummy, 0);
+    addToRunTimeSelectionTable(IATEsource, dummy, word);
+}
+}
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+Foam::tmp<Foam::volScalarField>
+Foam::diameterModels::IATEsources::dummy::R() const
+{
+    return tmp<volScalarField>
+    (
+        new volScalarField
+        (
+            IOobject
+            (
+                "R",
+                iate_.phase().U().time().timeName(),
+                iate_.phase().mesh()
+            ),
+            iate_.phase().U().mesh(),
+            dimensionedScalar("R", dimless/dimTime, 0)
+        )
+    );
+}
+
+
+// ************************************************************************* //
diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/diameterModels/IATE/IATEsources/dummy/dummy.H b/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/diameterModels/IATE/IATEsources/dummy/dummy.H
new file mode 100644
index 0000000000000000000000000000000000000000..9fcdbf53252527c7287fd27cecc7bf0885940323
--- /dev/null
+++ b/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/diameterModels/IATE/IATEsources/dummy/dummy.H
@@ -0,0 +1,96 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::diameterModels::IATEsources::dummy
+
+Description
+
+SourceFiles
+    dummy.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef dummy_H
+#define dummy_H
+
+#include "IATEsource.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace diameterModels
+{
+namespace IATEsources
+{
+
+/*---------------------------------------------------------------------------*\
+                           Class dummy Declaration
+\*---------------------------------------------------------------------------*/
+
+class dummy
+:
+    public IATEsource
+{
+
+public:
+
+    //- Runtime type information
+    TypeName("dummy");
+
+    // Constructors
+
+        dummy
+        (
+            const word& name,
+            const IATE& iate,
+            const dictionary& dict
+        )
+        :
+            IATEsource(iate)
+        {}
+
+
+    //- Destructor
+    virtual ~dummy()
+    {}
+
+
+    // Member Functions
+
+        virtual tmp<volScalarField> R() const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace IATEsources
+} // End namespace diameterModels
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/diameterModels/IATE/IATEsources/randomCoalescence/randomCoalescence.C b/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/diameterModels/IATE/IATEsources/randomCoalescence/randomCoalescence.C
new file mode 100644
index 0000000000000000000000000000000000000000..3ddfd70713ed29090cda85863d3bedcfe4f3bd92
--- /dev/null
+++ b/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/diameterModels/IATE/IATEsources/randomCoalescence/randomCoalescence.C
@@ -0,0 +1,109 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "randomCoalescence.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace diameterModels
+{
+namespace IATEsources
+{
+    defineTypeNameAndDebug(randomCoalescence, 0);
+    addToRunTimeSelectionTable(IATEsource, randomCoalescence, dictionary);
+}
+}
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::diameterModels::IATEsources::randomCoalescence::
+randomCoalescence
+(
+    const IATE& iate,
+    const dictionary& dict
+)
+:
+    IATEsource(iate),
+    Crc_("Crc", dimless, dict.lookup("Crc")),
+    C_("C", dimless, dict.lookup("C")),
+    alphaMax_("alphaMax", dimless, dict.lookup("alphaMax"))
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+Foam::tmp<Foam::volScalarField>
+Foam::diameterModels::IATEsources::randomCoalescence::R() const
+{
+    tmp<volScalarField> tR
+    (
+        new volScalarField
+        (
+            IOobject
+            (
+                "R",
+                iate_.phase().U().time().timeName(),
+                iate_.phase().mesh()
+            ),
+            iate_.phase().U().mesh(),
+            dimensionedScalar("R", dimless/dimTime, 0)
+        )
+    );
+
+    volScalarField R = tR();
+
+    scalar Crc = Crc_.value();
+    scalar C = C_.value();
+    scalar alphaMax = alphaMax_.value();
+    volScalarField Ut(this->Ut());
+    const volScalarField& alpha = phase();
+    const volScalarField& kappai = iate_.kappai();
+    scalar cbrtAlphaMax = cbrt(alphaMax);
+
+    forAll(R, celli)
+    {
+        if (alpha[celli] < alphaMax - SMALL)
+        {
+            scalar cbrtAlphaMaxMAlpha = cbrtAlphaMax - cbrt(alpha[celli]);
+
+            R[celli] =
+                12*phi()*kappai[celli]*alpha[celli]
+               *Crc
+               *Ut[celli]
+               *(1 - exp(-C*cbrt(alpha[celli]*alphaMax)/cbrtAlphaMaxMAlpha))
+               /(cbrtAlphaMax*cbrtAlphaMaxMAlpha);
+        }
+    }
+
+    return tR;
+}
+
+
+// ************************************************************************* //
diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/diameterModels/IATE/IATEsources/randomCoalescence/randomCoalescence.H b/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/diameterModels/IATE/IATEsources/randomCoalescence/randomCoalescence.H
new file mode 100644
index 0000000000000000000000000000000000000000..d8c03b47c3845eb663caad12839c3451e64fdab7
--- /dev/null
+++ b/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/diameterModels/IATE/IATEsources/randomCoalescence/randomCoalescence.H
@@ -0,0 +1,108 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::diameterModels::IATEsources::randomCoalescence
+
+Description
+    Random coalescence IATE source as defined in paper:
+
+    \verbatim
+        "Development of Interfacial Area Transport Equation"
+        M. Ishii,
+        S. Kim,
+        J Kelly,
+        Nuclear Engineering and Technology, Vol.37 No.6 December 2005
+    \endverbatim
+
+
+SourceFiles
+    randomCoalescence.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef randomCoalescence_H
+#define randomCoalescence_H
+
+#include "IATEsource.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace diameterModels
+{
+namespace IATEsources
+{
+
+/*---------------------------------------------------------------------------*\
+                           Class randomCoalescence Declaration
+\*---------------------------------------------------------------------------*/
+
+class randomCoalescence
+:
+    public IATEsource
+{
+    // Private data
+
+        dimensionedScalar Crc_;
+        dimensionedScalar C_;
+        dimensionedScalar alphaMax_;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("randomCoalescence");
+
+    // Constructors
+
+        randomCoalescence
+        (
+            const IATE& iate,
+            const dictionary& dict
+        );
+
+
+    //- Destructor
+    virtual ~randomCoalescence()
+    {}
+
+
+    // Member Functions
+
+        virtual tmp<volScalarField> R() const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace IATEsources
+} // End namespace diameterModels
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/diameterModels/IATE/IATEsources/turbulentBreakUp/turbulentBreakUp.C b/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/diameterModels/IATE/IATEsources/turbulentBreakUp/turbulentBreakUp.C
new file mode 100644
index 0000000000000000000000000000000000000000..39bd0daf89b8a10c212a31af11075435b96515e0
--- /dev/null
+++ b/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/diameterModels/IATE/IATEsources/turbulentBreakUp/turbulentBreakUp.C
@@ -0,0 +1,104 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "turbulentBreakUp.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace diameterModels
+{
+namespace IATEsources
+{
+    defineTypeNameAndDebug(turbulentBreakUp, 0);
+    addToRunTimeSelectionTable(IATEsource, turbulentBreakUp, dictionary);
+}
+}
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::diameterModels::IATEsources::turbulentBreakUp::
+turbulentBreakUp
+(
+    const IATE& iate,
+    const dictionary& dict
+)
+:
+    IATEsource(iate),
+    Cti_("Cti", dimless, dict.lookup("Cti")),
+    WeCr_("WeCr", dimless, dict.lookup("WeCr"))
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+Foam::tmp<Foam::volScalarField>
+Foam::diameterModels::IATEsources::turbulentBreakUp::R() const
+{
+    tmp<volScalarField> tR
+    (
+        new volScalarField
+        (
+            IOobject
+            (
+                "R",
+                iate_.phase().U().time().timeName(),
+                iate_.phase().mesh()
+            ),
+            iate_.phase().U().mesh(),
+            dimensionedScalar("R", dimless/dimTime, 0)
+        )
+    );
+
+    volScalarField R = tR();
+
+    scalar Cti = Cti_.value();
+    scalar WeCr = WeCr_.value();
+    volScalarField Ut(this->Ut());
+    volScalarField We(this->We());
+    const volScalarField& d(iate_.d()());
+
+    forAll(R, celli)
+    {
+        if (We[celli] > WeCr)
+        {
+            R[celli] =
+                (1.0/3.0)
+               *Cti/d[celli]
+               *Ut[celli]
+               *sqrt(1 - WeCr/We[celli])
+               *exp(-WeCr/We[celli]);
+        }
+    }
+
+    return tR;
+}
+
+
+// ************************************************************************* //
diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/diameterModels/IATE/IATEsources/turbulentBreakUp/turbulentBreakUp.H b/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/diameterModels/IATE/IATEsources/turbulentBreakUp/turbulentBreakUp.H
new file mode 100644
index 0000000000000000000000000000000000000000..c900b5d592d07a88eb551b33f8e803296916210c
--- /dev/null
+++ b/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/diameterModels/IATE/IATEsources/turbulentBreakUp/turbulentBreakUp.H
@@ -0,0 +1,106 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::diameterModels::IATEsources::turbulentBreakUp
+
+Description
+    Turbulence-induced break-up IATE source as defined in paper:
+
+    \verbatim
+        "Development of Interfacial Area Transport Equation"
+        M. Ishii,
+        S. Kim,
+        J Kelly,
+        Nuclear Engineering and Technology, Vol.37 No.6 December 2005
+    \endverbatim
+
+SourceFiles
+    turbulentBreakUp.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef turbulentBreakUp_H
+#define turbulentBreakUp_H
+
+#include "IATEsource.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace diameterModels
+{
+namespace IATEsources
+{
+
+/*---------------------------------------------------------------------------*\
+                           Class turbulentBreakUp Declaration
+\*---------------------------------------------------------------------------*/
+
+class turbulentBreakUp
+:
+    public IATEsource
+{
+    // Private data
+
+        dimensionedScalar Cti_;
+        dimensionedScalar WeCr_;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("turbulentBreakUp");
+
+    // Constructors
+
+        turbulentBreakUp
+        (
+            const IATE& iate,
+            const dictionary& dict
+        );
+
+
+    //- Destructor
+    virtual ~turbulentBreakUp()
+    {}
+
+
+    // Member Functions
+
+        virtual tmp<volScalarField> R() const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace IATEsources
+} // End namespace diameterModels
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/diameterModels/IATE/IATEsources/wakeEntrainmentCoalescence/wakeEntrainmentCoalescence.C b/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/diameterModels/IATE/IATEsources/wakeEntrainmentCoalescence/wakeEntrainmentCoalescence.C
new file mode 100644
index 0000000000000000000000000000000000000000..66c324c761f2c99af4db2afef7173c29c252d584
--- /dev/null
+++ b/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/diameterModels/IATE/IATEsources/wakeEntrainmentCoalescence/wakeEntrainmentCoalescence.C
@@ -0,0 +1,72 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "wakeEntrainmentCoalescence.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace diameterModels
+{
+namespace IATEsources
+{
+    defineTypeNameAndDebug(wakeEntrainmentCoalescence, 0);
+    addToRunTimeSelectionTable
+    (
+        IATEsource,
+        wakeEntrainmentCoalescence,
+        dictionary
+    );
+}
+}
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::diameterModels::IATEsources::wakeEntrainmentCoalescence::
+wakeEntrainmentCoalescence
+(
+    const IATE& iate,
+    const dictionary& dict
+)
+:
+    IATEsource(iate),
+    Cwe_("Cwe", dimless, dict.lookup("Cwe"))
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+Foam::tmp<Foam::volScalarField>
+Foam::diameterModels::IATEsources::wakeEntrainmentCoalescence::R() const
+{
+    return (-12)*phi()*Cwe_*cbrt(CD())*iate_.a()*Ur();
+}
+
+
+// ************************************************************************* //
diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/diameterModels/IATE/IATEsources/wakeEntrainmentCoalescence/wakeEntrainmentCoalescence.H b/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/diameterModels/IATE/IATEsources/wakeEntrainmentCoalescence/wakeEntrainmentCoalescence.H
new file mode 100644
index 0000000000000000000000000000000000000000..7c2df3991211cab13745d12b75661df5c04611c9
--- /dev/null
+++ b/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/diameterModels/IATE/IATEsources/wakeEntrainmentCoalescence/wakeEntrainmentCoalescence.H
@@ -0,0 +1,105 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::diameterModels::IATEsources::wakeEntrainmentCoalescence
+
+Description
+    Bubble coalescence due to wake entrainment IATE source as defined in paper:
+
+    \verbatim
+        "Development of Interfacial Area Transport Equation"
+        M. Ishii,
+        S. Kim,
+        J Kelly,
+        Nuclear Engineering and Technology, Vol.37 No.6 December 2005
+    \endverbatim
+
+SourceFiles
+    wakeEntrainmentCoalescence.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef wakeEntrainmentCoalescence_H
+#define wakeEntrainmentCoalescence_H
+
+#include "IATEsource.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace diameterModels
+{
+namespace IATEsources
+{
+
+/*---------------------------------------------------------------------------*\
+                           Class wakeEntrainmentCoalescence Declaration
+\*---------------------------------------------------------------------------*/
+
+class wakeEntrainmentCoalescence
+:
+    public IATEsource
+{
+    // Private data
+
+        dimensionedScalar Cwe_;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("wakeEntrainmentCoalescence");
+
+    // Constructors
+
+        wakeEntrainmentCoalescence
+        (
+            const IATE& iate,
+            const dictionary& dict
+        );
+
+
+    //- Destructor
+    virtual ~wakeEntrainmentCoalescence()
+    {}
+
+
+    // Member Functions
+
+        virtual tmp<volScalarField> R() const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace IATEsources
+} // End namespace diameterModels
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/diameterModels/constantDiameter/constantDiameter.C b/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/diameterModels/constantDiameter/constantDiameter.C
index 315f16089399f0b3f9d7a65f79df0549464e5d5a..233b9bf49ee90668f05df20057983fd37f4ab8ac 100644
--- a/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/diameterModels/constantDiameter/constantDiameter.C
+++ b/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/diameterModels/constantDiameter/constantDiameter.C
@@ -48,12 +48,12 @@ namespace diameterModels
 
 Foam::diameterModels::constant::constant
 (
-    const dictionary& dict,
+    const dictionary& diameterProperties,
     const phaseModel& phase
 )
 :
-    diameterModel(dict, phase),
-    d_("d", dimLength, dict.lookup("d"))
+    diameterModel(diameterProperties, phase),
+    d_("d", dimLength, diameterProperties_.lookup("d"))
 {}
 
 
@@ -84,4 +84,14 @@ Foam::tmp<Foam::volScalarField> Foam::diameterModels::constant::d() const
 }
 
 
+bool Foam::diameterModels::constant::read(const dictionary& phaseProperties)
+{
+    diameterModel::read(phaseProperties);
+
+    diameterProperties_.lookup("d") >> d_;
+
+    return true;
+}
+
+
 // ************************************************************************* //
diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/diameterModels/constantDiameter/constantDiameter.H b/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/diameterModels/constantDiameter/constantDiameter.H
index 071820c319f0202ca792e0531bd909aef52d8749..6f616c6f73df37fec0c840c918aa07f7ee1a2dd7 100644
--- a/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/diameterModels/constantDiameter/constantDiameter.H
+++ b/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/diameterModels/constantDiameter/constantDiameter.H
@@ -22,7 +22,7 @@ License
     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
 
 Class
-    Foam::constant
+    Foam::diameterModels::constant
 
 Description
     Constant dispersed-phase particle diameter model.
@@ -69,7 +69,7 @@ public:
         //- Construct from components
         constant
         (
-            const dictionary& dict,
+            const dictionary& diameterProperties,
             const phaseModel& phase
         );
 
@@ -80,7 +80,11 @@ public:
 
     // Member Functions
 
-        tmp<volScalarField> d() const;
+        //- Return the diameter as a field
+        virtual tmp<volScalarField> d() const;
+
+        //- Read diameterProperties dictionary
+        virtual bool read(const dictionary& diameterProperties);
 };
 
 
diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/diameterModels/diameterModel/diameterModel.C b/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/diameterModels/diameterModel/diameterModel.C
index 2618815d24d3e553a9a6e47f8452e352c3a7d6be..55225147acc6aa804e7aa87659621ac1bc85187d 100644
--- a/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/diameterModels/diameterModel/diameterModel.C
+++ b/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/diameterModels/diameterModel/diameterModel.C
@@ -38,11 +38,11 @@ namespace Foam
 
 Foam::diameterModel::diameterModel
 (
-    const dictionary& dict,
+    const dictionary& diameterProperties,
     const phaseModel& phase
 )
 :
-    dict_(dict),
+    diameterProperties_(diameterProperties),
     phase_(phase)
 {}
 
@@ -53,4 +53,18 @@ Foam::diameterModel::~diameterModel()
 {}
 
 
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void Foam::diameterModel::correct()
+{}
+
+
+bool Foam::diameterModel::read(const dictionary& phaseProperties)
+{
+    diameterProperties_ = phaseProperties.subDict(type() + "Coeffs");
+
+    return true;
+}
+
+
 // ************************************************************************* //
diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/diameterModels/diameterModel/diameterModel.H b/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/diameterModels/diameterModel/diameterModel.H
index c14f65e3645b255a235b3ecca78e1f44079ec6ed..b40a6cb8ca15fba63bd78a4fa9e2fb27550275d0 100644
--- a/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/diameterModels/diameterModel/diameterModel.H
+++ b/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/diameterModels/diameterModel/diameterModel.H
@@ -36,26 +36,27 @@ SourceFiles
 #ifndef diameterModel_H
 #define diameterModel_H
 
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
 #include "dictionary.H"
 #include "phaseModel.H"
 #include "runTimeSelectionTables.H"
 
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
 namespace Foam
 {
 
 /*---------------------------------------------------------------------------*\
-                           Class diameterModel Declaration
+                       Class diameterModel Declaration
 \*---------------------------------------------------------------------------*/
 
 class diameterModel
 {
+
 protected:
 
     // Protected data
 
-        const dictionary& dict_;
+        dictionary diameterProperties_;
         const phaseModel& phase_;
 
 
@@ -73,10 +74,10 @@ public:
             diameterModel,
             dictionary,
             (
-                const dictionary& dict,
+                const dictionary& diameterProperties,
                 const phaseModel& phase
             ),
-            (dict, phase)
+            (diameterProperties, phase)
         );
 
 
@@ -84,7 +85,7 @@ public:
 
         diameterModel
         (
-            const dictionary& dict,
+            const dictionary& diameterProperties,
             const phaseModel& phase
         );
 
@@ -97,15 +98,33 @@ public:
 
         static autoPtr<diameterModel> New
         (
-            const dictionary& dict,
+            const dictionary& diameterProperties,
             const phaseModel& phase
         );
 
 
     // Member Functions
 
+        //- Return the phase diameter properties dictionary
+        const dictionary& diameterProperties() const
+        {
+            return diameterProperties_;
+        }
+
+        //- Return the phase
+        const phaseModel& phase() const
+        {
+            return phase_;
+        }
+
         //- Return the phase mean diameter field
         virtual tmp<volScalarField> d() const = 0;
+
+        //- Correct the diameter field
+        virtual void correct();
+
+        //- Read phaseProperties dictionary
+        virtual bool read(const dictionary& phaseProperties) = 0;
 };
 
 
diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/diameterModels/isothermalDiameter/isothermalDiameter.C b/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/diameterModels/isothermalDiameter/isothermalDiameter.C
index 4b9e57150390d1412a9e284549923b657e5be686..aab9b469c7c36d054ca17a65d41aeb5277c0b33f 100644
--- a/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/diameterModels/isothermalDiameter/isothermalDiameter.C
+++ b/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/diameterModels/isothermalDiameter/isothermalDiameter.C
@@ -48,13 +48,13 @@ namespace diameterModels
 
 Foam::diameterModels::isothermal::isothermal
 (
-    const dictionary& dict,
+    const dictionary& diameterProperties,
     const phaseModel& phase
 )
 :
-    diameterModel(dict, phase),
-    d0_("d0", dimLength, dict.lookup("d0")),
-    p0_("p0", dimPressure, dict.lookup("p0"))
+    diameterModel(diameterProperties, phase),
+    d0_("d0", dimLength, diameterProperties_.lookup("d0")),
+    p0_("p0", dimPressure, diameterProperties_.lookup("p0"))
 {}
 
 
@@ -77,4 +77,15 @@ Foam::tmp<Foam::volScalarField> Foam::diameterModels::isothermal::d() const
 }
 
 
+bool Foam::diameterModels::isothermal::read(const dictionary& phaseProperties)
+{
+    diameterModel::read(phaseProperties);
+
+    diameterProperties_.lookup("d0") >> d0_;
+    diameterProperties_.lookup("p0") >> p0_;
+
+    return true;
+}
+
+
 // ************************************************************************* //
diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/diameterModels/isothermalDiameter/isothermalDiameter.H b/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/diameterModels/isothermalDiameter/isothermalDiameter.H
index 0d155ad16e0bb4a6aca3ab251e6c89c68c89cfbd..6e997a210aec2ffbc98797a8c0a94b43b1b8bfb6 100644
--- a/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/diameterModels/isothermalDiameter/isothermalDiameter.H
+++ b/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/diameterModels/isothermalDiameter/isothermalDiameter.H
@@ -22,7 +22,7 @@ License
     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
 
 Class
-    Foam::isothermal
+    Foam::diameterModels::isothermal
 
 Description
     Isothermal dispersed-phase particle diameter model.
@@ -72,7 +72,7 @@ public:
         //- Construct from components
         isothermal
         (
-            const dictionary& dict,
+            const dictionary& diameterProperties,
             const phaseModel& phase
         );
 
@@ -83,7 +83,11 @@ public:
 
     // Member Functions
 
-        tmp<volScalarField> d() const;
+        //- Return the diameter field
+        virtual tmp<volScalarField> d() const;
+
+        //- Read phaseProperties dictionary
+        virtual bool read(const dictionary& phaseProperties);
 };
 
 
diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/phaseModel/phaseModel.C b/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/phaseModel/phaseModel.C
index 697371c04440e6740780fb4d6a0edf7b7bc33ed3..da43a49881906051ca5926a17447e2037547b2ec 100644
--- a/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/phaseModel/phaseModel.C
+++ b/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/phaseModel/phaseModel.C
@@ -26,6 +26,8 @@ License
 #include "phaseModel.H"
 #include "twoPhaseSystem.H"
 #include "diameterModel.H"
+#include "fvMatrix.H"
+#include "PhaseIncompressibleTurbulenceModel.H"
 #include "dragModel.H"
 #include "heatTransferModel.H"
 #include "fixedValueFvPatchFields.H"
@@ -73,6 +75,17 @@ Foam::phaseModel::phaseModel
             IOobject::AUTO_WRITE
         ),
         fluid.mesh()
+    ),
+    phiAlpha_
+    (
+        IOobject
+        (
+            IOobject::groupName("alphaPhi", name_),
+            fluid.mesh().time().timeName(),
+            fluid.mesh()
+        ),
+        fluid.mesh(),
+        dimensionedScalar("0", dimensionSet(0, 3, -1, 0, 0), 0)
     )
 {
     thermo_->validate("phaseModel " + name_, "h", "e");
@@ -153,6 +166,16 @@ Foam::phaseModel::phaseModel
         phaseDict_,
         *this
     );
+
+    turbulence_ =
+        PhaseIncompressibleTurbulenceModel<phaseModel>::New
+        (
+            *this,
+            U_,
+            phiAlpha_,
+            phi(),
+            *this
+        );
 }
 
 
@@ -164,10 +187,39 @@ Foam::phaseModel::~phaseModel()
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
+const Foam::phaseModel& Foam::phaseModel::otherPhase() const
+{
+    return fluid_.otherPhase(*this);
+}
+
+
 Foam::tmp<Foam::volScalarField> Foam::phaseModel::d() const
 {
     return dPtr_().d();
 }
 
+Foam::PhaseIncompressibleTurbulenceModel<Foam::phaseModel>&
+Foam::phaseModel::turbulence()
+{
+    return turbulence_();
+}
+
+const Foam::PhaseIncompressibleTurbulenceModel<Foam::phaseModel>&
+Foam::phaseModel::turbulence() const
+{
+    return turbulence_();
+}
+
+void Foam::phaseModel::correct()
+{
+    return dPtr_->correct();
+}
+
+bool Foam::phaseModel::read(const dictionary& phaseProperties)
+{
+    phaseDict_ = phaseProperties.subDict(name_);
+    return dPtr_->read(phaseDict_);
+}
+
 
 // ************************************************************************* //
diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/phaseModel/phaseModel.H b/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/phaseModel/phaseModel.H
index ed2a6115462185851dc9f52857231cd65dfa8c46..7ca425520770823c0b6c1fd4cfdbf94e12423a4d 100644
--- a/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/phaseModel/phaseModel.H
+++ b/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/phaseModel/phaseModel.H
@@ -48,6 +48,9 @@ namespace Foam
 class twoPhaseSystem;
 class diameterModel;
 
+template<class Phase>
+class PhaseIncompressibleTurbulenceModel;
+
 
 /*---------------------------------------------------------------------------*\
                            Class phaseModel Declaration
@@ -74,12 +77,18 @@ class phaseModel
         //- Velocity
         volVectorField U_;
 
-        //- Fluxes
+        //- Volumetric flux of the phase
+        surfaceScalarField phiAlpha_;
+
+        //- Volumetric flux of the phase
         autoPtr<surfaceScalarField> phiPtr_;
 
         //- Diameter model
         autoPtr<diameterModel> dPtr_;
 
+        //- turbulence model
+        autoPtr<PhaseIncompressibleTurbulenceModel<phaseModel> > turbulence_;
+
 
 public:
 
@@ -99,19 +108,47 @@ public:
 
     // Member Functions
 
+        //- Return the name of this phase
+        const word& name() const
+        {
+            return name_;
+        }
+
         //- Return the twoPhaseSystem to which this phase belongs
         const twoPhaseSystem& fluid() const
         {
             return fluid_;
         }
 
-        const word& name() const
+        //- Return the other phase in this two-phase system
+        const phaseModel& otherPhase() const;
+
+        //- Return the Sauter-mean diameter
+        tmp<volScalarField> d() const;
+
+        //- Return the turbulence model
+        const PhaseIncompressibleTurbulenceModel<phaseModel>&
+            turbulence() const;
+
+        //- Return non-const access to the turbulence model
+        //  for correction
+        PhaseIncompressibleTurbulenceModel<phaseModel>&
+            turbulence();
+
+        //- Return the thermophysical model
+        const rhoThermo& thermo() const
         {
-            return name_;
+            return thermo_();
         }
 
-        tmp<volScalarField> d() const;
+        //- Return non-const access to the thermophysical model
+        //  for correction
+        rhoThermo& thermo()
+        {
+            return thermo_();
+        }
 
+        //- Return the laminar viscosity
         tmp<volScalarField> nu() const
         {
             return thermo_->nu();
@@ -123,57 +160,71 @@ public:
             return thermo_->nu(patchi);
         }
 
+        //- Return the thermal conductivity
         tmp<volScalarField> kappa() const
         {
             return thermo_->kappa();
         }
 
+        //- Return the specific heat capacity
         tmp<volScalarField> Cp() const
         {
             return thermo_->Cp();
         }
 
+        //- Return the density
         const volScalarField& rho() const
         {
             return thermo_->rho();
         }
 
-        const rhoThermo& thermo() const
-        {
-            return thermo_();
-        }
-
-        rhoThermo& thermo()
-        {
-            return thermo_();
-        }
-
+        //- Return the velocity
         const volVectorField& U() const
         {
             return U_;
         }
 
+        //- Return non-const access to the velocity
+        //  Used in the momentum equation
         volVectorField& U()
         {
             return U_;
         }
 
+        //- Return the volumetric flux
         const surfaceScalarField& phi() const
         {
             return phiPtr_();
         }
 
+        //- Return non-const access to the volumetric flux
         surfaceScalarField& phi()
         {
             return phiPtr_();
         }
 
-        //- Dummy correct
-        void correct()
-        {}
+        //- Return the volumetric flux of the phase
+        const surfaceScalarField& phiAlpha() const
+        {
+            return phiAlpha_;
+        }
+
+        //- Return non-const access to the volumetric flux of the phase
+        surfaceScalarField& phiAlpha()
+        {
+            return phiAlpha_;
+        }
+
+        //- Correct the phase properties
+        //  other than the thermodynamics and turbulence
+        //  which have special treatment
+        void correct();
+
+        //- Read phaseProperties dictionary
+        virtual bool read(const dictionary& phaseProperties);
 
-        //- Dummy read
-        bool read()
+        //- Dummy Read for transportModel
+        virtual bool read()
         {
             return true;
         }
diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/twoPhaseSystem.C b/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/twoPhaseSystem.C
index 83902a7b7d7385b1c5619d769e11287a05ca8c0f..8ba3382b86967b3d0d74522d8b0bcff64b8c215e 100644
--- a/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/twoPhaseSystem.C
+++ b/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/twoPhaseSystem.C
@@ -24,9 +24,19 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "twoPhaseSystem.H"
+#include "fvMatrix.H"
+#include "PhaseIncompressibleTurbulenceModel.H"
 #include "surfaceInterpolate.H"
-#include "fixedValueFvsPatchFields.H"
+#include "MULES.H"
+#include "subCycle.H"
+#include "fvcDdt.H"
+#include "fvcDiv.H"
+#include "fvcSnGrad.H"
+#include "fvcFlux.H"
 #include "fvcCurl.H"
+#include "fvmDdt.H"
+#include "fvmLaplacian.H"
+#include "fixedValueFvsPatchFields.H"
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
@@ -63,6 +73,37 @@ Foam::twoPhaseSystem::twoPhaseSystem
         wordList(lookup("phases"))[1]
     ),
 
+    phi_
+    (
+        IOobject
+        (
+            "phi",
+            mesh.time().timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        this->calcPhi()
+    ),
+
+    dgdt_
+    (
+        IOobject
+        (
+            "dgdt",
+            mesh.time().timeName(),
+            mesh
+        ),
+        pos(phase2_)*fvc::div(phi_)/max(phase2_, scalar(0.0001))
+    ),
+
+    sigma_
+    (
+        "sigma",
+        dimensionSet(1, 0, -2, 0, 0),
+        lookup("sigma")
+    ),
+
     Cvm_
     (
         "Cvm",
@@ -170,7 +211,7 @@ Foam::tmp<Foam::volVectorField> Foam::twoPhaseSystem::U() const
 }
 
 
-Foam::tmp<Foam::surfaceScalarField> Foam::twoPhaseSystem::phi() const
+Foam::tmp<Foam::surfaceScalarField> Foam::twoPhaseSystem::calcPhi() const
 {
     return
         fvc::interpolate(phase1_)*phase1_.phi()
@@ -366,18 +407,242 @@ Foam::tmp<Foam::volScalarField> Foam::twoPhaseSystem::heatTransferCoeff() const
 }
 
 
+void Foam::twoPhaseSystem::solve()
+{
+    const Time& runTime = mesh_.time();
+
+    volScalarField& alpha1 = phase1_;
+    volScalarField& alpha2 = phase2_;
+
+    const surfaceScalarField& phi1 = phase1_.phi();
+    const surfaceScalarField& phi2 = phase2_.phi();
+
+    const dictionary& alphaControls = mesh_.solverDict
+    (
+        alpha1.name()
+    );
+
+    label nAlphaSubCycles(readLabel(alphaControls.lookup("nAlphaSubCycles")));
+    label nAlphaCorr(readLabel(alphaControls.lookup("nAlphaCorr")));
+    Switch implicitPhasePressure
+    (
+        alphaControls.lookupOrDefault<Switch>("implicitPhasePressure", false)
+    );
+
+    word alphaScheme("div(phi," + alpha1.name() + ')');
+    word alpharScheme("div(phir," + alpha1.name() + ')');
+
+    alpha1.correctBoundaryConditions();
+
+
+    surfaceScalarField phic("phic", phi_);
+    surfaceScalarField phir("phir", phi1 - phi2);
+
+    surfaceScalarField alpha1f(fvc::interpolate(max(alpha1, scalar(0))));
+
+    tmp<surfaceScalarField> pPrimeByA;
+
+    if (implicitPhasePressure)
+    {
+        const volScalarField& rAU1 = mesh_.lookupObject<volScalarField>
+        (
+            IOobject::groupName("rAU", phase1_.name())
+        );
+        const volScalarField& rAU2 = mesh_.lookupObject<volScalarField>
+        (
+            IOobject::groupName("rAU", phase2_.name())
+        );
+
+        pPrimeByA =
+            fvc::interpolate((1.0/phase1_.rho())
+           *rAU1*phase1_.turbulence().pPrime())
+          + fvc::interpolate((1.0/phase2_.rho())
+           *rAU2*phase2_.turbulence().pPrime());
+
+        surfaceScalarField phiP
+        (
+            pPrimeByA()*fvc::snGrad(alpha1, "bounded")*mesh_.magSf()
+        );
+
+        phic += alpha1f*phiP;
+        phir += phiP;
+    }
+
+    for (int acorr=0; acorr<nAlphaCorr; acorr++)
+    {
+        volScalarField::DimensionedInternalField Sp
+        (
+            IOobject
+            (
+                "Sp",
+                runTime.timeName(),
+                mesh_
+            ),
+            mesh_,
+            dimensionedScalar("Sp", dgdt_.dimensions(), 0.0)
+        );
+
+        volScalarField::DimensionedInternalField Su
+        (
+            IOobject
+            (
+                "Su",
+                runTime.timeName(),
+                mesh_
+            ),
+            // Divergence term is handled explicitly to be
+            // consistent with the explicit transport solution
+            fvc::div(phi_)*min(alpha1, scalar(1))
+        );
+
+        forAll(dgdt_, celli)
+        {
+            if (dgdt_[celli] > 0.0 && alpha1[celli] > 0.0)
+            {
+                Sp[celli] -= dgdt_[celli]*alpha1[celli];
+                Su[celli] += dgdt_[celli]*alpha1[celli];
+            }
+            else if (dgdt_[celli] < 0.0 && alpha1[celli] < 1.0)
+            {
+                Sp[celli] += dgdt_[celli]*(1.0 - alpha1[celli]);
+            }
+        }
+
+        dimensionedScalar totalDeltaT = runTime.deltaT();
+        if (nAlphaSubCycles > 1)
+        {
+            phase1_.phiAlpha() =
+                dimensionedScalar("0", phase1_.phiAlpha().dimensions(), 0);
+        }
+
+        for
+        (
+            subCycle<volScalarField> alphaSubCycle(alpha1, nAlphaSubCycles);
+            !(++alphaSubCycle).end();
+        )
+        {
+            surfaceScalarField alphaPhic1
+            (
+                fvc::flux
+                (
+                    phic,
+                    alpha1,
+                    alphaScheme
+                )
+              + fvc::flux
+                (
+                    -fvc::flux(-phir, scalar(1) - alpha1, alpharScheme),
+                    alpha1,
+                    alpharScheme
+                )
+            );
+
+            // Ensure that the flux at inflow BCs is preserved
+            forAll(alphaPhic1.boundaryField(), patchi)
+            {
+                fvsPatchScalarField& alphaPhic1p =
+                    alphaPhic1.boundaryField()[patchi];
+
+                if (!alphaPhic1p.coupled())
+                {
+                    const scalarField& phi1p = phi1.boundaryField()[patchi];
+                    const scalarField& alpha1p = alpha1.boundaryField()[patchi];
+
+                    forAll(alphaPhic1p, facei)
+                    {
+                        if (phi1p[facei] < 0)
+                        {
+                            alphaPhic1p[facei] = alpha1p[facei]*phi1p[facei];
+                        }
+                    }
+                }
+            }
+
+            MULES::explicitSolve
+            (
+                geometricOneField(),
+                alpha1,
+                phi_,
+                alphaPhic1,
+                Sp,
+                Su,
+                1,
+                0
+            );
+
+            if (nAlphaSubCycles > 1)
+            {
+                phase1_.phiAlpha() += (runTime.deltaT()/totalDeltaT)*alphaPhic1;
+            }
+            else
+            {
+                phase1_.phiAlpha() = alphaPhic1;
+            }
+        }
+
+        if (implicitPhasePressure)
+        {
+            fvScalarMatrix alpha1Eqn
+            (
+                fvm::ddt(alpha1) - fvc::ddt(alpha1)
+              - fvm::laplacian(alpha1f*pPrimeByA, alpha1, "bounded")
+            );
+
+            alpha1Eqn.relax();
+            alpha1Eqn.solve();
+
+            phase1_.phiAlpha() += alpha1Eqn.flux();
+        }
+
+        phase2_.phiAlpha() = phi_ - phase1_.phiAlpha();
+        alpha2 = scalar(1) - alpha1;
+
+        Info<< alpha1.name() << " volume fraction = "
+            << alpha1.weightedAverage(mesh_.V()).value()
+            << "  Min(alpha1) = " << min(alpha1).value()
+            << "  Max(alpha1) = " << max(alpha1).value()
+            << endl;
+    }
+}
+
+
+void Foam::twoPhaseSystem::correct()
+{
+    phase1_.correct();
+    phase2_.correct();
+}
+
+
+void Foam::twoPhaseSystem::correctTurbulence()
+{
+    phase1_.turbulence().correct();
+    phase2_.turbulence().correct();
+}
+
+
 bool Foam::twoPhaseSystem::read()
 {
     if (regIOobject::read())
     {
         bool readOK = true;
 
-        readOK &= phase1_.read();
-        readOK &= phase2_.read();
+        readOK &= phase1_.read(*this);
+        readOK &= phase2_.read(*this);
 
+        lookup("sigma") >> sigma_;
         lookup("Cvm") >> Cvm_;
         lookup("Cl") >> Cl_;
 
+        // drag1_->read(*this);
+        // drag2_->read(*this);
+
+        // heatTransfer1_->read(*this);
+        // heatTransfer2_->read(*this);
+
+        lookup("dispersedPhase") >> dispersedPhase_;
+        lookup("residualPhaseFraction") >> residualPhaseFraction_;
+        lookup("residualSlip") >> residualSlip_;
+
         return readOK;
     }
     else
diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/twoPhaseSystem.H b/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/twoPhaseSystem.H
index 3a7896ea02b13a5e41563b6ee6bc29b39b508687..0b501152c15f410169d8617da1c1cff494dce192 100644
--- a/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/twoPhaseSystem.H
+++ b/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/twoPhaseSystem.H
@@ -74,6 +74,12 @@ class twoPhaseSystem
         phaseModel phase1_;
         phaseModel phase2_;
 
+        surfaceScalarField phi_;
+
+        volScalarField dgdt_;
+
+        dimensionedScalar sigma_;
+
         dimensionedScalar Cvm_;
         dimensionedScalar Cl_;
 
@@ -88,6 +94,10 @@ class twoPhaseSystem
         dimensionedScalar residualSlip_;
 
 
+        //- Return the mixture flux
+        tmp<surfaceScalarField> calcPhi() const;
+
+
 public:
 
     // Constructors
@@ -140,6 +150,28 @@ public:
             return phase2_;
         }
 
+        //- Return the mixture flux
+        const surfaceScalarField& phi() const
+        {
+            return phi_;
+        }
+
+        //- Return the mixture flux
+        surfaceScalarField& phi()
+        {
+            return phi_;
+        }
+
+        const volScalarField& dgdt() const
+        {
+            return dgdt_;
+        }
+
+        volScalarField& dgdt()
+        {
+            return dgdt_;
+        }
+
         const dragModel& drag1() const
         {
             return drag1_();
@@ -193,8 +225,11 @@ public:
         //- Return the mixture velocity
         tmp<volVectorField> U() const;
 
-        //- Return the mixture flux
-        tmp<surfaceScalarField> phi() const;
+        //- Return the surface tension coefficient
+        dimensionedScalar sigma() const
+        {
+            return sigma_;
+        }
 
         //- Return the virtual-mass coefficient
         dimensionedScalar Cvm() const
@@ -208,9 +243,14 @@ public:
             return Cl_;
         }
 
-        //- Dummy correct
-        void correct()
-        {}
+        //- Solve for the two-phase-fractions
+        void solve();
+
+        //- Correct two-phase properties other than turbulence
+        void correct();
+
+        //- Correct two-phase turbulence
+        void correctTurbulence();
 
         //- Read base phaseProperties dictionary
         bool read();
diff --git a/qq b/qq
new file mode 100644
index 0000000000000000000000000000000000000000..131428ac6a86bc5dd12b16299e95b038445f4ae9
--- /dev/null
+++ b/qq
@@ -0,0 +1,102 @@
+OSspecific/POSIX/POSIX.C:pid_t Foam::pid()
+OSspecific/POSIX/POSIX.C:bool Foam::ping
+OSspecific/POSIX/POSIX.C:            "Foam::ping(const string&, ...)"
+OSspecific/POSIX/POSIX.C:            "Foam::ping(const string&, const label)"
+OSspecific/POSIX/POSIX.C:bool Foam::ping(const string& hostname, const label timeOut)
+OSspecific/POSIX/fileStat.H:    So e.g. Foam::ping first and hope nfs is running.
+OpenFOAM/global/constants/atomic/atomicConstants.C:                4.0*mathematical::pi
+OpenFOAM/global/constants/atomic/atomicConstants.C:                4.0*mathematical::pi
+OpenFOAM/global/constants/electromagnetic/electromagneticConstants.C:        4.0*mathematical::pi*1e-07
+OpenFOAM/global/constants/electromagnetic/electromagneticConstants.C:            1.0/(4.0*mathematical::pi)
+OpenFOAM/global/constants/physicoChemical/physicoChemicalConstants.C:            Foam::sqr(mathematical::pi)/60.0
+OpenFOAM/global/unitConversion/unitConversion.H:    return (deg*constant::mathematical::pi/180.0);
+OpenFOAM/global/unitConversion/unitConversion.H:    return (rad*180.0/constant::mathematical::pi);
+OpenFOAM/meshes/meshShapes/face/face.C:            angle = constant::mathematical::pi + edgeAngle;
+OpenFOAM/meshes/meshShapes/face/face.C:            angle = constant::mathematical::pi - edgeAngle;
+OpenFOAM/meshes/meshShapes/face/face.C:        scalar minDiff = constant::mathematical::pi;
+OpenFOAM/primitives/transform/transform.H:        return (3.0 + cos)*constant::mathematical::piByTwo;
+OpenFOAM/primitives/transform/transform.H:        return (1.0 - cos)*constant::mathematical::piByTwo;
+engine/engineTime/engineTime.C:Foam::scalar Foam::engineTime::pistonPosition(const scalar theta) const
+engine/engineTime/engineTime.C:Foam::dimensionedScalar Foam::engineTime::pistonPosition() const
+engine/engineTime/engineTime.C:Foam::dimensionedScalar Foam::engineTime::pistonDisplacement() const
+engine/engineTime/engineTime.C:Foam::dimensionedScalar Foam::engineTime::pistonSpeed() const
+engine/include/StCorr.H:                    Ak = sphereFraction*4.0*constant::mathematical::pi
+engine/include/StCorr.H:                           /(sphereFraction*4.0*constant::mathematical::pi),
+engine/include/StCorr.H:                    Ak = circleFraction*constant::mathematical::pi*thickness
+engine/include/StCorr.H:                              *constant::mathematical::pi
+finiteVolume/cfdTools/general/solutionControl/pimpleControl/pimpleControl.C:void Foam::pimpleControl::read()
+finiteVolume/cfdTools/general/solutionControl/pimpleControl/pimpleControl.C:bool Foam::pimpleControl::criteriaSatisfied()
+finiteVolume/cfdTools/general/solutionControl/pimpleControl/pimpleControl.C:Foam::pimpleControl::pimpleControl(fvMesh& mesh)
+finiteVolume/cfdTools/general/solutionControl/pimpleControl/pimpleControl.C:Foam::pimpleControl::~pimpleControl()
+finiteVolume/cfdTools/general/solutionControl/pimpleControl/pimpleControl.C:bool Foam::pimpleControl::loop()
+finiteVolume/cfdTools/general/solutionControl/pimpleControl/pimpleControl.H:    Foam::pimpleControl
+finiteVolume/cfdTools/general/solutionControl/pimpleControl/pimpleControlI.H:inline Foam::label Foam::pimpleControl::nCorrPIMPLE() const
+finiteVolume/cfdTools/general/solutionControl/pimpleControl/pimpleControlI.H:inline Foam::label Foam::pimpleControl::nCorrPISO() const
+finiteVolume/cfdTools/general/solutionControl/pimpleControl/pimpleControlI.H:inline Foam::label Foam::pimpleControl::corrPISO() const
+finiteVolume/cfdTools/general/solutionControl/pimpleControl/pimpleControlI.H:inline bool Foam::pimpleControl::correct()
+finiteVolume/cfdTools/general/solutionControl/pimpleControl/pimpleControlI.H:inline bool Foam::pimpleControl::storeInitialResiduals() const
+finiteVolume/cfdTools/general/solutionControl/pimpleControl/pimpleControlI.H:inline bool Foam::pimpleControl::finalIter() const
+finiteVolume/cfdTools/general/solutionControl/pimpleControl/pimpleControlI.H:inline bool Foam::pimpleControl::finalInnerIter() const
+finiteVolume/cfdTools/general/solutionControl/pimpleControl/pimpleControlI.H:inline bool Foam::pimpleControl::turbCorr() const
+finiteVolume/fields/fvPatchFields/derived/cylindricalInletVelocity/cylindricalInletVelocityFvPatchVectorField.C:        (rpm*constant::mathematical::pi/30.0)*(hatAxis) ^ d
+finiteVolume/fields/fvPatchFields/derived/swirlFlowRateInletVelocity/swirlFlowRateInletVelocityFvPatchVectorField.C:            (rpm*constant::mathematical::pi/30.0)
+fvOptions/sources/derived/rotorDiskSource/rotorDiskSource.C:    const scalar diameter = Foam::sqrt(4.0*sumArea/mathematical::pi);
+fvOptions/sources/derived/rotorDiskSource/rotorDiskSource.C:                alphaGeom = mathematical::pi - alphaGeom;
+fvOptions/sources/derived/rotorDiskSource/trimModel/targetCoeff/targetCoeffTrim.C:    scalar coeff1 = alpha_*sqr(rotor_.omega())*mathematical::pi;
+lagrangian/coalCombustion/submodels/surfaceReactionModel/COxidationDiffusionLimitedRate/COxidationDiffusionLimitedRate.C:    scalar dmC = 4.0*mathematical::pi*d*D_*YO2*Tc*rhoc/(Sb_*(T + Tc))*dt;
+lagrangian/coalCombustion/submodels/surfaceReactionModel/COxidationHurtMitchell/COxidationHurtMitchell.C:    const scalar Ap = constant::mathematical::pi*sqr(d);
+lagrangian/coalCombustion/submodels/surfaceReactionModel/COxidationKineticDiffusionLimitedRate/COxidationKineticDiffusionLimitedRate.C:    const scalar Ap = constant::mathematical::pi*sqr(d);
+lagrangian/coalCombustion/submodels/surfaceReactionModel/COxidationMurphyShaddix/COxidationMurphyShaddix.C:    const scalar Ap = constant::mathematical::pi*sqr(d);
+lagrangian/distributionModels/normal/normal.C:    scalar k = 2.0/(constant::mathematical::pi*a_) +  0.5*log(1.0 - y*y);
+lagrangian/dsmc/parcels/Templates/DsmcParcel/DsmcParcelI.H:    return constant::mathematical::pi*d_*d_;
+lagrangian/intermediate/submodels/CloudFunctionObjects/ParticleCollector/ParticleCollector.C:            scalar theta = pCyl[1] + constant::mathematical::pi;
+lagrangian/intermediate/submodels/CloudFunctionObjects/ParticleErosion/ParticleErosion.C:        const scalar alpha = mathematical::pi/2.0 - acos(nw & Udir);
+lagrangian/intermediate/submodels/Kinematic/CollisionModel/PairCollision/PairModel/PairSpringSliderDashpot/PairSpringSliderDashpot.H:                mathematical::pi/4.0
+lagrangian/intermediate/submodels/Kinematic/CollisionModel/PairCollision/WallModel/WallLocalSpringSliderDashpot/WallLocalSpringSliderDashpot.C:           *mathematical::pi*(sqr(pREff) - sqr(r_PW_mag))
+lagrangian/intermediate/submodels/Kinematic/CollisionModel/PairCollision/WallModel/WallSpringSliderDashpot/WallSpringSliderDashpot.C:           *mathematical::pi*(sqr(pREff) - sqr(r_PW_mag))
+lagrangian/intermediate/submodels/Kinematic/InjectionModel/CellZoneInjection/CellZoneInjection.C:    this->volumeTotal_ = sum(pow3(diameters_))*constant::mathematical::pi/6.0;
+lagrangian/intermediate/submodels/Kinematic/InjectionModel/ConeNozzleInjection/ConeNozzleInjection.C:    const scalar deg2Rad = mathematical::pi/180.0;
+lagrangian/intermediate/submodels/Kinematic/InjectionModel/ConeNozzleInjection/ConeNozzleInjection.C:            scalar Ao = 0.25*mathematical::pi*outerDiameter_*outerDiameter_;
+lagrangian/intermediate/submodels/Kinematic/InjectionModel/ConeNozzleInjection/ConeNozzleInjection.C:            scalar Ai = 0.25*mathematical::pi*innerDiameter_*innerDiameter_;
+lagrangian/intermediate/submodels/Kinematic/SurfaceFilmModel/SurfaceFilmModel/SurfaceFilmModel.C:    scalar vol = mathematical::pi/6.0*pow3(diameterParcelPatch_[filmFaceI]);
+lagrangian/intermediate/submodels/Thermodynamic/ParticleForces/BrownianMotion/BrownianMotionForce.C:    scalar k = 2.0/(mathematical::pi*a) +  0.5*log(1.0 - y*y);
+lagrangian/intermediate/submodels/Thermodynamic/ParticleForces/BrownianMotion/BrownianMotionForce.C:        const scalar Dp = sigma*Tc*cc/(3*mathematical::pi*muc*dp);
+lagrangian/intermediate/submodels/Thermodynamic/ParticleForces/BrownianMotion/BrownianMotionForce.C:            216*muc*sigma*Tc/(sqr(mathematical::pi)*pow5(dp)*(rhoRatio)*cc);
+lagrangian/intermediate/submodels/Thermodynamic/ParticleForces/BrownianMotion/BrownianMotionForce.C:        f = eta*sqrt(mathematical::pi*s0/dt);
+lagrangian/molecularDynamics/molecule/molecule/moleculeI.H:inline const Foam::vector& Foam::molecule::pi() const
+lagrangian/molecularDynamics/molecule/molecule/moleculeI.H:inline Foam::vector& Foam::molecule::pi()
+lagrangian/molecularDynamics/potential/electrostaticPotential/electrostaticPotential.C:    prefactor(1.0/(4.0*constant::mathematical::pi*8.854187817e-12))
+lagrangian/molecularDynamics/potential/pairPotential/derived/coulomb/coulomb.C:    1.0/(4.0*constant::mathematical::pi*8.854187817e-12);
+lagrangian/molecularDynamics/potential/pairPotential/derived/dampedCoulomb/dampedCoulomb.C:    1.0/(4.0*constant::mathematical::pi*8.854187817e-12);
+lagrangian/molecularDynamics/potential/tetherPotential/derived/pitchForkRing/pitchForkRing.C:pitchForkRing::pitchForkRing
+lagrangian/molecularDynamics/potential/tetherPotential/derived/pitchForkRing/pitchForkRing.H:    Foam::tetherPotentials::pitchForkRing
+lagrangian/spray/submodels/AtomizationModel/BlobsSheetAtomization/BlobsSheetAtomization.C:            rho*sigma*d*cos(angle_*constant::mathematical::pi/360.0)
+lagrangian/spray/submodels/AtomizationModel/LISAAtomization/LISAAtomization.C:    scalar hSheet = volFlowRate/(constant::mathematical::pi*delta*Urel);
+lagrangian/spray/submodels/AtomizationModel/LISAAtomization/LISAAtomization.C:    scalar dD = cbrt(3.0*constant::mathematical::pi*sqr(dL)/kL);
+lagrangian/spray/submodels/AtomizationModel/LISAAtomization/LISASMDCalcMethod2.H:    // scalar factorGamma = 0.75*sqrt(mathematicalConstant::pi);   //nExp=2
+lagrangian/spray/submodels/BreakupModel/ReitzKHRT/ReitzKHRT.C:    scalar rhopi6 = rho*constant::mathematical::pi/6.0;
+lagrangian/spray/submodels/BreakupModel/ReitzKHRT/ReitzKHRT.C:            //scalar ms0 = rho*pow3(dc)*mathematicalConstant::pi/6.0;
+lagrangian/spray/submodels/BreakupModel/SHF/SHF.C:    scalar rhopi6 = rho*constant::mathematical::pi/6.0;
+lagrangian/spray/submodels/StochasticCollision/ORourkeCollision/ORourkeCollision.C:    scalar nu0 = 0.25*constant::mathematical::pi*sqr(sumD)*magUrel*dt/Vc;
+meshTools/coordinateSystems/coordinateRotation/EulerCoordinateRotation.C:        phi   *= constant::mathematical::pi/180.0;
+meshTools/coordinateSystems/coordinateRotation/EulerCoordinateRotation.C:        theta *= constant::mathematical::pi/180.0;
+meshTools/coordinateSystems/coordinateRotation/EulerCoordinateRotation.C:        psi   *= constant::mathematical::pi/180.0;
+meshTools/coordinateSystems/coordinateRotation/STARCDCoordinateRotation.C:        x *= constant::mathematical::pi/180.0;
+meshTools/coordinateSystems/coordinateRotation/STARCDCoordinateRotation.C:        y *= constant::mathematical::pi/180.0;
+meshTools/coordinateSystems/coordinateRotation/STARCDCoordinateRotation.C:        z *= constant::mathematical::pi/180.0;
+meshTools/coordinateSystems/cylindricalCS.C:        local.y()*(inDegrees_ ? constant::mathematical::pi/180.0 : 1.0)
+meshTools/coordinateSystems/cylindricalCS.C:       *(inDegrees_ ? constant::mathematical::pi/180.0 : 1.0)
+meshTools/coordinateSystems/cylindricalCS.C:        )*(inDegrees_ ? 180.0/constant::mathematical::pi : 1.0),
+meshTools/coordinateSystems/cylindricalCS.C:        )*(inDegrees_ ? 180.0/constant::mathematical::pi : 1.0)
+postProcessing/functionObjects/field/regionSizeDistribution/regionSizeDistribution.C:                   *6/constant::mathematical::pi
+postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionI.H:inline const Foam::vector& Foam::sixDoFRigidBodyMotion::pi() const
+postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionI.H:inline Foam::vector& Foam::sixDoFRigidBodyMotion::pi()
+postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionStateI.H:inline const Foam::vector& Foam::sixDoFRigidBodyMotionState::pi() const
+postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionStateI.H:inline Foam::vector& Foam::sixDoFRigidBodyMotionState::pi()
+randomProcesses/fft/kShellIntegration.C:    y *= sqr(x)*4.0*constant::mathematical::pi;
+randomProcesses/fft/kShellIntegration.C:    scalar factor = pow((l0/(2.0*constant::mathematical::pi)),3.0);
+regionModels/surfaceFilmModels/submodels/kinematic/injectionModel/drippingInjection/drippingInjection.C:    const scalar pi = constant::mathematical::pi;
+sampling/sampledSet/circle/circleSet.C:    const scalar alpha = constant::mathematical::pi/180.0*dTheta_;
+thermophysicalModels/radiationModels/radiationModel/fvDOM/radiativeIntensityRay/radiativeIntensityRay.C:                1.0/constant::mathematical::pi*omega_
+thermophysicalModels/radiationModels/radiationModel/fvDOM/radiativeIntensityRay/radiativeIntensityRay.C:               1.0/constant::mathematical::pi*omega_
+transportModels/interfaceProperties/interfaceProperties.C:    Foam::constant::mathematical::pi/180.0;
diff --git a/tutorials/compressible/rhoCentralFoam/LadenburgJet60psi/0.org/p b/tutorials/compressible/rhoCentralFoam/LadenburgJet60psi/0.org/p
index 19974b798e3a022a09c16ebe68ab8535bbc71719..bf63a32f63b4218ce1616722a73f1367215000a5 100644
--- a/tutorials/compressible/rhoCentralFoam/LadenburgJet60psi/0.org/p
+++ b/tutorials/compressible/rhoCentralFoam/LadenburgJet60psi/0.org/p
@@ -32,7 +32,7 @@ boundaryField
         field           p;
         phi             phi;
         rho             rho;
-        psi             psi;
+        psi             thermo:psi;
         fieldInf        101325;
         gamma           1.4;
         lInf            0.025;
@@ -47,7 +47,7 @@ boundaryField
         U               U;
         phi             phi;
         rho             none;
-        psi             psi;
+        psi             thermo:psi;
         gamma           1.4;
     }
 
diff --git a/tutorials/multiphase/twoPhaseEulerFoam/LES/bubbleColumn/constant/phaseProperties b/tutorials/multiphase/twoPhaseEulerFoam/LES/bubbleColumn/constant/phaseProperties
index 445b31ed27b4c066e571fb74975243250710e8f6..591ebc005ac39ad60b1c023764ac7d66914caaf1 100644
--- a/tutorials/multiphase/twoPhaseEulerFoam/LES/bubbleColumn/constant/phaseProperties
+++ b/tutorials/multiphase/twoPhaseEulerFoam/LES/bubbleColumn/constant/phaseProperties
@@ -36,6 +36,9 @@ water
     }
 }
 
+// Surface tension coefficient
+sigma           0.07;
+
 drag
 {
     air     SchillerNaumann;
diff --git a/tutorials/multiphase/twoPhaseEulerFoam/RAS/bubbleColumn/constant/phaseProperties b/tutorials/multiphase/twoPhaseEulerFoam/RAS/bubbleColumn/constant/phaseProperties
index 445b31ed27b4c066e571fb74975243250710e8f6..591ebc005ac39ad60b1c023764ac7d66914caaf1 100644
--- a/tutorials/multiphase/twoPhaseEulerFoam/RAS/bubbleColumn/constant/phaseProperties
+++ b/tutorials/multiphase/twoPhaseEulerFoam/RAS/bubbleColumn/constant/phaseProperties
@@ -36,6 +36,9 @@ water
     }
 }
 
+// Surface tension coefficient
+sigma           0.07;
+
 drag
 {
     air     SchillerNaumann;
diff --git a/tutorials/multiphase/twoPhaseEulerFoam/RAS/fluidisedBed/constant/phaseProperties b/tutorials/multiphase/twoPhaseEulerFoam/RAS/fluidisedBed/constant/phaseProperties
index 7d994cb3fd15234bf6cd8674e841b546b691b1c1..cb95638812c8df4060622783826d56924dcc110f 100644
--- a/tutorials/multiphase/twoPhaseEulerFoam/RAS/fluidisedBed/constant/phaseProperties
+++ b/tutorials/multiphase/twoPhaseEulerFoam/RAS/fluidisedBed/constant/phaseProperties
@@ -35,6 +35,9 @@ air
     }
 }
 
+// Surface tension coefficient
+sigma           0;
+
 drag
 {
     particles   GidaspowErgunWenYu;
diff --git a/tutorials/multiphase/twoPhaseEulerFoam/laminar/bubbleColumn/constant/phaseProperties b/tutorials/multiphase/twoPhaseEulerFoam/laminar/bubbleColumn/constant/phaseProperties
index 8e7701bf4cbb55811305ff6b022674cfd2bc6e95..b7325116e31f9bcf61fd7177b45779058dd024ce 100644
--- a/tutorials/multiphase/twoPhaseEulerFoam/laminar/bubbleColumn/constant/phaseProperties
+++ b/tutorials/multiphase/twoPhaseEulerFoam/laminar/bubbleColumn/constant/phaseProperties
@@ -36,6 +36,9 @@ water
     }
 }
 
+// Surface tension coefficient
+sigma           0.07;
+
 drag
 {
     air     SchillerNaumann;
diff --git a/tutorials/multiphase/twoPhaseEulerFoam/laminar/bubbleColumnIATE/0/T.air b/tutorials/multiphase/twoPhaseEulerFoam/laminar/bubbleColumnIATE/0/T.air
new file mode 100644
index 0000000000000000000000000000000000000000..9c0fd4206c19f34a05d8c05da28254b22882ebcf
--- /dev/null
+++ b/tutorials/multiphase/twoPhaseEulerFoam/laminar/bubbleColumnIATE/0/T.air
@@ -0,0 +1,45 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    object      Tair;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions          [0 0 0 1 0 0 0];
+
+internalField       uniform 300;
+
+boundaryField
+{
+    walls
+    {
+        type               zeroGradient;
+    }
+    outlet
+    {
+        type               inletOutlet;
+        phi                phi.air;
+        inletValue         $internalField;
+        value              $internalField;
+    }
+    inlet
+    {
+        type               fixedValue;
+        value              $internalField;
+    }
+    frontAndBackPlanes
+    {
+        type               empty;
+    }
+}
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
diff --git a/tutorials/multiphase/twoPhaseEulerFoam/laminar/bubbleColumnIATE/0/T.water b/tutorials/multiphase/twoPhaseEulerFoam/laminar/bubbleColumnIATE/0/T.water
new file mode 100644
index 0000000000000000000000000000000000000000..1cfd38f926516878085292090cf2e55699fef0cb
--- /dev/null
+++ b/tutorials/multiphase/twoPhaseEulerFoam/laminar/bubbleColumnIATE/0/T.water
@@ -0,0 +1,45 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    object      Twater;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions          [0 0 0 1 0 0 0];
+
+internalField       uniform 350;
+
+boundaryField
+{
+    walls
+    {
+        type               zeroGradient;
+    }
+    outlet
+    {
+        type               inletOutlet;
+        phi                phi.water;
+        inletValue         uniform 300;
+        value              $internalField;
+    }
+    inlet
+    {
+        type               fixedValue;
+        value              $internalField;
+    }
+    frontAndBackPlanes
+    {
+        type               empty;
+    }
+}
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
diff --git a/tutorials/multiphase/twoPhaseEulerFoam/laminar/bubbleColumnIATE/0/Theta b/tutorials/multiphase/twoPhaseEulerFoam/laminar/bubbleColumnIATE/0/Theta
new file mode 100644
index 0000000000000000000000000000000000000000..e45304b83464ea6a9568531b35b570abae0d768f
--- /dev/null
+++ b/tutorials/multiphase/twoPhaseEulerFoam/laminar/bubbleColumnIATE/0/Theta
@@ -0,0 +1,47 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    object      Theta;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions          [ 0 2 -2 0 0 0 0 ];
+
+internalField   uniform 0.0;
+
+boundaryField
+{
+    inlet
+    {
+        type            fixedValue;
+        value           uniform 1.0e-7;
+    }
+
+    outlet
+    {
+        type            inletOutlet;
+        inletValue      uniform 1.0e-7;
+        value           uniform 1.0e-7;
+    }
+
+    walls
+    {
+        type            zeroGradient;
+    }
+
+    defaultFaces
+    {
+        type            empty;
+    }
+}
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/twoPhaseEulerFoam/laminar/bubbleColumnIATE/0/U.air b/tutorials/multiphase/twoPhaseEulerFoam/laminar/bubbleColumnIATE/0/U.air
new file mode 100644
index 0000000000000000000000000000000000000000..e81fffac0612e0a20d26f1cb9be4211ba4be9a16
--- /dev/null
+++ b/tutorials/multiphase/twoPhaseEulerFoam/laminar/bubbleColumnIATE/0/U.air
@@ -0,0 +1,41 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      binary;
+    class       volVectorField;
+    object      U.air;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 1 -1 0 0 0 0];
+
+internalField   uniform (0 0.1 0);
+
+boundaryField
+{
+    inlet
+    {
+        type               fixedValue;
+        value              $internalField;
+    }
+    outlet
+    {
+        type               pressureInletOutletVelocity;
+        phi                phi.air;
+        value              $internalField;
+    }
+    walls
+    {
+        type               fixedValue;
+        value              uniform (0 0 0);
+    }
+}
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/twoPhaseEulerFoam/laminar/bubbleColumnIATE/0/U.water b/tutorials/multiphase/twoPhaseEulerFoam/laminar/bubbleColumnIATE/0/U.water
new file mode 100644
index 0000000000000000000000000000000000000000..aab00fd78bb2a097604737fb5b1b77d2ae159967
--- /dev/null
+++ b/tutorials/multiphase/twoPhaseEulerFoam/laminar/bubbleColumnIATE/0/U.water
@@ -0,0 +1,41 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      binary;
+    class       volVectorField;
+    object      U.water;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 1 -1 0 0 0 0];
+
+internalField   uniform (0 0 0);
+
+boundaryField
+{
+    inlet
+    {
+        type               fixedValue;
+        value              $internalField;
+    }
+    outlet
+    {
+        type               pressureInletOutletVelocity;
+        phi                phi.water;
+        value              $internalField;
+    }
+    walls
+    {
+        type               fixedValue;
+        value              uniform (0 0 0);
+    }
+}
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/twoPhaseEulerFoam/laminar/bubbleColumnIATE/0/alpha.air b/tutorials/multiphase/twoPhaseEulerFoam/laminar/bubbleColumnIATE/0/alpha.air
new file mode 100644
index 0000000000000000000000000000000000000000..1b1a35684a06787a1eb3cb75f388563aad1443a1
--- /dev/null
+++ b/tutorials/multiphase/twoPhaseEulerFoam/laminar/bubbleColumnIATE/0/alpha.air
@@ -0,0 +1,1926 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    location    "0";
+    object      alpha.air;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 0 0 0 0 0 0];
+
+internalField   nonuniform List<scalar>
+1875
+(
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+)
+;
+
+boundaryField
+{
+    inlet
+    {
+        type            fixedValue;
+        value           uniform 0.5;
+    }
+    outlet
+    {
+        type            inletOutlet;
+        phi             phi.air;
+        inletValue      uniform 1;
+        value           uniform 1;
+    }
+    walls
+    {
+        type            zeroGradient;
+    }
+    defaultFaces
+    {
+        type            empty;
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/twoPhaseEulerFoam/laminar/bubbleColumnIATE/0/alpha.air.org b/tutorials/multiphase/twoPhaseEulerFoam/laminar/bubbleColumnIATE/0/alpha.air.org
new file mode 100644
index 0000000000000000000000000000000000000000..4472b0c63598b7f95acf8d79fac10b3e9a08ebf7
--- /dev/null
+++ b/tutorials/multiphase/twoPhaseEulerFoam/laminar/bubbleColumnIATE/0/alpha.air.org
@@ -0,0 +1,42 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    location    "0";
+    object      alpha.air;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 0 0 0 0 0 0];
+
+internalField   uniform 0;
+
+boundaryField
+{
+    inlet
+    {
+        type            fixedValue;
+        value           uniform 0.5;
+    }
+    outlet
+    {
+        type            inletOutlet;
+        phi             phi.air;
+        inletValue      uniform 1;
+        value           uniform 1;
+    }
+    walls
+    {
+        type            zeroGradient;
+    }
+}
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/twoPhaseEulerFoam/laminar/bubbleColumnIATE/0/kappai.air b/tutorials/multiphase/twoPhaseEulerFoam/laminar/bubbleColumnIATE/0/kappai.air
new file mode 100644
index 0000000000000000000000000000000000000000..cd5560fb91c54e17ae70b91f32c52ef8fde3aa56
--- /dev/null
+++ b/tutorials/multiphase/twoPhaseEulerFoam/laminar/bubbleColumnIATE/0/kappai.air
@@ -0,0 +1,45 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    object      kappai.air;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions          [0 -1 0 0 0 0 0];
+
+internalField       uniform 2000;
+
+boundaryField
+{
+    walls
+    {
+        type               zeroGradient;
+    }
+    outlet
+    {
+        type               inletOutlet;
+        phi                phi.air;
+        inletValue         $internalField;
+        value              $internalField;
+    }
+    inlet
+    {
+        type               fixedValue;
+        value              $internalField;
+    }
+    frontAndBackPlanes
+    {
+        type               empty;
+    }
+}
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
diff --git a/tutorials/multiphase/twoPhaseEulerFoam/laminar/bubbleColumnIATE/0/p b/tutorials/multiphase/twoPhaseEulerFoam/laminar/bubbleColumnIATE/0/p
new file mode 100644
index 0000000000000000000000000000000000000000..ae2c14b7460e91d82d3717d540b6a9d33dc087e8
--- /dev/null
+++ b/tutorials/multiphase/twoPhaseEulerFoam/laminar/bubbleColumnIATE/0/p
@@ -0,0 +1,40 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    object      p;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions          [ 1 -1 -2 0 0 0 0 ];
+
+internalField       uniform 1e5;
+
+boundaryField
+{
+    inlet
+    {
+        type               fixedFluxPressure;
+        value              $internalField;
+    }
+    outlet
+    {
+        type               fixedValue;
+        value              $internalField;
+    }
+    walls
+    {
+        type               fixedFluxPressure;
+        value              $internalField;
+    }
+}
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
diff --git a/tutorials/multiphase/twoPhaseEulerFoam/laminar/bubbleColumnIATE/constant/g b/tutorials/multiphase/twoPhaseEulerFoam/laminar/bubbleColumnIATE/constant/g
new file mode 100644
index 0000000000000000000000000000000000000000..0cc222ca3457ed24bf9753d0926fbee84359e624
--- /dev/null
+++ b/tutorials/multiphase/twoPhaseEulerFoam/laminar/bubbleColumnIATE/constant/g
@@ -0,0 +1,22 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       uniformDimensionedVectorField;
+    location    "constant";
+    object      g;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 1 -2 0 0 0 0];
+value           (0 -9.81 0);
+
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/twoPhaseEulerFoam/laminar/bubbleColumnIATE/constant/phaseProperties b/tutorials/multiphase/twoPhaseEulerFoam/laminar/bubbleColumnIATE/constant/phaseProperties
new file mode 100644
index 0000000000000000000000000000000000000000..d9304132f75d573a0a262e64d6cba880d083b95c
--- /dev/null
+++ b/tutorials/multiphase/twoPhaseEulerFoam/laminar/bubbleColumnIATE/constant/phaseProperties
@@ -0,0 +1,90 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "constant";
+    object      phaseProperties;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+phases (air water);
+
+air
+{
+    diameterModel   IATE;
+
+    IATECoeffs
+    {
+        dMax 1e-2;
+        dMin 1e-4;
+
+        sources
+        (
+            wakeEntrainmentCoalescence
+            {
+                Cwe         0.002;
+            }
+            randomCoalescence
+            {
+                Crc         0.04;
+                C           3;
+                alphaMax    0.75;
+            }
+            turbulentBreakUp
+            {
+                Cti         0.085;
+                WeCr        6;
+            }
+        );
+    }
+}
+
+water
+{
+    diameterModel constant;
+    constantCoeffs
+    {
+        d               1e-4;
+    }
+}
+
+// Surface tension coefficient
+sigma           0.07;
+
+drag
+{
+    air     SchillerNaumann;
+    water   SchillerNaumann;
+}
+
+heatTransfer
+{
+    air     RanzMarshall;
+    water   RanzMarshall;
+}
+
+dispersedPhase          both;
+
+residualPhaseFraction   1e-3;
+residualSlip            1e-2;
+
+// Virtual-mass ceofficient
+Cvm             0.5;
+
+// Lift coefficient
+Cl              0;
+
+// Minimum allowable pressure
+pMin            10000;
+
+
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/twoPhaseEulerFoam/laminar/bubbleColumnIATE/constant/polyMesh/blockMeshDict b/tutorials/multiphase/twoPhaseEulerFoam/laminar/bubbleColumnIATE/constant/polyMesh/blockMeshDict
new file mode 100644
index 0000000000000000000000000000000000000000..d03967afdc2ccf7afbecdf32d50159c309a475e3
--- /dev/null
+++ b/tutorials/multiphase/twoPhaseEulerFoam/laminar/bubbleColumnIATE/constant/polyMesh/blockMeshDict
@@ -0,0 +1,61 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    object      blockMeshDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+convertToMeters 1;
+
+vertices
+(
+    (0 0 0)
+    (0.15 0 0)
+    (0.15 1 0)
+    (0 1 0)
+    (0 0 0.1)
+    (0.15 0 0.1)
+    (0.15 1 0.1)
+    (0 1 0.1)
+);
+
+blocks
+(
+    hex (0 1 2 3 4 5 6 7) (25 75 1) simpleGrading (1 1 1)
+);
+
+edges
+(
+);
+
+patches
+(
+    patch inlet
+    (
+        (1 5 4 0)
+    )
+    patch outlet
+    (
+        (3 7 6 2)
+    )
+    wall walls
+    (
+        (0 4 7 3)
+        (2 6 5 1)
+    )
+);
+
+mergePatchPairs
+(
+);
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/twoPhaseEulerFoam/laminar/bubbleColumnIATE/constant/polyMesh/boundary b/tutorials/multiphase/twoPhaseEulerFoam/laminar/bubbleColumnIATE/constant/polyMesh/boundary
new file mode 100644
index 0000000000000000000000000000000000000000..bf47f69643c9925d3a1ef19c6b4ddc67cf604e0a
--- /dev/null
+++ b/tutorials/multiphase/twoPhaseEulerFoam/laminar/bubbleColumnIATE/constant/polyMesh/boundary
@@ -0,0 +1,47 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       polyBoundaryMesh;
+    location    "constant/polyMesh";
+    object      boundary;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+4
+(
+    inlet
+    {
+        type            patch;
+        nFaces          25;
+        startFace       3650;
+    }
+    outlet
+    {
+        type            patch;
+        nFaces          25;
+        startFace       3675;
+    }
+    walls
+    {
+        type            wall;
+        nFaces          150;
+        startFace       3700;
+    }
+    defaultFaces
+    {
+        type            empty;
+        inGroups        1(empty);
+        nFaces          3750;
+        startFace       3850;
+    }
+)
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/twoPhaseEulerFoam/laminar/bubbleColumnIATE/constant/thermophysicalProperties.air b/tutorials/multiphase/twoPhaseEulerFoam/laminar/bubbleColumnIATE/constant/thermophysicalProperties.air
new file mode 100644
index 0000000000000000000000000000000000000000..3cac781ab1bf649afe060310f5ea7c898123efd9
--- /dev/null
+++ b/tutorials/multiphase/twoPhaseEulerFoam/laminar/bubbleColumnIATE/constant/thermophysicalProperties.air
@@ -0,0 +1,49 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "constant";
+    object      thermophysicalProperties.air;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+thermoType
+{
+    type            heRhoThermo;
+    mixture         pureMixture;
+    transport       const;
+    thermo          hConst;
+    equationOfState perfectGas;
+    specie          specie;
+    energy          sensibleEnthalpy;
+}
+
+mixture
+{
+    specie
+    {
+        nMoles      1;
+        molWeight   28.9;
+    }
+    thermodynamics
+    {
+        Cp          1007;
+        Hf          0;
+    }
+    transport
+    {
+        mu          1.84e-05;
+        Pr          0.7;
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/twoPhaseEulerFoam/laminar/bubbleColumnIATE/constant/thermophysicalProperties.water b/tutorials/multiphase/twoPhaseEulerFoam/laminar/bubbleColumnIATE/constant/thermophysicalProperties.water
new file mode 100644
index 0000000000000000000000000000000000000000..394eb31679cf250cd4f43a40924f257c503e52c1
--- /dev/null
+++ b/tutorials/multiphase/twoPhaseEulerFoam/laminar/bubbleColumnIATE/constant/thermophysicalProperties.water
@@ -0,0 +1,54 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "constant";
+    object      thermophysicalProperties.water;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+thermoType
+{
+    type            heRhoThermo;
+    mixture         pureMixture;
+    transport       const;
+    thermo          hConst;
+    equationOfState perfectFluid;
+    specie          specie;
+    energy          sensibleEnthalpy;
+}
+
+mixture
+{
+    specie
+    {
+        nMoles      1;
+        molWeight   18;
+    }
+    equationOfState
+    {
+        R           3000;
+        rho0        1027;
+    }
+    thermodynamics
+    {
+        Cp          4195;
+        Hf          0;
+    }
+    transport
+    {
+        mu          3.645e-4;
+        Pr          2.289;
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/twoPhaseEulerFoam/laminar/bubbleColumnIATE/constant/turbulenceProperties.air b/tutorials/multiphase/twoPhaseEulerFoam/laminar/bubbleColumnIATE/constant/turbulenceProperties.air
new file mode 100644
index 0000000000000000000000000000000000000000..1296429b72a21953def920b08774aa75e1d048b1
--- /dev/null
+++ b/tutorials/multiphase/twoPhaseEulerFoam/laminar/bubbleColumnIATE/constant/turbulenceProperties.air
@@ -0,0 +1,20 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "constant";
+    object      turbulenceProperties.air;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+simulationType  laminar;
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/twoPhaseEulerFoam/laminar/bubbleColumnIATE/constant/turbulenceProperties.water b/tutorials/multiphase/twoPhaseEulerFoam/laminar/bubbleColumnIATE/constant/turbulenceProperties.water
new file mode 100644
index 0000000000000000000000000000000000000000..7f0d75a82fcdc99677fa0be8a4111cfe91e4a82c
--- /dev/null
+++ b/tutorials/multiphase/twoPhaseEulerFoam/laminar/bubbleColumnIATE/constant/turbulenceProperties.water
@@ -0,0 +1,20 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "constant";
+    object      turbulenceProperties.water;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+simulationType  laminar;
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/twoPhaseEulerFoam/laminar/bubbleColumnIATE/system/controlDict b/tutorials/multiphase/twoPhaseEulerFoam/laminar/bubbleColumnIATE/system/controlDict
new file mode 100644
index 0000000000000000000000000000000000000000..0a0bfc89cf3126321d001724e2eeae4bcb69e8aa
--- /dev/null
+++ b/tutorials/multiphase/twoPhaseEulerFoam/laminar/bubbleColumnIATE/system/controlDict
@@ -0,0 +1,95 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "system";
+    object      controlDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+application     twoPhaseEulerFoam;
+
+startFrom       startTime;
+
+startTime       0;
+
+stopAt          endTime;
+
+endTime         100;
+
+deltaT          0.005;
+
+writeControl    runTime;
+
+writeInterval   1;
+
+purgeWrite      0;
+
+writeFormat     ascii;
+
+writePrecision  6;
+
+writeCompression uncompressed;
+
+timeFormat      general;
+
+timePrecision   6;
+
+runTimeModifiable yes;
+
+adjustTimeStep  no;
+
+maxCo           0.5;
+
+maxDeltaT       1;
+
+functions
+{
+    fieldAverage1
+    {
+        type            fieldAverage;
+        functionObjectLibs ( "libfieldFunctionObjects.so" );
+        outputControl   outputTime;
+        fields
+        (
+            U.air
+            {
+                 mean        on;
+                 prime2Mean  off;
+                 base        time;
+            }
+
+            U.water
+            {
+                 mean        on;
+                 prime2Mean  off;
+                 base        time;
+            }
+
+            alpha.air
+            {
+                 mean        on;
+                 prime2Mean  off;
+                 base        time;
+            }
+
+            p
+            {
+                 mean        on;
+                 prime2Mean  off;
+                 base        time;
+            }
+        );
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/twoPhaseEulerFoam/laminar/bubbleColumnIATE/system/fvSchemes b/tutorials/multiphase/twoPhaseEulerFoam/laminar/bubbleColumnIATE/system/fvSchemes
new file mode 100644
index 0000000000000000000000000000000000000000..3fd8876cb600ec235534eb44522c437e989fb8df
--- /dev/null
+++ b/tutorials/multiphase/twoPhaseEulerFoam/laminar/bubbleColumnIATE/system/fvSchemes
@@ -0,0 +1,69 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "system";
+    object      fvSchemes;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+ddtSchemes
+{
+    default         Euler;
+}
+
+gradSchemes
+{
+    default         Gauss linear;
+}
+
+divSchemes
+{
+    default                     none;
+
+    div(phi,alpha.air)          Gauss vanLeer;
+    div(phir,alpha.air)         Gauss vanLeer;
+
+    div(phi.air,kappai.air)     Gauss vanLeer;
+
+    "div\(alphaPhi.*,U.*\)"     Gauss limitedLinearV 1;
+    "div\(phi.*,U.*\)"          Gauss limitedLinearV 1;
+    "div\(phi.*,.*rho.*\)"      Gauss linear;
+
+    "div\(alphaPhi.*,(h|e).*\)" Gauss limitedLinear 1;
+    "div\(alphaPhi.*,(K.*|p)\)" Gauss limitedLinear 1;
+
+    "div\(\(\(alpha.*nuEff.*\)*dev2\(T\(grad\(U.*\)\)\)\)\)" Gauss linear;
+}
+
+laplacianSchemes
+{
+    default         Gauss linear uncorrected;
+}
+
+interpolationSchemes
+{
+    default         linear;
+}
+
+snGradSchemes
+{
+    default         uncorrected;
+}
+
+fluxRequired
+{
+    default         no;
+    p               ;
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/twoPhaseEulerFoam/laminar/bubbleColumnIATE/system/fvSolution b/tutorials/multiphase/twoPhaseEulerFoam/laminar/bubbleColumnIATE/system/fvSolution
new file mode 100644
index 0000000000000000000000000000000000000000..0a266db563c41fbb5c61dc5c9ad62e8713d31c06
--- /dev/null
+++ b/tutorials/multiphase/twoPhaseEulerFoam/laminar/bubbleColumnIATE/system/fvSolution
@@ -0,0 +1,84 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "system";
+    object      fvSolution;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+solvers
+{
+    alpha.air
+    {
+        nAlphaCorr      1;
+        nAlphaSubCycles 2;
+    }
+
+    p
+    {
+        solver          GAMG;
+        smoother        DIC;
+        nPreSweeps      0;
+        nPostSweeps     2;
+        nFinestSweeps   2;
+        cacheAgglomeration true;
+        nCellsInCoarsestLevel 10;
+        agglomerator    faceAreaPair;
+        mergeLevels     1;
+        tolerance       1e-08;
+        relTol          0.01;
+    }
+
+    pFinal
+    {
+        $p;
+        tolerance       1e-08;
+        relTol          0;
+    }
+
+    "(U|kappai).*"
+    {
+        solver          PBiCG;
+        preconditioner  DILU;
+        tolerance       1e-05;
+        relTol          0;
+    }
+
+    "h.*"
+    {
+        solver          PBiCG;
+        preconditioner  DILU;
+        tolerance       1e-8;
+        relTol          0;
+    }
+}
+
+PIMPLE
+{
+    nOuterCorrectors 1;
+    nCorrectors     2;
+    nNonOrthogonalCorrectors 0;
+}
+
+relaxationFactors
+{
+    fields
+    {
+    }
+    equations
+    {
+        ".*"            1;
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/twoPhaseEulerFoam/laminar/bubbleColumnIATE/system/setFieldsDict b/tutorials/multiphase/twoPhaseEulerFoam/laminar/bubbleColumnIATE/system/setFieldsDict
new file mode 100644
index 0000000000000000000000000000000000000000..85996cf966762c6ca3a7c37a1eaa8ae462ecdb19
--- /dev/null
+++ b/tutorials/multiphase/twoPhaseEulerFoam/laminar/bubbleColumnIATE/system/setFieldsDict
@@ -0,0 +1,36 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "system";
+    object      setFieldsDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+defaultFieldValues
+(
+    volScalarFieldValue alpha1 1
+);
+
+regions
+(
+    boxToCell
+    {
+        box (0 0 -0.1) (0.15 0.701 0.1);
+        fieldValues
+        (
+            volScalarFieldValue alphaair 0
+        );
+    }
+);
+
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/twoPhaseEulerFoam/laminar/fluidisedBed/constant/phaseProperties b/tutorials/multiphase/twoPhaseEulerFoam/laminar/fluidisedBed/constant/phaseProperties
index 7d994cb3fd15234bf6cd8674e841b546b691b1c1..cb95638812c8df4060622783826d56924dcc110f 100644
--- a/tutorials/multiphase/twoPhaseEulerFoam/laminar/fluidisedBed/constant/phaseProperties
+++ b/tutorials/multiphase/twoPhaseEulerFoam/laminar/fluidisedBed/constant/phaseProperties
@@ -35,6 +35,9 @@ air
     }
 }
 
+// Surface tension coefficient
+sigma           0;
+
 drag
 {
     particles   GidaspowErgunWenYu;
diff --git a/tutorials/multiphase/twoPhaseEulerFoam/laminar/mixerVessel2D/constant/phaseProperties b/tutorials/multiphase/twoPhaseEulerFoam/laminar/mixerVessel2D/constant/phaseProperties
index 0046495f912257f66dc73da800577ee976f73f42..f867970c9afab3871b6b24bf0450c43f712df6d7 100644
--- a/tutorials/multiphase/twoPhaseEulerFoam/laminar/mixerVessel2D/constant/phaseProperties
+++ b/tutorials/multiphase/twoPhaseEulerFoam/laminar/mixerVessel2D/constant/phaseProperties
@@ -36,6 +36,9 @@ water
     }
 }
 
+// Surface tension coefficient
+sigma           0.07;
+
 drag
 {
     air     SchillerNaumann;