diff --git a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/EEqns.H b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/EEqns.H
new file mode 100644
index 0000000000000000000000000000000000000000..c894b36b41c4e48a07a650590f9c6fe6fdc7edc2
--- /dev/null
+++ b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/EEqns.H
@@ -0,0 +1,49 @@
+{
+    volScalarField k1
+    (
+        "k1",
+        alpha1*(thermo1.alpha()/rho1 + sqr(Ct)*nut2*thermo1.CpByCpv()/Prt)
+    );
+
+    volScalarField k2
+    (
+        "k2",
+        alpha2*(thermo2.alpha()/rho2 + nut2*thermo2.CpByCpv()/Prt)
+    );
+
+    volScalarField& he1 = thermo1.he();
+    volScalarField& he2 = thermo2.he();
+
+    Info<< max(he1) << min(he1) << endl;
+
+    fvScalarMatrix he1Eqn
+    (
+        fvm::ddt(alpha1, he1)
+      + fvm::div(alphaPhi1, he1)
+      - fvm::laplacian(k1, he1)
+     ==
+        heatTransferCoeff*(he1/thermo1.Cp())/rho1
+      - fvm::Sp(heatTransferCoeff/thermo1.Cp()/rho1, he1)
+      + alpha1*(dpdt/rho1 - (fvc::ddt(K1) + fvc::div(phi1, K1)))
+    );
+
+    fvScalarMatrix he2Eqn
+    (
+        fvm::ddt(alpha2, he2)
+      + fvm::div(alphaPhi2, he2)
+      - fvm::laplacian(k2, he2)
+     ==
+        heatTransferCoeff*(he2/thermo2.Cp())/rho2
+      - fvm::Sp(heatTransferCoeff/thermo2.Cp()/rho2, he2)
+      + alpha2*(dpdt/rho2 - (fvc::ddt(K2) + fvc::div(phi2, K2)))
+    );
+
+    he1Eqn.relax();
+    he1Eqn.solve();
+
+    he2Eqn.relax();
+    he2Eqn.solve();
+
+    thermo1.correct();
+    thermo2.correct();
+}
diff --git a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/Make/options b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/Make/options
index 16ce2e3674b971c19291acbbfb292aeb2da3795e..d90ba7bbbdf308e0b43489c0e0efd851b5367e28 100644
--- a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/Make/options
+++ b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/Make/options
@@ -1,4 +1,5 @@
 EXE_INC = \
+    -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
     -I$(LIB_SRC)/transportModels/incompressible/lnInclude \
     -IturbulenceModel \
@@ -8,6 +9,8 @@ EXE_INC = \
     -Iaveraging
 
 EXE_LIBS = \
+    -lfluidThermophysicalModels \
+    -lspecie \
     -lincompressibleTransportModels \
     -lcompressiblePhaseModel \
     -lcompressibleEulerianInterfacialModels \
diff --git a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/TEqns.H b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/TEqns.H
deleted file mode 100644
index 8f38ca9d91fe910cfad67e53608612d277564b3a..0000000000000000000000000000000000000000
--- a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/TEqns.H
+++ /dev/null
@@ -1,36 +0,0 @@
-{
-    volScalarField kByCp1("kByCp1", alpha1*(k1/Cp1/rho1 + sqr(Ct)*nut2/Prt));
-    volScalarField kByCp2("kByCp2", alpha2*(k2/Cp2/rho2 + nut2/Prt));
-
-    fvScalarMatrix T1Eqn
-    (
-        fvm::ddt(alpha1, T1)
-      + fvm::div(alphaPhi1, T1)
-      - fvm::laplacian(kByCp1, T1)
-     ==
-        heatTransferCoeff*T2/Cp1/rho1
-      - fvm::Sp(heatTransferCoeff/Cp1/rho1, T1)
-      + alpha1*(dpdt/rho1 - (fvc::ddt(K1) + fvc::div(phi1, K1)))/Cp1
-    );
-
-    fvScalarMatrix T2Eqn
-    (
-        fvm::ddt(alpha2, T2)
-      + fvm::div(alphaPhi2, T2)
-      - fvm::laplacian(kByCp2, T2)
-     ==
-        heatTransferCoeff*T1/Cp2/rho2
-      - fvm::Sp(heatTransferCoeff/Cp2/rho2, T2)
-      + alpha2*(dpdt/rho2 - (fvc::ddt(K2) + fvc::div(phi2, K2)))/Cp2
-    );
-
-    T1Eqn.relax();
-    T1Eqn.solve();
-
-    T2Eqn.relax();
-    T2Eqn.solve();
-
-    // Update compressibilities
-    psi1 = 1.0/(R1*T1);
-    psi2 = 1.0/(R2*T2);
-}
diff --git a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/UEqns.H b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/UEqns.H
index 606b4573e47988c921954a011626b7ce80ecfb4b..a243db6300e555bfc8b909c74cd11352d7d66abd 100644
--- a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/UEqns.H
+++ b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/UEqns.H
@@ -12,12 +12,12 @@ fvVectorMatrix U2Eqn(U2, U2.dimensions()*dimVol/dimTime);
         }
         else // If not using kinetic theory is using Ct model
         {
-            nuEff1 = sqr(Ct)*nut2 + nu1;
+            nuEff1 = sqr(Ct)*nut2 + thermo1.mu()/rho1;
         }
 
         volTensorField Rc1
         (
-            "Rc1",
+            "Rc",
             (((2.0/3.0)*I)*nuEff1)*tr(gradU1T) - nuEff1*gradU1T
         );
 
@@ -52,7 +52,7 @@ fvVectorMatrix U2Eqn(U2, U2.dimensions()*dimVol/dimTime);
         volTensorField gradU2T(fvc::grad(U2)().T());
         volTensorField Rc2
         (
-            "Rc2",
+            "Rc",
             (((2.0/3.0)*I)*nuEff2)*tr(gradU2T) - nuEff2*gradU2T
         );
 
diff --git a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/alphaEqn.H b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/alphaEqn.H
index a818ee2e9ea7707a7707b25f40026a2bb08ad591..86d9203dce190cfe520a07149dccbe7ef66eab83 100644
--- a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/alphaEqn.H
+++ b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/alphaEqn.H
@@ -1,9 +1,9 @@
-surfaceScalarField alphaPhi1("alphaPhi1", phi1);
-surfaceScalarField alphaPhi2("alphaPhi2", phi2);
+surfaceScalarField alphaPhi1("alphaPhi", phi1);
+surfaceScalarField alphaPhi2("alphaPhi", phi2);
 
 {
-    word scheme("div(phi,alpha1)");
-    word schemer("div(phir,alpha1)");
+    word scheme("div(phi,alpha)");
+    word schemer("div(phir,alpha)");
 
     surfaceScalarField phic("phic", phi);
     surfaceScalarField phir("phir", phi1 - phi2);
@@ -13,7 +13,7 @@ surfaceScalarField alphaPhi2("alphaPhi2", phi2);
         surfaceScalarField alpha1f(fvc::interpolate(alpha1));
         surfaceScalarField phipp(ppMagf*fvc::snGrad(alpha1)*mesh.magSf());
         phir += phipp;
-        phic += fvc::interpolate(alpha1)*phipp;
+        phic += alpha1f*phipp;
     }
 
     for (int acorr=0; acorr<nAlphaCorr; acorr++)
@@ -68,16 +68,22 @@ surfaceScalarField alphaPhi2("alphaPhi2", phi2);
 
         if (g0.value() > 0.0)
         {
+            surfaceScalarField alpha1f(fvc::interpolate(alpha1));
+
             ppMagf =
                 fvc::interpolate((1.0/rho1)*rAU1)
-               *fvc::interpolate
-                (
-                    g0*min(exp(preAlphaExp*(alpha1 - alphaMax)), expMax)
-                );
+               *g0*min(exp(preAlphaExp*(alpha1f - alphaMax)), expMax);
+
+            // ppMagf =
+            //     fvc::interpolate((1.0/rho1)*rAU1)
+            //    *fvc::interpolate
+            //     (
+            //         g0*min(exp(preAlphaExp*(alpha1 - alphaMax)), expMax)
+            //     );
 
             alpha1Eqn -= fvm::laplacian
             (
-                (fvc::interpolate(alpha1) + scalar(0.0001))*ppMagf,
+                alpha1f*ppMagf,
                 alpha1,
                 "laplacian(alphaPpMag,alpha1)"
             );
diff --git a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/compressibleTwoPhaseEulerFoam.C b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/compressibleTwoPhaseEulerFoam.C
index 855a862b6ba51ac11c1fdadbae38e37f5fa35dc2..65372612b68532228120394c2bcc6e0925e44c2d 100644
--- a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/compressibleTwoPhaseEulerFoam.C
+++ b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/compressibleTwoPhaseEulerFoam.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -31,6 +31,7 @@ Description
 \*---------------------------------------------------------------------------*/
 
 #include "fvCFD.H"
+#include "rhoThermo.H"
 #include "nearWallDist.H"
 #include "wallFvPatch.H"
 #include "fixedValueFvsPatchFields.H"
@@ -86,7 +87,7 @@ int main(int argc, char *argv[])
             #include "alphaEqn.H"
             #include "kEpsilon.H"
             #include "interfacialCoeffs.H"
-            #include "TEqns.H"
+            #include "EEqns.H"
             #include "UEqns.H"
 
             // --- Pressure corrector loop
diff --git a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/createFields.H b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/createFields.H
index e0cd5328c98a5585040b805c5257cf676ab24996..6095ad6845e3a6b08767950eb9298bc27db2e416 100644
--- a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/createFields.H
+++ b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/createFields.H
@@ -12,55 +12,43 @@
         )
     );
 
+    word phase1Name
+    (
+        transportProperties.found("phases")
+      ? wordList(transportProperties.lookup("phases"))[0]
+      : "phase1"
+    );
+
+    word phase2Name
+    (
+        transportProperties.found("phases")
+      ? wordList(transportProperties.lookup("phases"))[1]
+      : "phase2"
+    );
+
     autoPtr<phaseModel> phase1 = phaseModel::New
     (
         mesh,
         transportProperties,
-        "1"
+        phase1Name
     );
 
     autoPtr<phaseModel> phase2 = phaseModel::New
     (
         mesh,
         transportProperties,
-        "2"
+        phase2Name
     );
 
+    volScalarField& alpha1 = phase1();
+    volScalarField& alpha2 = phase2();
+    alpha2 = scalar(1) - alpha1;
+
     volVectorField& U1 = phase1->U();
     surfaceScalarField& phi1 = phase1->phi();
-    const dimensionedScalar& nu1 = phase1->nu();
-    const dimensionedScalar& k1 = phase1->kappa();
-    const dimensionedScalar& Cp1 = phase1->Cp();
-    dimensionedScalar rho10
-    (
-        "rho0",
-        dimDensity,
-        transportProperties.subDict("phase1").lookup("rho0")
-    );
-    dimensionedScalar R1
-    (
-        "R",
-        dimensionSet(0, 2, -2, -1, 0),
-        transportProperties.subDict("phase1").lookup("R")
-    );
 
     volVectorField& U2 = phase2->U();
     surfaceScalarField& phi2 = phase2->phi();
-    const dimensionedScalar& nu2 = phase2->nu();
-    const dimensionedScalar& k2 = phase2->kappa();
-    const dimensionedScalar& Cp2 = phase2->Cp();
-    dimensionedScalar rho20
-    (
-        "rho0",
-        dimDensity,
-        transportProperties.subDict("phase2").lookup("rho0")
-    );
-    dimensionedScalar R2
-    (
-        "R",
-        dimensionSet(0, 2, -2, -1, 0),
-        transportProperties.subDict("phase2").lookup("R")
-    );
 
     dimensionedScalar pMin
     (
@@ -69,102 +57,16 @@
         transportProperties.lookup("pMin")
     );
 
+    rhoThermo& thermo1 = phase1->thermo();
+    rhoThermo& thermo2 = phase2->thermo();
 
-    Info<< "Reading field T1" << endl;
-    volScalarField T1
-    (
-        IOobject
-        (
-            "T1",
-            runTime.timeName(),
-            mesh,
-            IOobject::MUST_READ,
-            IOobject::AUTO_WRITE
-        ),
-        mesh
-    );
+    volScalarField& p = thermo1.p();
 
-    Info<< "Reading field T2" << endl;
-    volScalarField T2
-    (
-        IOobject
-        (
-            "T2",
-            runTime.timeName(),
-            mesh,
-            IOobject::MUST_READ,
-            IOobject::AUTO_WRITE
-        ),
-        mesh
-    );
-
-    volScalarField psi1
-    (
-        IOobject
-        (
-            "psi1",
-            runTime.timeName(),
-            mesh
-        ),
-        1.0/(R1*T1)
-    );
-
-    volScalarField psi2
-    (
-        IOobject
-        (
-            "psi2",
-            runTime.timeName(),
-            mesh
-        ),
-        1.0/(R2*T2)
-    );
-    psi2.oldTime();
-
-    Info<< "Reading field p\n" << endl;
-    volScalarField p
-    (
-        IOobject
-        (
-            "p",
-            runTime.timeName(),
-            mesh,
-            IOobject::MUST_READ,
-            IOobject::AUTO_WRITE
-        ),
-        mesh
-    );
-
-    volScalarField rho1("rho1", rho10 + psi1*p);
-    volScalarField rho2("rho2", rho20 + psi2*p);
-
-
-    Info<< "Reading field alpha1\n" << endl;
-    volScalarField alpha1
-    (
-        IOobject
-        (
-            "alpha1",
-            runTime.timeName(),
-            mesh,
-            IOobject::MUST_READ,
-            IOobject::AUTO_WRITE
-        ),
-        mesh
-    );
+    volScalarField rho1("rho" + phase1Name, thermo1.rho());
+    const volScalarField& psi1 = thermo1.psi();
 
-    volScalarField alpha2
-    (
-        IOobject
-        (
-            "alpha2",
-            runTime.timeName(),
-            mesh,
-            IOobject::NO_READ,
-            IOobject::NO_WRITE
-        ),
-        scalar(1) - alpha1
-    );
+    volScalarField rho2("rho" + phase2Name, thermo2.rho());
+    const volScalarField& psi2 = thermo2.psi();
 
     volVectorField U
     (
@@ -179,27 +81,6 @@
         alpha1*U1 + alpha2*U2
     );
 
-    dimensionedScalar Cvm
-    (
-        "Cvm",
-        dimless,
-        transportProperties.lookup("Cvm")
-    );
-
-    dimensionedScalar Cl
-    (
-        "Cl",
-        dimless,
-        transportProperties.lookup("Cl")
-    );
-
-    dimensionedScalar Ct
-    (
-        "Ct",
-        dimless,
-        transportProperties.lookup("Ct")
-    );
-
     surfaceScalarField phi
     (
         IOobject
@@ -226,8 +107,6 @@
         alpha1*rho1 + alpha2*rho2
     );
 
-    #include "createRASTurbulence.H"
-
     Info<< "Calculating field DDtU1 and DDtU2\n" << endl;
 
     volVectorField DDtU1
@@ -248,6 +127,29 @@
     Info<< "Calculating field g.h\n" << endl;
     volScalarField gh("gh", g & mesh.C());
 
+    dimensionedScalar Cvm
+    (
+        "Cvm",
+        dimless,
+        transportProperties.lookup("Cvm")
+    );
+
+    dimensionedScalar Cl
+    (
+        "Cl",
+        dimless,
+        transportProperties.lookup("Cl")
+    );
+
+    dimensionedScalar Ct
+    (
+        "Ct",
+        dimless,
+        transportProperties.lookup("Ct")
+    );
+
+    #include "createRASTurbulence.H"
+
     IOdictionary interfacialProperties
     (
         IOobject
@@ -297,8 +199,8 @@
     if
     (
         !(
-            dispersedPhase == "1"
-         || dispersedPhase == "2"
+            dispersedPhase == phase1Name
+         || dispersedPhase == phase2Name
          || dispersedPhase == "both"
         )
     )
@@ -337,7 +239,7 @@
     (
         IOobject
         (
-            "rAU1",
+            "rAU" + phase1Name,
             runTime.timeName(),
             mesh,
             IOobject::NO_READ,
@@ -387,5 +289,5 @@
     );
 
     Info<< "Creating field kinetic energy K\n" << endl;
-    volScalarField K1("K1", 0.5*magSqr(U1));
-    volScalarField K2("K2", 0.5*magSqr(U2));
+    volScalarField K1("K" + phase1Name, 0.5*magSqr(U1));
+    volScalarField K2("K" + phase2Name, 0.5*magSqr(U2));
diff --git a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/createRASTurbulence.H b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/createRASTurbulence.H
index 8f05a1d1b373ccde772bfba4cee1b575574ae9f3..0a782ef51edf99f32eaf986325fd926db596e410 100644
--- a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/createRASTurbulence.H
+++ b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/createRASTurbulence.H
@@ -151,7 +151,7 @@
     (
         IOobject
         (
-            "nut2",
+            "nut" + phase2Name,
             runTime.timeName(),
             mesh,
             IOobject::NO_READ,
@@ -165,13 +165,13 @@
     (
         IOobject
         (
-            "nuEff1",
+            "nuEff" + phase1Name,
             runTime.timeName(),
             mesh,
             IOobject::NO_READ,
             IOobject::NO_WRITE
         ),
-        sqr(Ct)*nut2 + nu1
+        sqr(Ct)*nut2 + thermo1.mu()/rho1
     );
 
     Info<< "Calculating field nuEff2\n" << endl;
@@ -179,11 +179,11 @@
     (
         IOobject
         (
-            "nuEff2",
+            "nuEff" + phase2Name,
             runTime.timeName(),
             mesh,
             IOobject::NO_READ,
             IOobject::NO_WRITE
         ),
-        nut2 + nu2
+        nut2 + thermo2.mu()/rho2
     );
diff --git a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/interfacialCoeffs.H b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/interfacialCoeffs.H
index d6ccf90289ad6891248e665609408a68d4e7c766..d53bec5ea4f5ed2dc8393425bb4c4ddfa816c39d 100644
--- a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/interfacialCoeffs.H
+++ b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/interfacialCoeffs.H
@@ -38,12 +38,12 @@ volScalarField heatTransferCoeff
     volVectorField Ur(U1 - U2);
     volScalarField magUr(mag(Ur) + residualSlip);
 
-    if (dispersedPhase == "1")
+    if (dispersedPhase == phase1Name)
     {
         dragCoeff = drag1->K(magUr);
         heatTransferCoeff = heatTransfer1->K(magUr);
     }
-    else if (dispersedPhase == "2")
+    else if (dispersedPhase == phase2Name)
     {
         dragCoeff = drag2->K(magUr);
         heatTransferCoeff = heatTransfer2->K(magUr);
diff --git a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/interfacialModels/Make/options b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/interfacialModels/Make/options
index 92fe0b0d74ab657d08fb917834ecd7b4afa9783d..f031e058986a9c4ac3c37c83fafeaf2f592f28df 100644
--- a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/interfacialModels/Make/options
+++ b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/interfacialModels/Make/options
@@ -1,7 +1,9 @@
 EXE_INC = \
+    -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
     -I../phaseModel/lnInclude
 
 LIB_LIBS = \
-    -lcompressiblePhaseModel
-
+    -lcompressiblePhaseModel \
+    -lfluidThermophysicalModels \
+    -lspecie
diff --git a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/interfacialModels/heatTransferModels/RanzMarshall/RanzMarshall.C b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/interfacialModels/heatTransferModels/RanzMarshall/RanzMarshall.C
index 9bc6cadaf66495c673b2e7c69781aeabdaaf850d..95ade6ba05944fc306d26ce7c7ed1ffff4d97d2d 100644
--- a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/interfacialModels/heatTransferModels/RanzMarshall/RanzMarshall.C
+++ b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/interfacialModels/heatTransferModels/RanzMarshall/RanzMarshall.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -72,8 +72,7 @@ Foam::tmp<Foam::volScalarField> Foam::heatTransferModels::RanzMarshall::K
 ) const
 {
     volScalarField Re(max(Ur*phase1_.d()/phase2_.nu(), scalar(1.0e-3)));
-    dimensionedScalar Prb =
-        phase2_.rho()*phase2_.nu()*phase2_.Cp()/phase2_.kappa();
+    volScalarField Prb(phase2_.rho()*phase2_.nu()*phase2_.Cp()/phase2_.kappa());
     volScalarField Nu(scalar(2) + 0.6*sqrt(Re)*cbrt(Prb));
 
     return 6.0*phase2_.kappa()*Nu/sqr(phase1_.d());
diff --git a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/kineticTheoryModels/Make/options b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/kineticTheoryModels/Make/options
index 2fcce9913d9408482098065ad511e280e366c83f..7fdabf573092bd4219e75fa41044236e02b7d50e 100644
--- a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/kineticTheoryModels/Make/options
+++ b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/kineticTheoryModels/Make/options
@@ -1,5 +1,6 @@
 EXE_INC = \
     -I$(LIB_SRC)/foam/lnInclude \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
     -I../phaseModel/lnInclude \
     -I../interfacialModels/lnInclude
diff --git a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/kineticTheoryModels/kineticTheoryModel/kineticTheoryModel.C b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/kineticTheoryModels/kineticTheoryModel/kineticTheoryModel.C
index 7d0767914ea7ef3d5ef2fa6f11bd9b1d8ef4abda..18db3a70df7a5265075b96ee74e3bc87a073576d 100644
--- a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/kineticTheoryModels/kineticTheoryModel/kineticTheoryModel.C
+++ b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/kineticTheoryModels/kineticTheoryModel/kineticTheoryModel.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -45,8 +45,8 @@ Foam::kineticTheoryModel::kineticTheoryModel
     phi1_(phase1.phi()),
     draga_(draga),
 
-    rho1_(phase1.rho()),
-    nu1_(phase1.nu()),
+    rho1_(phase1.rho()[0]), //***HGW
+    nu1_(phase1.nu()()[0]), //***HGW
 
     kineticTheoryProperties_
     (
@@ -120,7 +120,7 @@ Foam::kineticTheoryModel::kineticTheoryModel
     (
         IOobject
         (
-            "mu1",
+            "mu" + phase1.name(),
             U1_.time().timeName(),
             U1_.mesh(),
             IOobject::NO_READ,
diff --git a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/pEqn.H b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/pEqn.H
index f7010b1c459cc8419b658797c7d7e4998ba1ab3a..f9a7807d00d46f51c83c8ee2f21679866b898d68 100644
--- a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/pEqn.H
+++ b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/pEqn.H
@@ -1,6 +1,6 @@
 {
-    rho1 = rho10 + psi1*p;
-    rho2 = rho20 + psi2*p;
+    rho1 = thermo1.rho();
+    rho2 = thermo2.rho();
 
     surfaceScalarField alpha1f(fvc::interpolate(alpha1));
     surfaceScalarField alpha2f(scalar(1) - alpha1f);
@@ -11,10 +11,10 @@
     surfaceScalarField rAlphaAU1f(fvc::interpolate(alpha1*rAU1));
     surfaceScalarField rAlphaAU2f(fvc::interpolate(alpha2*rAU2));
 
-    volVectorField HbyA1("HbyA1", U1);
+    volVectorField HbyA1("HbyA" + phase1Name, U1);
     HbyA1 = rAU1*U1Eqn.H();
 
-    volVectorField HbyA2("HbyA2", U2);
+    volVectorField HbyA2("HbyA" + phase2Name, U2);
     HbyA2 = rAU2*U2Eqn.H();
 
     mrfZones.absoluteFlux(phi1.oldTime());
@@ -38,14 +38,14 @@
 
     surfaceScalarField phiHbyA1
     (
-        "phiHbyA1",
+        "phiHbyA" + phase1Name,
         (fvc::interpolate(HbyA1) & mesh.Sf())
       + fvc::ddtPhiCorr(rAU1, alpha1, U1, phi1)
     );
 
     surfaceScalarField phiHbyA2
     (
-        "phiHbyA2",
+        "phiHbyA" + phase2Name,
         (fvc::interpolate(HbyA2) & mesh.Sf())
       + fvc::ddtPhiCorr(rAU2, alpha2, U2, phi2)
     );
@@ -96,8 +96,16 @@
     //}
     //else
     {
-        surfaceScalarField phid1("phid1", fvc::interpolate(psi1)*phi1);
-        surfaceScalarField phid2("phid2", fvc::interpolate(psi2)*phi2);
+        surfaceScalarField phid1
+        (
+            "phid" + phase1Name,
+            fvc::interpolate(psi1)*phi1
+        );
+        surfaceScalarField phid2
+        (
+            "phid" + phase2Name,
+            fvc::interpolate(psi2)*phi2
+        );
 
         pEqnComp1 =
             fvc::ddt(rho1) + psi1*correction(fvm::ddt(p))
@@ -174,13 +182,15 @@
 
     p = max(p, pMin);
 
-    rho1 = rho10 + psi1*p;
-    rho2 = rho20 + psi2*p;
+    thermo1.correct();
+    thermo2.correct();
+    rho1 = thermo1.rho();
+    rho2 = thermo2.rho();
 
     K1 = 0.5*magSqr(U1);
     K2 = 0.5*magSqr(U2);
 
-    //***HGW if (thermo.dpdt())
+    if (thermo1.dpdt())
     {
         dpdt = fvc::ddt(p);
     }
diff --git a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/phaseModel/Make/options b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/phaseModel/Make/options
index 0ec11392098bd862e78e2790e25cf1f49ec557fd..e441b0417bb82b021ee7e0ae470e030900e43a71 100644
--- a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/phaseModel/Make/options
+++ b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/phaseModel/Make/options
@@ -1,6 +1,9 @@
 EXE_INC = \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
     -I$(LIB_SRC)/transportModels/incompressible/lnInclude
 
 LIB_LIBS = \
-    -lincompressibleTransportModels
+    -lincompressibleTransportModels \
+    -lfluidThermophysicalModels \
+    -lspecie
diff --git a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/phaseModel/phaseModel/phaseModel.C b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/phaseModel/phaseModel/phaseModel.C
index d80758ac8c7d6cb754af8fe5aeaf9d2e7e2f7b78..64d59b436d5f2e85d2bd263227878101b3b51c76 100644
--- a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/phaseModel/phaseModel/phaseModel.C
+++ b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/phaseModel/phaseModel/phaseModel.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -37,35 +37,22 @@ Foam::phaseModel::phaseModel
     const word& phaseName
 )
 :
-    dict_
+    volScalarField
     (
-        transportProperties.subDict("phase" + phaseName)
+        IOobject
+        (
+            "alpha" + phaseName,
+            mesh.time().timeName(),
+            mesh,
+            IOobject::READ_IF_PRESENT,
+            IOobject::AUTO_WRITE
+        ),
+        mesh,
+        dimensionedScalar("alpha", dimless, 0)
     ),
     name_(phaseName),
-    nu_
-    (
-        "nu",
-        dimensionSet(0, 2, -1, 0, 0),
-        dict_.lookup("nu")
-    ),
-    kappa_
-    (
-        "kappa",
-        dimensionSet(1, 1, -3, -1, 0),
-        dict_.lookup("kappa")
-    ),
-    Cp_
-    (
-        "Cp",
-        dimensionSet(0, 2, -2, -1, 0),
-        dict_.lookup("Cp")
-    ),
-    rho_
-    (
-        "rho",
-        dimDensity,
-        dict_.lookup("rho")
-    ),
+    phaseDict_(transportProperties.subDict(phaseName)),
+    thermo_(rhoThermo::New(mesh, phaseName)),
     U_
     (
         IOobject
@@ -79,6 +66,8 @@ Foam::phaseModel::phaseModel
         mesh
     )
 {
+    thermo_->validate("phaseModel " + phaseName, "h", "e");
+
     const word phiName = "phi" + phaseName;
 
     IOobject phiHeader
@@ -147,7 +136,7 @@ Foam::phaseModel::phaseModel
 
     dPtr_ = diameterModel::New
     (
-        dict_,
+        phaseDict_,
         *this
     );
 }
diff --git a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/phaseModel/phaseModel/phaseModel.H b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/phaseModel/phaseModel/phaseModel.H
index 90b2d2184e0d3f971beeb959fe46864805073e89..18f5a47e46950dbec54b8a364b170aa5a3d5bef2 100644
--- a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/phaseModel/phaseModel/phaseModel.H
+++ b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/phaseModel/phaseModel/phaseModel.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -36,6 +36,7 @@ SourceFiles
 #include "dimensionedScalar.H"
 #include "volFields.H"
 #include "surfaceFields.H"
+#include "rhoThermo.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -50,25 +51,18 @@ class diameterModel;
 \*---------------------------------------------------------------------------*/
 
 class phaseModel
+:
+    public volScalarField
 {
     // Private data
 
-        dictionary dict_;
-
         //- Name of phase
         word name_;
 
-        //- Kinematic viscosity
-        dimensionedScalar nu_;
-
-        //- Thermal conductivity
-        dimensionedScalar kappa_;
-
-        //- Heat capacity
-        dimensionedScalar Cp_;
+        dictionary phaseDict_;
 
-        //- Density
-        dimensionedScalar rho_;
+        //- Thermophysical properties
+        autoPtr<rhoThermo> thermo_;
 
         //- Velocity
         volVectorField U_;
@@ -116,24 +110,34 @@ public:
 
         tmp<volScalarField> d() const;
 
-        const dimensionedScalar& nu() const
+        tmp<volScalarField> nu() const
+        {
+            return thermo_->mu()/thermo_->rho();
+        }
+
+        tmp<volScalarField> kappa() const
+        {
+            return thermo_->kappa();
+        }
+
+        tmp<volScalarField> Cp() const
         {
-            return nu_;
+            return thermo_->Cp();
         }
 
-        const dimensionedScalar& kappa() const
+        const volScalarField& rho() const
         {
-            return kappa_;
+            return thermo_->rho();
         }
 
-        const dimensionedScalar& Cp() const
+        const rhoThermo& thermo() const
         {
-            return Cp_;
+            return thermo_();
         }
 
-        const dimensionedScalar& rho() const
+        rhoThermo& thermo()
         {
-            return rho_;
+            return thermo_();
         }
 
         const volVectorField& U() const
diff --git a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/turbulenceModel/kEpsilon.H b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/turbulenceModel/kEpsilon.H
index 7f53411529a31a855f706d13d801f5624d1b7c45..bc9a07b0a81105cdbacc53336c054253844f2d2b 100644
--- a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/turbulenceModel/kEpsilon.H
+++ b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/turbulenceModel/kEpsilon.H
@@ -61,4 +61,4 @@ if (turbulence)
     #include "wallViscosity.H"
 }
 
-nuEff2 = nut2 + nu2;
+nuEff2 = nut2 + thermo2.mu()/rho2;
diff --git a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/turbulenceModel/wallFunctions.H b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/turbulenceModel/wallFunctions.H
index c91cce9a00e2df5ca912a75a57a02d6f813d72c4..d85181cba2e020e57fb2063d09d4af4f14b01807 100644
--- a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/turbulenceModel/wallFunctions.H
+++ b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/turbulenceModel/wallFunctions.H
@@ -4,7 +4,6 @@
     scalar Cmu25 = ::pow(Cmu.value(), 0.25);
     scalar Cmu75 = ::pow(Cmu.value(), 0.75);
     scalar kappa_ = kappa.value();
-    scalar nu2_ = nu2.value();
 
     const fvPatchList& patches = mesh.boundary();
 
@@ -30,6 +29,8 @@
     forAll(patches, patchi)
     {
         const fvPatch& currPatch = patches[patchi];
+        const scalarField& mu2_ = thermo2.mu().boundaryField()[patchi];
+        const scalarField& rho2_ = rho2.boundaryField()[patchi];
 
         if (isA<wallFvPatch>(currPatch))
         {
@@ -52,7 +53,8 @@
                     /(kappa_*y[patchi][facei]);
 
                 G[faceCelli] +=
-                    (nut2w[facei] + nu2_)*magFaceGradU[facei]
+                    (nut2w[facei] + mu2_[facei]/rho2_[facei])
+                   *magFaceGradU[facei]
                    *Cmu25*::sqrt(k[faceCelli])
                    /(kappa_*y[patchi][facei]);
             }
diff --git a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/turbulenceModel/wallViscosity.H b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/turbulenceModel/wallViscosity.H
index b153a360146deefed6b9b4a9c2a8a8ecb66075d9..9aa032149ca51f65d51d412c4dc13b7d61b9c46c 100644
--- a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/turbulenceModel/wallViscosity.H
+++ b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/turbulenceModel/wallViscosity.H
@@ -2,13 +2,14 @@
     scalar Cmu25 = ::pow(Cmu.value(), 0.25);
     scalar kappa_ = kappa.value();
     scalar E_ = E.value();
-    scalar nu2_ = nu2.value();
 
     const fvPatchList& patches = mesh.boundary();
 
     forAll(patches, patchi)
     {
         const fvPatch& currPatch = patches[patchi];
+        const scalarField& mu2_ = thermo2.mu().boundaryField()[patchi];
+        const scalarField& rho2_ = rho2.boundaryField()[patchi];
 
         if (isA<wallFvPatch>(currPatch))
         {
@@ -20,11 +21,14 @@
 
                 // calculate yPlus
                 scalar yPlus =
-                    Cmu25*y[patchi][facei]*::sqrt(k[faceCelli])/nu2_;
+                    Cmu25*y[patchi][facei]*::sqrt(k[faceCelli])
+                   /(mu2_[facei]/rho2_[facei]);
 
                 if (yPlus > 11.6)
                 {
-                    nutw[facei] = nu2_*(yPlus*kappa_/::log(E_*yPlus) - 1);
+                    nutw[facei] =
+                        (mu2_[facei]/rho2_[facei])
+                       *(yPlus*kappa_/::log(E_*yPlus) - 1);
                 }
                 else
                 {
diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/alphaEqn.H b/applications/solvers/multiphase/twoPhaseEulerFoam/alphaEqn.H
index 34fdc24afc532d56c0fd0ddf0c08479237d2cdb9..ddb923cf8726096119eac8c6b62e4b763877c642 100644
--- a/applications/solvers/multiphase/twoPhaseEulerFoam/alphaEqn.H
+++ b/applications/solvers/multiphase/twoPhaseEulerFoam/alphaEqn.H
@@ -13,7 +13,7 @@
         surfaceScalarField alpha1f(fvc::interpolate(alpha1));
         surfaceScalarField phipp(ppMagf*fvc::snGrad(alpha1)*mesh.magSf());
         phir += phipp;
-        phic += fvc::interpolate(alpha1)*phipp;
+        phic += alpha1f*phipp;
     }
 
     for (int acorr=0; acorr<nAlphaCorr; acorr++)
@@ -52,18 +52,23 @@
 
         if (g0.value() > 0)
         {
-            ppMagf = rAU1f*fvc::interpolate
-            (
-                (1.0/(rho1*(alpha1 + scalar(0.0001))))
-               *g0*min(exp(preAlphaExp*(alpha1 - alphaMax)), expMax)
-            );
+            surfaceScalarField alpha1f(fvc::interpolate(alpha1));
+
+            // ppMagf = rAU1f*fvc::interpolate
+            // (
+            //     (1.0/(rho1*(alpha1 + scalar(0.0001))))
+            //    *g0*min(exp(preAlphaExp*(alpha1 - alphaMax)), expMax)
+            // );
+            ppMagf =
+                rAU1f/(alpha1f + scalar(0.0001))
+               *(g0/rho1)*min(exp(preAlphaExp*(alpha1f - alphaMax)), expMax);
 
             fvScalarMatrix alpha1Eqn
             (
                 fvm::ddt(alpha1) - fvc::ddt(alpha1)
               - fvm::laplacian
                 (
-                    (fvc::interpolate(alpha1) + scalar(0.0001))*ppMagf,
+                    alpha1f*ppMagf,
                     alpha1,
                     "laplacian(alpha1PpMag,alpha1)"
                 )
diff --git a/src/OpenFOAM/matrices/lduMatrix/smoothers/DILU/DILUSmoother.H b/src/OpenFOAM/matrices/lduMatrix/smoothers/DILU/DILUSmoother.H
index 5aec38af8684fc145149df546d0e029e5691e969..843b386e77fd49bda9cec52ae59a69970c1da600 100644
--- a/src/OpenFOAM/matrices/lduMatrix/smoothers/DILU/DILUSmoother.H
+++ b/src/OpenFOAM/matrices/lduMatrix/smoothers/DILU/DILUSmoother.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -27,9 +27,6 @@ Class
 Description
     Simplified diagonal-based incomplete LU smoother for asymmetric matrices.
 
-    To improve efficiency, the residual is evaluated after every nSweeps
-    sweeps.
-
 SourceFiles
     DILUSmoother.C
 
diff --git a/src/finiteVolume/finiteVolume/fvc/fvcReconstruct.C b/src/finiteVolume/finiteVolume/fvc/fvcReconstruct.C
index ac3762db048d5c5476ed42b99e9608d5176ad65b..0329782315c4d4f4af6ef77d1c4da98602b6a6dd 100644
--- a/src/finiteVolume/finiteVolume/fvc/fvcReconstruct.C
+++ b/src/finiteVolume/finiteVolume/fvc/fvcReconstruct.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -56,6 +56,48 @@ reconstruct
 
     const fvMesh& mesh = ssf.mesh();
 
+    surfaceVectorField faceVols
+    (
+        mesh.Sf()/(mesh.magSf()*mesh.nonOrthDeltaCoeffs())
+    );
+
+    faceVols.internalField() *= (1.0 -  mesh.weights().internalField());
+    forAll(faceVols.boundaryField(), patchi)
+    {
+        if (faceVols.boundaryField()[patchi].coupled())
+        {
+            faceVols.boundaryField()[patchi] *=
+                (1.0 -  mesh.weights().boundaryField()[patchi]);
+        }
+    }
+
+    tmp<GeometricField<GradType, fvPatchField, volMesh> > treconField
+    (
+        new GeometricField<GradType, fvPatchField, volMesh>
+        (
+            IOobject
+            (
+                "volIntegrate("+ssf.name()+')',
+                ssf.instance(),
+                mesh,
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            inv(surfaceSum(mesh.Sf()*faceVols))&surfaceSum(faceVols*ssf),
+            zeroGradientFvPatchField<GradType>::typeName
+        )
+    );
+
+    treconField().correctBoundaryConditions();
+
+    return treconField;
+}
+/*
+{
+    typedef typename outerProduct<vector, Type>::type GradType;
+
+    const fvMesh& mesh = ssf.mesh();
+
     tmp<GeometricField<GradType, fvPatchField, volMesh> > treconField
     (
         new GeometricField<GradType, fvPatchField, volMesh>
@@ -78,6 +120,7 @@ reconstruct
 
     return treconField;
 }
+*/
 
 
 template<class Type>
diff --git a/tutorials/multiphase/compressibleTwoPhaseEulerFoam/bubbleColumn/0/T1 b/tutorials/multiphase/compressibleTwoPhaseEulerFoam/bubbleColumn/0/Tair
similarity index 95%
rename from tutorials/multiphase/compressibleTwoPhaseEulerFoam/bubbleColumn/0/T1
rename to tutorials/multiphase/compressibleTwoPhaseEulerFoam/bubbleColumn/0/Tair
index b07ba8e4080e6d222a80228aa49f6880ee0f895e..25bc08d671014e05362550e22ecd36921d254728 100644
--- a/tutorials/multiphase/compressibleTwoPhaseEulerFoam/bubbleColumn/0/T1
+++ b/tutorials/multiphase/compressibleTwoPhaseEulerFoam/bubbleColumn/0/Tair
@@ -10,7 +10,7 @@ FoamFile
     version     2.0;
     format      ascii;
     class       volScalarField;
-    object      T1;
+    object      Tair;
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -27,7 +27,7 @@ boundaryField
     outlet
     {
         type               inletOutlet;
-        phi                phi1;
+        phi                phiair;
         inletValue         $internalField;
         value              $internalField;
     }
diff --git a/tutorials/multiphase/compressibleTwoPhaseEulerFoam/bubbleColumn/0/T2 b/tutorials/multiphase/compressibleTwoPhaseEulerFoam/bubbleColumn/0/Twater
similarity index 95%
rename from tutorials/multiphase/compressibleTwoPhaseEulerFoam/bubbleColumn/0/T2
rename to tutorials/multiphase/compressibleTwoPhaseEulerFoam/bubbleColumn/0/Twater
index cc4307a53689ef43412a13c49df25ccc731623c6..f0f366c7dc27735b261b92c9eb7e2570dcaff136 100644
--- a/tutorials/multiphase/compressibleTwoPhaseEulerFoam/bubbleColumn/0/T2
+++ b/tutorials/multiphase/compressibleTwoPhaseEulerFoam/bubbleColumn/0/Twater
@@ -10,7 +10,7 @@ FoamFile
     version     2.0;
     format      ascii;
     class       volScalarField;
-    object      T2;
+    object      Twater;
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -27,7 +27,7 @@ boundaryField
     outlet
     {
         type               inletOutlet;
-        phi                phi2;
+        phi                phiwater;
         inletValue         uniform 300;
         value              $internalField;
     }
diff --git a/tutorials/multiphase/compressibleTwoPhaseEulerFoam/bubbleColumn/0/U1 b/tutorials/multiphase/compressibleTwoPhaseEulerFoam/bubbleColumn/0/Uair
similarity index 95%
rename from tutorials/multiphase/compressibleTwoPhaseEulerFoam/bubbleColumn/0/U1
rename to tutorials/multiphase/compressibleTwoPhaseEulerFoam/bubbleColumn/0/Uair
index 2cf5dec2bc7111419b65e3a82852df94d550df37..ac6020947e16d1537f7ed05c0878183567bac15c 100644
--- a/tutorials/multiphase/compressibleTwoPhaseEulerFoam/bubbleColumn/0/U1
+++ b/tutorials/multiphase/compressibleTwoPhaseEulerFoam/bubbleColumn/0/Uair
@@ -10,7 +10,7 @@ FoamFile
     version     2.0;
     format      binary;
     class       volVectorField;
-    object      U1;
+    object      Uair;
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -28,7 +28,7 @@ boundaryField
     outlet
     {
         type               pressureInletOutletVelocity;
-        phi                phi1;
+        phi                phiair;
         value              $internalField;
     }
     walls
diff --git a/tutorials/multiphase/compressibleTwoPhaseEulerFoam/bubbleColumn/0/U2 b/tutorials/multiphase/compressibleTwoPhaseEulerFoam/bubbleColumn/0/Uwater
similarity index 95%
rename from tutorials/multiphase/compressibleTwoPhaseEulerFoam/bubbleColumn/0/U2
rename to tutorials/multiphase/compressibleTwoPhaseEulerFoam/bubbleColumn/0/Uwater
index 68e9f354a6e5cab2f825bfd1e3b28d601d791e4d..22ed59a0ef3e78143b36a455a5c2d9f41f09df3a 100644
--- a/tutorials/multiphase/compressibleTwoPhaseEulerFoam/bubbleColumn/0/U2
+++ b/tutorials/multiphase/compressibleTwoPhaseEulerFoam/bubbleColumn/0/Uwater
@@ -10,7 +10,7 @@ FoamFile
     version     2.0;
     format      binary;
     class       volVectorField;
-    object      U2;
+    object      Uwater;
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -28,7 +28,7 @@ boundaryField
     outlet
     {
         type               pressureInletOutletVelocity;
-        phi                phi2;
+        phi                phiwater;
         value              $internalField;
     }
     walls
diff --git a/tutorials/multiphase/compressibleTwoPhaseEulerFoam/bubbleColumn/0/alpha1 b/tutorials/multiphase/compressibleTwoPhaseEulerFoam/bubbleColumn/0/alphaair
similarity index 99%
rename from tutorials/multiphase/compressibleTwoPhaseEulerFoam/bubbleColumn/0/alpha1
rename to tutorials/multiphase/compressibleTwoPhaseEulerFoam/bubbleColumn/0/alphaair
index 9a514b78074968c1df0cd3b7baa6ae041996ef8a..a2224f95c3c32f38ce0e3e0dd361b3fabc5bf60f 100644
--- a/tutorials/multiphase/compressibleTwoPhaseEulerFoam/bubbleColumn/0/alpha1
+++ b/tutorials/multiphase/compressibleTwoPhaseEulerFoam/bubbleColumn/0/alphaair
@@ -1908,7 +1908,7 @@ boundaryField
     outlet
     {
         type            inletOutlet;
-        phi             phi1;
+        phi             phiair;
         inletValue      uniform 1;
         value           uniform 1;
     }
diff --git a/tutorials/multiphase/compressibleTwoPhaseEulerFoam/bubbleColumn/0/alpha1.org b/tutorials/multiphase/compressibleTwoPhaseEulerFoam/bubbleColumn/0/alphaair.org
similarity index 95%
rename from tutorials/multiphase/compressibleTwoPhaseEulerFoam/bubbleColumn/0/alpha1.org
rename to tutorials/multiphase/compressibleTwoPhaseEulerFoam/bubbleColumn/0/alphaair.org
index 4b871e21b9da4e1ddd6159564e28dbf7f303b297..62088eca2323cb4b66162a3131679877f60bb945 100644
--- a/tutorials/multiphase/compressibleTwoPhaseEulerFoam/bubbleColumn/0/alpha1.org
+++ b/tutorials/multiphase/compressibleTwoPhaseEulerFoam/bubbleColumn/0/alphaair.org
@@ -10,7 +10,7 @@ FoamFile
     version     2.0;
     format      ascii;
     class       volScalarField;
-    object      alpha1;
+    object      alphaair;
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -28,7 +28,7 @@ boundaryField
     outlet
     {
         type            inletOutlet;
-        phi             phi1;
+        phi             phiair;
         inletValue      uniform 1;
         value           uniform 1;
     }
diff --git a/tutorials/multiphase/compressibleTwoPhaseEulerFoam/bubbleColumn/constant/interfacialProperties b/tutorials/multiphase/compressibleTwoPhaseEulerFoam/bubbleColumn/constant/interfacialProperties
index 03a3a667a3d8c1ee3e013766f338b7928053cabc..5964adcedc85ce4cfb8757b6dd2f652ba6cd9b88 100644
--- a/tutorials/multiphase/compressibleTwoPhaseEulerFoam/bubbleColumn/constant/interfacialProperties
+++ b/tutorials/multiphase/compressibleTwoPhaseEulerFoam/bubbleColumn/constant/interfacialProperties
@@ -15,13 +15,13 @@ FoamFile
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-dragModel1          SchillerNaumann;
-dragModel2          SchillerNaumann;
+dragModelair            SchillerNaumann;
+dragModelwater          SchillerNaumann;
 
-heatTransferModel1  RanzMarshall;
-heatTransferModel2  RanzMarshall;
+heatTransferModelair    RanzMarshall;
+heatTransferModelwater  RanzMarshall;
 
-dispersedPhase      both;
+dispersedPhase          both;
 
 residualPhaseFraction   1e-3;
 residualSlip            1e-2;
diff --git a/tutorials/multiphase/compressibleTwoPhaseEulerFoam/bubbleColumn/constant/polyMesh/boundary b/tutorials/multiphase/compressibleTwoPhaseEulerFoam/bubbleColumn/constant/polyMesh/boundary
index 56e0a545c15839416fcd3beca616ce0f609c5879..bf47f69643c9925d3a1ef19c6b4ddc67cf604e0a 100644
--- a/tutorials/multiphase/compressibleTwoPhaseEulerFoam/bubbleColumn/constant/polyMesh/boundary
+++ b/tutorials/multiphase/compressibleTwoPhaseEulerFoam/bubbleColumn/constant/polyMesh/boundary
@@ -38,6 +38,7 @@ FoamFile
     defaultFaces
     {
         type            empty;
+        inGroups        1(empty);
         nFaces          3750;
         startFace       3850;
     }
diff --git a/tutorials/multiphase/compressibleTwoPhaseEulerFoam/bubbleColumn/constant/thermophysicalPropertiesair b/tutorials/multiphase/compressibleTwoPhaseEulerFoam/bubbleColumn/constant/thermophysicalPropertiesair
new file mode 100644
index 0000000000000000000000000000000000000000..11c033af59bbfbf08d2b11244e50f0fd00c1542f
--- /dev/null
+++ b/tutorials/multiphase/compressibleTwoPhaseEulerFoam/bubbleColumn/constant/thermophysicalPropertiesair
@@ -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;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+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/compressibleTwoPhaseEulerFoam/bubbleColumn/constant/thermophysicalPropertieswater b/tutorials/multiphase/compressibleTwoPhaseEulerFoam/bubbleColumn/constant/thermophysicalPropertieswater
new file mode 100644
index 0000000000000000000000000000000000000000..672b24a98b365f170f3f2567723294befd79a06e
--- /dev/null
+++ b/tutorials/multiphase/compressibleTwoPhaseEulerFoam/bubbleColumn/constant/thermophysicalPropertieswater
@@ -0,0 +1,53 @@
+/*--------------------------------*- 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;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+thermoType
+{
+    type            heRhoThermo;
+    mixture         pureMixture;
+    transport       const;
+    thermo          hConst;
+    equationOfState perfectFluid;
+    specie          specie;
+    energy          sensibleEnthalpy;
+}
+
+mixture
+{
+    specie
+    {
+        nMoles      1;
+        molWeight   28.9;
+    }
+    equationOfState
+    {
+        rho0        1027;
+    }
+    thermodynamics
+    {
+        Cp          4195;
+        Hf          0;
+    }
+    transport
+    {
+        mu          3.645e-4;
+        Pr          2.289;
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/compressibleTwoPhaseEulerFoam/bubbleColumn/constant/transportProperties b/tutorials/multiphase/compressibleTwoPhaseEulerFoam/bubbleColumn/constant/transportProperties
index d469dc706f932ef77a08bd0aa00f6cd9ec867c91..408d5f97d66a726605f89ae7c8b5a3e18711f2a6 100644
--- a/tutorials/multiphase/compressibleTwoPhaseEulerFoam/bubbleColumn/constant/transportProperties
+++ b/tutorials/multiphase/compressibleTwoPhaseEulerFoam/bubbleColumn/constant/transportProperties
@@ -15,16 +15,10 @@ FoamFile
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-phase1
-{
-    rho             1;
-    rho0            0;
-    R               287;
-    Cp              1007;
-
-    nu              1.589e-05;
-    kappa           2.63e-2;
+phases (air water);
 
+air
+{
     diameterModel   isothermal;
     isothermalCoeffs
     {
@@ -33,15 +27,9 @@ phase1
     }
 }
 
-phase2
+water
 {
-    rho             1027;
-    rho0            1027;
-    R               1e10;
-    Cp              4195;
-
-    nu              3.55e-7;
-    kappa           0.668;
+    //R               1e10;
 
     diameterModel constant;
     constantCoeffs
diff --git a/tutorials/multiphase/compressibleTwoPhaseEulerFoam/bubbleColumn/system/controlDict b/tutorials/multiphase/compressibleTwoPhaseEulerFoam/bubbleColumn/system/controlDict
index f28ade7e84c4a5d839fed48c74b21ffbdca6b371..24d66f91c7d44dfeb7b24b508c2593511edb6ed9 100644
--- a/tutorials/multiphase/compressibleTwoPhaseEulerFoam/bubbleColumn/system/controlDict
+++ b/tutorials/multiphase/compressibleTwoPhaseEulerFoam/bubbleColumn/system/controlDict
@@ -60,21 +60,21 @@ functions
         outputControl   outputTime;
         fields
         (
-            U1
+            Uair
             {
                  mean        on;
                  prime2Mean  off;
                  base        time;
             }
 
-            U2
+            Uwater
             {
                  mean        on;
                  prime2Mean  off;
                  base        time;
             }
 
-            alpha1
+            alphaair
             {
                  mean        on;
                  prime2Mean  off;
diff --git a/tutorials/multiphase/compressibleTwoPhaseEulerFoam/bubbleColumn/system/fvSchemes b/tutorials/multiphase/compressibleTwoPhaseEulerFoam/bubbleColumn/system/fvSchemes
index 34b14edfcce7676bc1e185e87f785fc2ab086c56..d53ef42583a230181a0d5bc0a16d73acdc017ce2 100644
--- a/tutorials/multiphase/compressibleTwoPhaseEulerFoam/bubbleColumn/system/fvSchemes
+++ b/tutorials/multiphase/compressibleTwoPhaseEulerFoam/bubbleColumn/system/fvSchemes
@@ -29,23 +29,23 @@ divSchemes
 {
     default             none;
 
-    div(phi,alpha1)     Gauss limitedLinear01 1;
-    div(phir,alpha1)    Gauss limitedLinear01 1;
-    div(alphaPhi1,U1)   Gauss limitedLinearV 1;
-    div(alphaPhi2,U2)   Gauss limitedLinearV 1;
-    div(phi1,U1)        Gauss limitedLinearV 1;
-    div(phi2,U2)        Gauss limitedLinearV 1;
-    div((alpha1*Rc1))   Gauss linear;
-    div((alpha2*Rc2))   Gauss linear;
-    div(alphaPhi1,T1)   Gauss limitedLinear 1;
-    div(alphaPhi2,T2)   Gauss limitedLinear 1;
-    div(alphaPhi2,k)    Gauss limitedLinear 1;
-    div(alphaPhi2,epsilon) Gauss limitedLinear 1;
+    div(phi,alpha)     Gauss limitedLinear01 1;
+    div(phir,alpha)    Gauss limitedLinear01 1;
+    div(alphaPhi,Uair)   Gauss limitedLinearV 1;
+    div(alphaPhi,Uwater)   Gauss limitedLinearV 1;
+    div(phiair,Uair)        Gauss limitedLinearV 1;
+    div(phiwater,Uwater)        Gauss limitedLinearV 1;
+    div((alphaair*Rc))   Gauss linear;
+    div((alphawater*Rc))   Gauss linear;
+    div(alphaPhi,hair)   Gauss limitedLinear 1;
+    div(alphaPhi,hwater)   Gauss limitedLinear 1;
+    div(alphaPhiwater,k)    Gauss limitedLinear 1;
+    div(alphaPhiwater,epsilon) Gauss limitedLinear 1;
     div(phi,Theta)      Gauss limitedLinear 1;
-    div(phid1,p)        Gauss upwind;
-    div(phid2,p)        Gauss upwind;
-    div(phi1,K1)        Gauss limitedLinear 1;
-    div(phi2,K2)        Gauss limitedLinear 1;
+    div(phidair,p)        Gauss upwind;
+    div(phidwater,p)        Gauss upwind;
+    div(phiair,Kair)        Gauss limitedLinear 1;
+    div(phiwater,Kwater)        Gauss limitedLinear 1;
 }
 
 laplacianSchemes
@@ -67,7 +67,7 @@ fluxRequired
 {
     default         no;
     p               ;
-    alpha1          ;
+    alphaair        ;
 }
 
 
diff --git a/tutorials/multiphase/compressibleTwoPhaseEulerFoam/bubbleColumn/system/fvSolution b/tutorials/multiphase/compressibleTwoPhaseEulerFoam/bubbleColumn/system/fvSolution
index d09194f954dde0ac71cc6201ac1cb1f475455c51..96e43b3c1d253411d40e8e71b59a2c071ad2f124 100644
--- a/tutorials/multiphase/compressibleTwoPhaseEulerFoam/bubbleColumn/system/fvSolution
+++ b/tutorials/multiphase/compressibleTwoPhaseEulerFoam/bubbleColumn/system/fvSolution
@@ -47,10 +47,10 @@ solvers
         relTol          0;
     }
 
-    "T.*"
+    "h.*"
     {
-        solver          PBiCG;
-        preconditioner  DILU;
+        solver          PCG; //PBiCG;
+        preconditioner  DIC; //DILU;
         tolerance       1e-8;
         relTol          0;
     }
@@ -105,7 +105,7 @@ relaxationFactors
     equations
     {
         "U.*"           1;
-        "T.*"           1;
+        "h.*"           1;
         "alpha.*"       1;
         "Theta.*"       1;
         "k.*"           1;
diff --git a/tutorials/multiphase/compressibleTwoPhaseEulerFoam/bubbleColumn/system/setFieldsDict b/tutorials/multiphase/compressibleTwoPhaseEulerFoam/bubbleColumn/system/setFieldsDict
index 5b6d65bbecb7870ba8626ebd335c4a8f8d75d8d9..85996cf966762c6ca3a7c37a1eaa8ae462ecdb19 100644
--- a/tutorials/multiphase/compressibleTwoPhaseEulerFoam/bubbleColumn/system/setFieldsDict
+++ b/tutorials/multiphase/compressibleTwoPhaseEulerFoam/bubbleColumn/system/setFieldsDict
@@ -27,7 +27,7 @@ regions
         box (0 0 -0.1) (0.15 0.701 0.1);
         fieldValues
         (
-            volScalarFieldValue alpha1 0
+            volScalarFieldValue alphaair 0
         );
     }
 );