Merge branch 'master' of ssh://opencfd:8007/home/dm4/OpenFOAM/OpenFOAM-dev

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

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)

applications/solvers/multiphase/twoPhaseEulerFoam/alphaEqn.H

@@ -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();

applications/solvers/multiphase/twoPhaseEulerFoam/createFields.H

@@ -10,19 +10,13 @@
 
     volVectorField& U1 = phase1.U();
     surfaceScalarField& phi1 = phase1.phi();
-    surfaceScalarField alphaPhi1
+    surfaceScalarField& alphaPhi1 = phase1.phiAlpha();
-    (
-        IOobject::groupName("alphaPhi", phase1.name()),
-        fvc::interpolate(alpha1)*phi1
-    );
 
     volVectorField& U2 = phase2.U();
     surfaceScalarField& phi2 = phase2.phi();
-    surfaceScalarField alphaPhi2
+    surfaceScalarField& alphaPhi2 = phase2.phiAlpha();
-    (
+
-        IOobject::groupName("alphaPhi", phase2.name()),
+    surfaceScalarField& phi = fluid.phi();
-        fvc::interpolate(alpha2)*phi2
-    );
 
     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
-        )
-    );

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

applications/solvers/multiphase/twoPhaseEulerFoam/readTwoPhaseEulerFoamControls.H

@@ -1,7 +0,0 @@
-    #include "readTimeControls.H"
-    #include "alphaControls.H"
-
-    Switch implicitPhasePressure
-    (
-        alphaControls.lookupOrDefault<Switch>("implicitPhasePressure", false)
-    );

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();
+                fluid.correctTurbulence();
-                turbulence2->correct();
             }
         }
 

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

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 \

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;
+}
+
+
+// ************************************************************************* //

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
+
+// ************************************************************************* //

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();
+}
+
+
+// ************************************************************************* //

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
+
+// ************************************************************************* //

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)
+        )
+    );
+}
+
+
+// ************************************************************************* //

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
+
+// ************************************************************************* //

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;
+}
+
+
+// ************************************************************************* //

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
+
+// ************************************************************************* //

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;
+}
+
+
+// ************************************************************************* //

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
+
+// ************************************************************************* //

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();
+}
+
+
+// ************************************************************************* //

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
+
+// ************************************************************************* //

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),
+    diameterModel(diameterProperties, phase),
-    d_("d", dimLength, dict.lookup("d"))
+    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;
+}
+
+
 // ************************************************************************* //

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);
 };
 
 

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;
+}
+
+
 // ************************************************************************* //

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;
 };
 
 

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),
+    diameterModel(diameterProperties, phase),
-    d0_("d0", dimLength, dict.lookup("d0")),
+    d0_("d0", dimLength, diameterProperties_.lookup("d0")),
-    p0_("p0", dimPressure, dict.lookup("p0"))
+    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;
+}
+
+
 // ************************************************************************* //

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);
 };
 
 

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_);
+}
+
 
 // ************************************************************************* //

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 name_;
-        }
 
+        //- 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 thermo_();
+        }
+
+        //- 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 the velocity
-        {
-            return thermo_();
-        }
-
-        rhoThermo& thermo()
-        {
-            return thermo_();
-        }
-
         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
+        //- Return the volumetric flux of the phase
-        void correct()
+        const surfaceScalarField& phiAlpha() const
-        {}
+        {
+            return phiAlpha_;
+        }
 
-        //- Dummy read
+        //- Return non-const access to the volumetric flux of the phase
-        bool read()
+        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 for transportModel
+        virtual bool read()
         {
             return true;
         }

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 &= phase1_.read(*this);
-        readOK &= phase2_.read();
+        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

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
+        //- Return the surface tension coefficient
-        tmp<surfaceScalarField> phi() const;
+        dimensionedScalar sigma() const
+        {
+            return sigma_;
+        }
 
         //- Return the virtual-mass coefficient
         dimensionedScalar Cvm() const
@@ -208,9 +243,14 @@ public:
             return Cl_;
         }
 
-        //- Dummy correct
+        //- Solve for the two-phase-fractions
-        void correct()
+        void solve();
-        {}
+
+        //- Correct two-phase properties other than turbulence
+        void correct();
+
+        //- Correct two-phase turbulence
+        void correctTurbulence();
 
         //- Read base phaseProperties dictionary
         bool read();

applications/utilities/mesh/manipulation/checkMesh/checkTopology.C

@@ -1,3 +1,28 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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 "checkTopology.H"
 #include "polyMesh.H"
 #include "Time.H"
@@ -9,6 +34,8 @@
 #include "emptyPolyPatch.H"
 #include "processorPolyPatch.H"
 
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
 Foam::label Foam::checkTopology
 (
     const polyMesh& mesh,
@@ -287,6 +314,37 @@ Foam::label Foam::checkTopology
                 rs
             );
             ctr.write();
+
+
+            // write cellSet for each region
+            PtrList<cellSet> cellRegions(rs.nRegions());
+            for (label i = 0; i < rs.nRegions(); i++)
+            {
+                cellRegions.set
+                (
+                    i,
+                    new cellSet
+                    (
+                        mesh,
+                        "region" + Foam::name(i),
+                        mesh.nCells()/100
+                    )
+                );
+            }
+
+            forAll(rs, i)
+            {
+                cellRegions[rs[i]].insert(i);
+            }
+
+            for (label i = 0; i < rs.nRegions(); i++)
+            {
+                Info<< "  <<Writing region " << i << " with "
+                    << returnReduce(cellRegions[i].size(), sumOp<scalar>())
+                    << " cells to cellSet " << cellRegions[i].name() << endl;
+
+                cellRegions[i].write();
+            }
         }
     }
 
@@ -437,3 +495,6 @@ Foam::label Foam::checkTopology
 
     return noFailedChecks;
 }
+
+
+// ************************************************************************* //

src/OpenFOAM/fields/GeometricFields/GeometricField/GeometricField.C

@@ -934,8 +934,8 @@ void Foam::GeometricField<Type, PatchField, GeoMesh>::writeMinMax
 ) const
 {
     os  << "min/max(" << this->name() << ") = "
-        << min(*this).value() << ", "
+        << Foam::min(*this).value() << ", "
-        << max(*this).value()
+        << Foam::max(*this).value()
         << endl;
 }
 

src/dynamicMesh/layerAdditionRemoval/addCellLayer.C

@@ -161,6 +161,9 @@ void Foam::layerAdditionRemoval::addCellLayer
 
     forAll(mf, faceI)
     {
+        label cellI = mc[faceI];
+        label zoneI =  mesh.cellZones().whichZone(cellI);
+
         addedCells[faceI] =
             ref.setAction
             (
@@ -170,7 +173,7 @@ void Foam::layerAdditionRemoval::addCellLayer
                     -1,           // master edge
                     mf[faceI],    // master face
                     -1,           // master cell
-                    -1            // zone for cell
+                    zoneI         // zone for cell
                 )
             );
     }

src/finiteVolume/interpolation/surfaceInterpolation/limitedSchemes/LimitedScheme/LimitedScheme.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
@@ -28,33 +28,17 @@ License
 #include "fvcGrad.H"
 #include "coupledFvPatchFields.H"
 
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+// * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * * //
 
 template<class Type, class Limiter, template<class> class LimitFunc>
-Foam::tmp<Foam::surfaceScalarField>
+void Foam::LimitedScheme<Type, Limiter, LimitFunc>::calcLimiter
-Foam::LimitedScheme<Type, Limiter, LimitFunc>::limiter
 (
-    const GeometricField<Type, fvPatchField, volMesh>& phi
+    const GeometricField<Type, fvPatchField, volMesh>& phi,
+    surfaceScalarField& limiterField
 ) const
 {
     const fvMesh& mesh = this->mesh();
 
-    tmp<surfaceScalarField> tLimiter
-    (
-        new surfaceScalarField
-        (
-            IOobject
-            (
-                type() + "Limiter(" + phi.name() + ')',
-                mesh.time().timeName(),
-                mesh
-            ),
-            mesh,
-            dimless
-        )
-    );
-    surfaceScalarField& lim = tLimiter();
-
     tmp<GeometricField<typename Limiter::phiType, fvPatchField, volMesh> >
         tlPhi = LimitFunc<Type>()(phi);
 
@@ -73,7 +57,7 @@ Foam::LimitedScheme<Type, Limiter, LimitFunc>::limiter
 
     const vectorField& C = mesh.C();
 
-    scalarField& pLim = lim.internalField();
+    scalarField& pLim = limiterField.internalField();
 
     forAll(pLim, face)
     {
@@ -92,7 +76,8 @@ Foam::LimitedScheme<Type, Limiter, LimitFunc>::limiter
         );
     }
 
-    surfaceScalarField::GeometricBoundaryField& bLim = lim.boundaryField();
+    surfaceScalarField::GeometricBoundaryField& bLim =
+        limiterField.boundaryField();
 
     forAll(bLim, patchi)
     {
@@ -143,8 +128,80 @@ Foam::LimitedScheme<Type, Limiter, LimitFunc>::limiter
             pLim = 1.0;
         }
     }
+}
 
-    return tLimiter;
+
+// * * * * * * * * * * * * Public Member Functions  * * * * * * * * * * * * //
+
+template<class Type, class Limiter, template<class> class LimitFunc>
+Foam::tmp<Foam::surfaceScalarField>
+Foam::LimitedScheme<Type, Limiter, LimitFunc>::limiter
+(
+    const GeometricField<Type, fvPatchField, volMesh>& phi
+) const
+{
+    const fvMesh& mesh = this->mesh();
+
+    const word limiterFieldName(type() + "Limiter(" + phi.name() + ')');
+
+    if (this->mesh().cache("limiter"))
+    {
+        if (!mesh.foundObject<surfaceScalarField>(limiterFieldName))
+        {
+            surfaceScalarField* limiterField
+            (
+                new surfaceScalarField
+                (
+                    IOobject
+                    (
+                        limiterFieldName,
+                        mesh.time().timeName(),
+                        mesh,
+                        IOobject::NO_READ,
+                        IOobject::NO_WRITE
+                    ),
+                    mesh,
+                    dimless
+                )
+            );
+
+            mesh.objectRegistry::store(limiterField);
+        }
+
+        surfaceScalarField& limiterField =
+            const_cast<surfaceScalarField&>
+            (
+                mesh.lookupObject<surfaceScalarField>
+                (
+                    limiterFieldName
+                )
+            );
+
+        calcLimiter(phi, limiterField);
+
+        return limiterField;
+    }
+    else
+    {
+        tmp<surfaceScalarField> tlimiterField
+        (
+            new surfaceScalarField
+            (
+                IOobject
+                (
+                    limiterFieldName,
+                    mesh.time().timeName(),
+                    mesh
+                ),
+                mesh,
+                dimless
+            )
+        );
+
+        calcLimiter(phi, tlimiterField());
+
+        return tlimiterField;
+    }
 }
 
 

src/finiteVolume/interpolation/surfaceInterpolation/limitedSchemes/LimitedScheme/LimitedScheme.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
@@ -57,7 +57,7 @@ namespace Foam
 {
 
 /*---------------------------------------------------------------------------*\
-                           Class LimitedScheme Declaration
+                        Class LimitedScheme Declaration
 \*---------------------------------------------------------------------------*/
 
 template<class Type, class Limiter, template<class> class LimitFunc>
@@ -68,6 +68,13 @@ class LimitedScheme
 {
     // Private Member Functions
 
+        //- Calculate the limiter
+        void calcLimiter
+        (
+            const GeometricField<Type, fvPatchField, volMesh>& phi,
+            surfaceScalarField& limiterField
+        ) const;
+
         //- Disallow default bitwise copy construct
         LimitedScheme(const LimitedScheme&);
 

src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/InjectionModel/InjectionModel.C

@@ -550,7 +550,7 @@ void Foam::InjectionModel<CloudType>::inject(TrackData& td)
                 if (cellI > -1)
                 {
                     // Lagrangian timestep
-                    scalar dt = time - timeInj;
+                    const scalar dt = time - timeInj;
 
                     // Apply corrections to position for 2-D cases
                     meshTools::constrainToMeshCentre(mesh, pos);
@@ -586,8 +586,11 @@ void Foam::InjectionModel<CloudType>::inject(TrackData& td)
                             pPtr->rho()
                         );
 
+                    const scalar mParcel0 = pPtr->nParticle()*pPtr->mass();
+
                     if (!pPtr->move(td, dt))
                     {
+                        massAdded += mParcel0;
                         delete pPtr;
                     }
                     else
@@ -595,7 +598,7 @@ void Foam::InjectionModel<CloudType>::inject(TrackData& td)
                         if (pPtr->nParticle() >= 1.0)
                         {
                             td.cloud().addParticle(pPtr);
-                            massAdded += pPtr->nParticle()*pPtr->mass();
+                            massAdded += mParcel0;
                             parcelsAdded++;
                         }
                         else

src/postProcessing/functionObjects/IO/writeRegisteredObject/writeRegisteredObject.C

@@ -91,25 +91,14 @@ void Foam::writeRegisteredObject::write()
 {
     Info<< type() << " " << name_ << " output:" << nl;
 
+    DynamicList<word> allNames(obr_.toc().size());
     forAll(objectNames_, i)
     {
-        if (obr_.foundObject<regIOobject>(objectNames_[i]))
+        wordList names(obr_.names<regIOobject>(objectNames_[i]));
+
+        if (names.size())
         {
-            regIOobject& obj =
+            allNames.append(names);
-                const_cast<regIOobject&>
-                (
-                    obr_.lookupObject<regIOobject>(objectNames_[i])
-                );
-
-            if (exclusiveWriting_)
-            {
-                // Switch off automatic writing to prevent double write
-                obj.writeOpt() = IOobject::NO_WRITE;
-            }
-
-            Info<< "    writing object " << obj.name() << nl << endl;
-
-            obj.write();
         }
         else
         {
@@ -119,6 +108,25 @@ void Foam::writeRegisteredObject::write()
                 << endl;
         }
     }
+
+    forAll(allNames, i)
+    {
+        regIOobject& obj =
+            const_cast<regIOobject&>
+            (
+                obr_.lookupObject<regIOobject>(allNames[i])
+            );
+
+        if (exclusiveWriting_)
+        {
+            // Switch off automatic writing to prevent double write
+            obj.writeOpt() = IOobject::NO_WRITE;
+        }
+
+        Info<< "    writing object " << obj.name() << nl << endl;
+
+        obj.write();
+    }
 }
 
 

src/postProcessing/functionObjects/IO/writeRegisteredObject/writeRegisteredObject.H

@@ -74,7 +74,7 @@ SourceFiles
 #ifndef writeRegisteredObject_H
 #define writeRegisteredObject_H
 
-#include "wordList.H"
+#include "wordReList.H"
 #include "runTimeSelectionTables.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -110,7 +110,7 @@ protected:
         // Read from dictionary
 
             //- Names of objects to control
-            wordList objectNames_;
+            wordReList objectNames_;
 
 
     // Private Member Functions

src/regionModels/surfaceFilmModels/Make/files

@@ -14,7 +14,6 @@ $(KINEMATICMODELS)/force/force/force.C
 $(KINEMATICMODELS)/force/force/forceNew.C
 $(KINEMATICMODELS)/force/forceList/forceList.C
 $(KINEMATICMODELS)/force/contactAngleForce/contactAngleForce.C
-$(KINEMATICMODELS)/force/surfaceShearForce/surfaceShearForce.C
 $(KINEMATICMODELS)/force/thermocapillaryForce/thermocapillaryForce.C
 
 $(KINEMATICMODELS)/injectionModel/injectionModel/injectionModel.C
@@ -24,6 +23,11 @@ $(KINEMATICMODELS)/injectionModel/drippingInjection/drippingInjection.C
 $(KINEMATICMODELS)/injectionModel/removeInjection/removeInjection.C
 $(KINEMATICMODELS)/injectionModel/curvatureSeparation/curvatureSeparation.C
 
+$(KINEMATICMODELS)/filmThermoModel/filmThermoModel/filmThermoModel.C
+$(KINEMATICMODELS)/filmThermoModel/filmThermoModel/filmThermoModelNew.C
+$(KINEMATICMODELS)/filmThermoModel/constantFilmThermo/constantFilmThermo.C
+$(KINEMATICMODELS)/filmThermoModel/liquidFilmThermo/liquidFilmThermo.C
+
 $(KINEMATICMODELS)/filmTurbulenceModel/filmTurbulenceModel/filmTurbulenceModel.C
 $(KINEMATICMODELS)/filmTurbulenceModel/filmTurbulenceModel/filmTurbulenceModelNew.C
 $(KINEMATICMODELS)/filmTurbulenceModel/laminar/laminar.C
@@ -33,6 +37,7 @@ $(THERMOMODELS)/phaseChangeModel/phaseChangeModel/phaseChangeModel.C
 $(THERMOMODELS)/phaseChangeModel/phaseChangeModel/phaseChangeModelNew.C
 $(THERMOMODELS)/phaseChangeModel/noPhaseChange/noPhaseChange.C
 $(THERMOMODELS)/phaseChangeModel/standardPhaseChange/standardPhaseChange.C
+$(THERMOMODELS)/phaseChangeModel/solidification/solidification.C
 
 $(THERMOMODELS)/heatTransferModel/heatTransferModel/heatTransferModel.C
 $(THERMOMODELS)/heatTransferModel/heatTransferModel/heatTransferModelNew.C
@@ -46,6 +51,12 @@ $(THERMOMODELS)/filmRadiationModel/constantRadiation/constantRadiation.C
 $(THERMOMODELS)/filmRadiationModel/primaryRadiation/primaryRadiation.C
 $(THERMOMODELS)/filmRadiationModel/standardRadiation/standardRadiation.C
 
+$(THERMOMODELS)/filmViscosityModel/filmViscosityModel/filmViscosityModel.C
+$(THERMOMODELS)/filmViscosityModel/filmViscosityModel/filmViscosityModelNew.C
+$(THERMOMODELS)/filmViscosityModel/constantViscosity/constantViscosity.C
+$(THERMOMODELS)/filmViscosityModel/liquidViscosity/liquidViscosity.C
+$(THERMOMODELS)/filmViscosityModel/thixotropicViscosity/thixotropicViscosity.C
+
 
 /* Boundary conditions */
 PATCHFIELDS=derivedFvPatchFields

src/regionModels/surfaceFilmModels/kinematicSingleLayer/kinematicSingleLayer.C

@@ -33,6 +33,7 @@ License
 #include "addToRunTimeSelectionTable.H"
 #include "mappedWallPolyPatch.H"
 #include "mapDistribute.H"
+#include "filmThermoModel.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -72,22 +73,9 @@ bool kinematicSingleLayer::read()
 
 void kinematicSingleLayer::correctThermoFields()
 {
-    if (thermoModel_ == tmConstant)
+    rho_ == filmThermo_->rho();
-    {
+    mu_ == filmThermo_->mu();
-        const dictionary& constDict(coeffs_.subDict("constantThermoCoeffs"));
+    sigma_ == filmThermo_->sigma();
-        rho_ == dimensionedScalar(constDict.lookup("rho0"));
-        mu_ == dimensionedScalar(constDict.lookup("mu0"));
-        sigma_ == dimensionedScalar(constDict.lookup("sigma0"));
-    }
-    else
-    {
-        FatalErrorIn
-        (
-            "void Foam::surfaceFilmModels::kinematicSingleLayer::"
-            "correctThermo()"
-        )   << "Kinematic surface film must use "
-            << thermoModelTypeNames_[tmConstant] << "thermodynamics" << endl;
-    }
 }
 
 
@@ -195,7 +183,7 @@ tmp<volScalarField> kinematicSingleLayer::pp()
 
 void kinematicSingleLayer::correctAlpha()
 {
-    alpha_ == pos(delta_ - dimensionedScalar("SMALL", dimLength, SMALL));
+    alpha_ == pos(delta_ - deltaSmall_);
 }
 
 
@@ -450,6 +438,8 @@ kinematicSingleLayer::kinematicSingleLayer
 
     cumulativeContErr_(0.0),
 
+    deltaSmall_("deltaSmall", dimLength, SMALL),
+
     rho_
     (
         IOobject
@@ -772,6 +762,8 @@ kinematicSingleLayer::kinematicSingleLayer
         this->mappedFieldAndInternalPatchTypes<scalar>()
     ),
 
+    filmThermo_(filmThermoModel::New(*this, coeffs_)),
+
     availableMass_(regionMesh().nCells(), 0.0),
 
     injection_(*this, coeffs_),
@@ -786,6 +778,8 @@ kinematicSingleLayer::kinematicSingleLayer
     {
         transferPrimaryRegionThermoFields();
 
+        correctAlpha();
+
         correctThermoFields();
 
         deltaRho_ == delta_*rho_;
@@ -944,6 +938,12 @@ const volVectorField& kinematicSingleLayer::Uw() const
 }
 
 
+const volScalarField& kinematicSingleLayer::deltaRho() const
+{
+    return deltaRho_;
+}
+
+
 const surfaceScalarField& kinematicSingleLayer::phi() const
 {
     return phi_;

src/regionModels/surfaceFilmModels/kinematicSingleLayer/kinematicSingleLayer.H

@@ -54,6 +54,9 @@ namespace regionModels
 namespace surfaceFilmModels
 {
 
+// Forward class declarations
+class filmThermoModel;
+
 /*---------------------------------------------------------------------------*\
                    Class kinematicSingleLayer Declaration
 \*---------------------------------------------------------------------------*/
@@ -94,6 +97,9 @@ protected:
             //- Cumulative continuity error
             scalar cumulativeContErr_;
 
+            //- Small delta
+            const dimensionedScalar deltaSmall_;
+
 
         // Thermo properties
 
@@ -192,6 +198,9 @@ protected:
 
         // Sub-models
 
+            //- Film thermo model
+            autoPtr<filmThermoModel> filmThermo_;
+
             //- Available mass for transfer via sub-models
             scalarField availableMass_;
 
@@ -318,6 +327,9 @@ public:
             //- Return the number of non-orthogonal correctors
             inline label nNonOrthCorr() const;
 
+            //- Return small delta
+            inline const dimensionedScalar& deltaSmall() const;
+
 
         // Thermo properties
 
@@ -345,6 +357,9 @@ public:
             //- Return the film wall velocity [m/s]
             virtual const volVectorField& Uw() const;
 
+            //- Return the film thickness*density (helper field) [kg/m3]
+            virtual const volScalarField& deltaRho() const;
+
             //- Return the film flux [kg.m/s]
             virtual const surfaceScalarField& phi() const;
 
@@ -445,6 +460,9 @@ public:
 
         // Sub-models
 
+            //- Film thermo
+            inline const filmThermoModel& filmThermo() const;
+
             //- Injection
             inline injectionModelList& injection();
 
@@ -460,6 +478,9 @@ public:
             //- Return the net film mass available over the next integration
             inline tmp<volScalarField> netMass() const;
 
+            //- Return the change in film mass due to sources/sinks
+            inline tmp<volScalarField> deltaMass() const;
+
             //- Return the gravity normal-to-patch component contribution
             inline tmp<volScalarField> gNorm() const;
 

src/regionModels/surfaceFilmModels/kinematicSingleLayer/kinematicSingleLayerI.H

@@ -23,6 +23,7 @@ License
 
 \*---------------------------------------------------------------------------*/
 
+#include "filmThermoModel.H"
 #include "surfaceInterpolate.H"
 #include "fvcSurfaceIntegrate.H"
 
@@ -61,6 +62,12 @@ inline label kinematicSingleLayer::nNonOrthCorr() const
 }
 
 
+inline const dimensionedScalar& kinematicSingleLayer::deltaSmall() const
+{
+    return deltaSmall_;
+}
+
+
 inline const volScalarField& kinematicSingleLayer::mu() const
 {
     return mu_;
@@ -163,6 +170,12 @@ inline const volScalarField& kinematicSingleLayer::muPrimary() const
 }
 
 
+inline const filmThermoModel& kinematicSingleLayer::filmThermo() const
+{
+    return filmThermo_();
+}
+
+
 inline injectionModelList& kinematicSingleLayer::injection()
 {
     return injection_;
@@ -183,14 +196,19 @@ inline tmp<volScalarField> kinematicSingleLayer::mass() const
 
 inline tmp<volScalarField> kinematicSingleLayer::netMass() const
 {
-    dimensionedScalar d0("SMALL", dimLength, ROOTVSMALL);
     return
-        fvc::surfaceSum(pos(phi_)*phi_/(fvc::interpolate(delta_) + d0))
+        fvc::surfaceSum(pos(phi_)*phi_/(fvc::interpolate(delta_) + deltaSmall_))
        *time().deltaT()
       + rho_*delta_*magSf();
 }
 
 
+inline tmp<volScalarField> kinematicSingleLayer::deltaMass() const
+{
+    return rhoSp_*magSf()*time().deltaT();
+}
+
+
 inline tmp<volScalarField> kinematicSingleLayer::gNorm() const
 {
     tmp<volScalarField> tgNorm

src/regionModels/surfaceFilmModels/submodels/kinematic/filmThermoModel/constantFilmThermo/constantFilmThermo.C

@@ -0,0 +1,398 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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 "constantFilmThermo.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace regionModels
+{
+namespace surfaceFilmModels
+{
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+defineTypeNameAndDebug(constantFilmThermo, 0);
+
+addToRunTimeSelectionTable
+(
+    filmThermoModel,
+    constantFilmThermo,
+    dictionary
+);
+
+
+void constantFilmThermo::init(thermoData& td)
+{
+    if (coeffs_.readIfPresent(td.name_, td.value_))
+    {
+        td.set_ = true;
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+constantFilmThermo::constantFilmThermo
+(
+    const surfaceFilmModel& owner,
+    const dictionary& dict
+)
+:
+    filmThermoModel(typeName, owner, dict),
+    name_(coeffs_.lookup("specieName")),
+    rho0_("rho0"),
+    mu0_("mu0"),
+    sigma0_("sigma0"),
+    Cp0_("Cp0"),
+    kappa0_("kappa0"),
+    hl0_("hl0"),
+    pv0_("pv0"),
+    W0_("W0"),
+    Tb0_("Tb0")
+{
+    init(rho0_);
+    init(mu0_);
+    init(sigma0_);
+    init(Cp0_);
+    init(kappa0_);
+    init(hl0_);
+    init(pv0_);
+    init(W0_);
+    init(Tb0_);
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+constantFilmThermo::~constantFilmThermo()
+{}
+
+
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+
+const word& constantFilmThermo::name() const
+{
+    return name_;
+}
+
+
+scalar constantFilmThermo::rho
+(
+    const scalar p,
+    const scalar T
+) const
+{
+    if (!rho0_.set_)
+    {
+        coeffs_.lookup(rho0_.name_) >> rho0_.value_;
+        rho0_.set_ = true;
+    }
+
+    return rho0_.value_;
+}
+
+
+scalar constantFilmThermo::mu
+(
+    const scalar p,
+    const scalar T
+) const
+{
+    if (!mu0_.set_)
+    {
+        coeffs_.lookup(mu0_.name_) >> mu0_.value_;
+        mu0_.set_ = true;
+    }
+
+    return mu0_.value_;
+}
+
+
+scalar constantFilmThermo::sigma
+(
+    const scalar p,
+    const scalar T
+) const
+{
+    if (!sigma0_.set_)
+    {
+        coeffs_.lookup(sigma0_.name_) >> sigma0_.value_;
+        sigma0_.set_ = true;
+    }
+
+    return sigma0_.value_;
+}
+
+
+scalar constantFilmThermo::Cp
+(
+    const scalar p,
+    const scalar T
+) const
+{
+    if (!Cp0_.set_)
+    {
+        coeffs_.lookup(Cp0_.name_) >> Cp0_.value_;
+        Cp0_.set_ = true;
+    }
+
+    return Cp0_.value_;
+}
+
+
+scalar constantFilmThermo::kappa
+(
+    const scalar p,
+    const scalar T
+) const
+{
+    if (!kappa0_.set_)
+    {
+        coeffs_.lookup(kappa0_.name_) >> kappa0_.value_;
+        kappa0_.set_ = true;
+    }
+
+    return kappa0_.value_;
+}
+
+
+scalar constantFilmThermo::D
+(
+    const scalar p,
+    const scalar T
+) const
+{
+    if (!D0_.set_)
+    {
+        coeffs_.lookup(D0_.name_) >> D0_.value_;
+        D0_.set_ = true;
+    }
+
+    return D0_.value_;
+}
+
+
+scalar constantFilmThermo::hl
+(
+    const scalar p,
+    const scalar T
+) const
+{
+    if (!hl0_.set_)
+    {
+        coeffs_.lookup(hl0_.name_) >> hl0_.value_;
+        hl0_.set_ = true;
+    }
+
+    return hl0_.value_;
+}
+
+
+scalar constantFilmThermo::pv
+(
+    const scalar p,
+    const scalar T
+) const
+{
+    if (!pv0_.set_)
+    {
+        coeffs_.lookup(pv0_.name_) >> pv0_.value_;
+        pv0_.set_ = true;
+    }
+
+    return pv0_.value_;
+}
+
+
+scalar constantFilmThermo::W() const
+{
+    if (!W0_.set_)
+    {
+        coeffs_.lookup(W0_.name_) >> W0_.value_;
+        W0_.set_ = true;
+    }
+
+    return W0_.value_;
+}
+
+
+scalar constantFilmThermo::Tb(const scalar p) const
+{
+    if (!Tb0_.set_)
+    {
+        coeffs_.lookup(Tb0_.name_) >> Tb0_.value_;
+        Tb0_.set_ = true;
+    }
+
+    return Tb0_.value_;
+}
+
+
+tmp<volScalarField> constantFilmThermo::rho() const
+{
+    tmp<volScalarField> trho
+    (
+        new volScalarField
+        (
+            IOobject
+            (
+                type() + ':' + rho0_.name_,
+                owner().time().timeName(),
+                owner().regionMesh(),
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            owner().regionMesh(),
+            dimensionedScalar("0", dimDensity, 0.0),
+            zeroGradientFvPatchScalarField::typeName
+        )
+    );
+
+    trho().internalField() = this->rho(0, 0);
+    trho().correctBoundaryConditions();
+
+    return trho;
+}
+
+
+tmp<volScalarField> constantFilmThermo::mu() const
+{
+    tmp<volScalarField> tmu
+    (
+        new volScalarField
+        (
+            IOobject
+            (
+                type() + ':' + mu0_.name_,
+                owner().time().timeName(),
+                owner().regionMesh(),
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            owner().regionMesh(),
+            dimensionedScalar("0", dimPressure*dimTime, 0.0),
+            zeroGradientFvPatchScalarField::typeName
+        )
+    );
+
+    tmu().internalField() = this->mu(0, 0);
+    tmu().correctBoundaryConditions();
+
+    return tmu;
+}
+
+
+tmp<volScalarField> constantFilmThermo::sigma() const
+{
+    tmp<volScalarField> tsigma
+    (
+        new volScalarField
+        (
+            IOobject
+            (
+                type() + ':' + sigma0_.name_,
+                owner().time().timeName(),
+                owner().regionMesh(),
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            owner().regionMesh(),
+            dimensionedScalar("0", dimMass/sqr(dimTime), 0.0),
+            zeroGradientFvPatchScalarField::typeName
+        )
+    );
+
+    tsigma().internalField() = this->sigma(0, 0);
+    tsigma().correctBoundaryConditions();
+
+    return tsigma;
+}
+
+
+tmp<volScalarField> constantFilmThermo::Cp() const
+{
+    tmp<volScalarField> tCp
+    (
+        new volScalarField
+        (
+            IOobject
+            (
+                type() + ':' + Cp0_.name_,
+                owner().time().timeName(),
+                owner().regionMesh(),
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            owner().regionMesh(),
+            dimensionedScalar("0", dimEnergy/dimMass/dimTemperature, 0.0),
+            zeroGradientFvPatchScalarField::typeName
+        )
+    );
+
+    tCp().internalField() = this->Cp(0, 0);
+    tCp().correctBoundaryConditions();
+
+    return tCp;
+}
+
+
+tmp<volScalarField> constantFilmThermo::kappa() const
+{
+    tmp<volScalarField> tkappa
+    (
+        new volScalarField
+        (
+            IOobject
+            (
+                type() + ':' + kappa0_.name_,
+                owner().time().timeName(),
+                owner().regionMesh(),
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            owner().regionMesh(),
+            dimensionedScalar("0", dimPower/dimLength/dimTemperature, 0.0),
+            zeroGradientFvPatchScalarField::typeName
+        )
+    );
+
+    tkappa().internalField() = this->kappa(0, 0);
+    tkappa().correctBoundaryConditions();
+
+    return tkappa;
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace surfaceFilmModels
+} // End namespace regionModels
+} // End namespace Foam
+
+// ************************************************************************* //

src/regionModels/surfaceFilmModels/submodels/kinematic/filmThermoModel/constantFilmThermo/constantFilmThermo.H

@@ -0,0 +1,224 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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::constantFilmThermo
+
+Description
+    Constant thermo model
+
+SourceFiles
+    constantFilmThermo.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef constantFilmThermo_H
+#define constantFilmThermo_H
+
+#include "filmThermoModel.H"
+#include "dimensionSet.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace regionModels
+{
+namespace surfaceFilmModels
+{
+
+/*---------------------------------------------------------------------------*\
+                     Class constantFilmThermo Declaration
+\*---------------------------------------------------------------------------*/
+
+class constantFilmThermo
+:
+    public filmThermoModel
+{
+public:
+
+    struct thermoData
+    {
+//    private:
+        word name_;
+        scalar value_;
+        bool set_;
+
+//    public:
+        thermoData()
+        :
+            name_("unknown"),
+            value_(0.0),
+            set_(false)
+        {}
+        thermoData(const word& n)
+        :
+            name_(n),
+            value_(0.0),
+            set_(false)
+        {}
+
+//        virtual ~thermoData()
+//        {}
+    };
+
+
+private:
+
+    // Private data
+
+        //- Specie name
+        word name_;
+
+        //- Density [kg/m3]
+        mutable thermoData rho0_;
+
+        //- Dynamic viscosity [Pa.s]
+        mutable thermoData mu0_;
+
+        //- Surface tension [kg/s2]
+        mutable thermoData sigma0_;
+
+        //- Specific heat capacity [J/kg/K]
+        mutable thermoData Cp0_;
+
+        //- Thermal conductivity [W/m/K]
+        mutable thermoData kappa0_;
+
+        //- Diffusivity [m2/s]
+        mutable thermoData D0_;
+
+        //- Latent heat [J/kg]
+        mutable thermoData hl0_;
+
+        //- Vapour pressure [Pa]
+        mutable thermoData pv0_;
+
+        //- Molecular weight [kg/kmol]
+        mutable thermoData W0_;
+
+        //- Boiling temperature [K]
+        mutable thermoData Tb0_;
+
+
+    // Private member functions
+
+        //- Initialise thermoData object
+        void init(thermoData& td);
+
+        //- Disallow default bitwise copy construct
+        constantFilmThermo(const constantFilmThermo&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const constantFilmThermo&);
+
+
+public:
+
+    //- Runtime type information
+    TypeName("constant");
+
+
+    // Constructors
+
+        //- Construct from surface film model and dictionary
+        constantFilmThermo
+        (
+            const surfaceFilmModel& owner,
+            const dictionary& dict
+        );
+
+
+    //- Destructor
+    virtual ~constantFilmThermo();
+
+
+    // Member Functions
+
+        //- Return the specie name
+        virtual const word& name() const;
+
+
+        // Elemental access
+
+            //- Return density [kg/m3]
+            virtual scalar rho(const scalar p, const scalar T) const;
+
+            //- Return dynamic viscosity [Pa.s]
+            virtual scalar mu(const scalar p, const scalar T) const;
+
+            //- Return surface tension [kg/s2]
+            virtual scalar sigma(const scalar p, const scalar T) const;
+
+            //- Return specific heat capacity [J/kg/K]
+            virtual scalar Cp(const scalar p, const scalar T) const;
+
+            //- Return thermal conductivity [W/m/K]
+            virtual scalar kappa(const scalar p, const scalar T) const;
+
+            //- Return diffusivity [m2/s]
+            virtual scalar D(const scalar p, const scalar T) const;
+
+            //- Return latent heat [J/kg]
+            virtual scalar hl(const scalar p, const scalar T) const;
+
+            //- Return vapour pressure [Pa]
+            virtual scalar pv(const scalar p, const scalar T) const;
+
+            //- Return molecular weight [kg/kmol]
+            virtual scalar W() const;
+
+            //- Return boiling temperature [K]
+            virtual scalar Tb(const scalar p) const;
+
+
+        // Field access
+
+            //- Return density [kg/m3]
+            virtual tmp<volScalarField> rho() const;
+
+            //- Return dynamic viscosity [Pa.s]
+            virtual tmp<volScalarField> mu() const;
+
+            //- Return surface tension [kg/s2]
+            virtual tmp<volScalarField> sigma() const;
+
+            //- Return specific heat capacity [J/kg/K]
+            virtual tmp<volScalarField> Cp() const;
+
+            //- Return thermal conductivity [W/m/K]
+            virtual tmp<volScalarField> kappa() const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace surfaceFilmModels
+} // End namespace regionModels
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/regionModels/surfaceFilmModels/submodels/kinematic/filmThermoModel/filmThermoModel/filmThermoModel.C

@@ -0,0 +1,76 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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 "filmThermoModel.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace regionModels
+{
+namespace surfaceFilmModels
+{
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+defineTypeNameAndDebug(filmThermoModel, 0);
+defineRunTimeSelectionTable(filmThermoModel, dictionary);
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+filmThermoModel::filmThermoModel
+(
+    const surfaceFilmModel& owner
+)
+:
+    subModelBase(owner)
+{}
+
+
+filmThermoModel::filmThermoModel
+(
+    const word& type,
+    const surfaceFilmModel& owner,
+    const dictionary& dict
+)
+:
+    subModelBase(type, owner, dict)
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+filmThermoModel::~filmThermoModel()
+{}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace surfaceFilmModels
+} // End namespace regionModels
+} // End namespace Foam
+
+// ************************************************************************* //

src/regionModels/surfaceFilmModels/submodels/kinematic/filmThermoModel/filmThermoModel/filmThermoModel.H

@@ -0,0 +1,186 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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::filmThermoModel
+
+Description
+    Base class for film thermo models
+
+SourceFiles
+    filmThermoModel.C
+    filmThermoModelNew.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef filmThermoModel_H
+#define filmThermoModel_H
+
+#include "subModelBase.H"
+#include "runTimeSelectionTables.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace regionModels
+{
+namespace surfaceFilmModels
+{
+
+/*---------------------------------------------------------------------------*\
+                       Class filmThermoModel Declaration
+\*---------------------------------------------------------------------------*/
+
+class filmThermoModel
+:
+    public subModelBase
+{
+private:
+
+    // Private Member Functions
+
+        //- Disallow default bitwise copy construct
+        filmThermoModel(const filmThermoModel&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const filmThermoModel&);
+
+
+public:
+
+    //- Runtime type information
+    TypeName("filmThermoModel");
+
+
+    // Declare runtime constructor selection table
+
+         declareRunTimeSelectionTable
+         (
+             autoPtr,
+             filmThermoModel,
+             dictionary,
+             (
+                 const surfaceFilmModel& owner,
+                 const dictionary& dict
+             ),
+             (owner, dict)
+         );
+
+    // Constructors
+
+        //- Construct null
+        filmThermoModel(const surfaceFilmModel& owner);
+
+        //- Construct from type name, dictionary and surface film model
+        filmThermoModel
+        (
+            const word& type,
+            const surfaceFilmModel& owner,
+            const dictionary& dict
+        );
+
+
+    // Selectors
+
+        //- Return a reference to the selected phase change model
+        static autoPtr<filmThermoModel> New
+        (
+            const surfaceFilmModel& owner,
+            const dictionary& dict
+        );
+
+
+    //- Destructor
+    virtual ~filmThermoModel();
+
+
+    // Member Functions
+
+        //- Return the specie name
+        virtual const word& name() const = 0;
+
+
+        // Elemental access
+
+            //- Return density [kg/m3]
+            virtual scalar rho(const scalar p, const scalar T) const = 0;
+
+            //- Return dynamic viscosity [Pa.s]
+            virtual scalar mu(const scalar p, const scalar T) const = 0;
+
+            //- Return surface tension [kg/s2]
+            virtual scalar sigma(const scalar p, const scalar T) const = 0;
+
+            //- Return specific heat capacity [J/kg/K]
+            virtual scalar Cp(const scalar p, const scalar T) const = 0;
+
+            //- Return thermal conductivity [W/m/K]
+            virtual scalar kappa(const scalar p, const scalar T) const = 0;
+
+            //- Return diffusivity [m2/s]
+            virtual scalar D(const scalar p, const scalar T) const = 0;
+
+            //- Return latent heat [J/kg]
+            virtual scalar hl(const scalar p, const scalar T) const = 0;
+
+            //- Return vapour pressure [Pa]
+            virtual scalar pv(const scalar p, const scalar T) const = 0;
+
+            //- Return molecular weight [kg/kmol]
+            virtual scalar W() const = 0;
+
+            //- Return boiling temperature [K]
+            virtual scalar Tb(const scalar p) const = 0;
+
+
+        // Field access
+
+            //- Return density [kg/m3]
+            virtual tmp<volScalarField> rho() const = 0;
+
+            //- Return dynamic viscosity [Pa.s]
+            virtual tmp<volScalarField> mu() const = 0;
+
+            //- Return surface tension [kg/s2]
+            virtual tmp<volScalarField> sigma() const = 0;
+
+            //- Return specific heat capacity [J/kg/K]
+            virtual tmp<volScalarField> Cp() const = 0;
+
+            //- Return thermal conductivity [W/m/K]
+            virtual tmp<volScalarField> kappa() const = 0;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace surfaceFilmModels
+} // End namespace regionModels
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/regionModels/surfaceFilmModels/submodels/kinematic/filmThermoModel/filmThermoModel/filmThermoModelNew.C

@@ -0,0 +1,73 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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 "filmThermoModel.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace regionModels
+{
+namespace surfaceFilmModels
+{
+
+// * * * * * * * * * * * * * * * * Selectors * * * * * * * * * * * * * * * * //
+
+autoPtr<filmThermoModel> filmThermoModel::New
+(
+    const surfaceFilmModel& model,
+    const dictionary& dict
+)
+{
+    word modelType(dict.lookup("filmThermoModel"));
+
+    Info<< "    Selecting filmThermoModel " << modelType << endl;
+
+    dictionaryConstructorTable::iterator cstrIter =
+        dictionaryConstructorTablePtr_->find(modelType);
+
+    if (cstrIter == dictionaryConstructorTablePtr_->end())
+    {
+        FatalErrorIn
+        (
+            "filmThermoModel::New(const surfaceFilmModel&, const dictionary&)"
+        )   << "Unknown filmThermoModel type " << modelType << nl << nl
+            << "Valid filmThermoModel types are:" << nl
+            << dictionaryConstructorTablePtr_->toc()
+            << exit(FatalError);
+    }
+
+    return autoPtr<filmThermoModel>(cstrIter()(model, dict));
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace surfaceFilmModels
+} // End namespace regionModels
+} // End namespace Foam
+
+// ************************************************************************* //

src/regionModels/surfaceFilmModels/submodels/kinematic/filmThermoModel/liquidFilmThermo/liquidFilmThermo.C

@@ -0,0 +1,486 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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 "liquidFilmThermo.H"
+#include "demandDrivenData.H"
+#include "thermoSingleLayer.H"
+#include "SLGThermo.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace regionModels
+{
+namespace surfaceFilmModels
+{
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+defineTypeNameAndDebug(liquidFilmThermo, 0);
+
+addToRunTimeSelectionTable
+(
+    filmThermoModel,
+    liquidFilmThermo,
+    dictionary
+);
+
+
+// * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
+
+const thermoSingleLayer& liquidFilmThermo::thermoFilm() const
+{
+    if (!isA<thermoSingleLayer>(owner_))
+    {
+        FatalErrorIn
+        (
+            "const thermoSingleLayer& liquidFilmThermo::thermoFilm() const"
+        )
+            << "Thermo model requires a " << thermoSingleLayer::typeName
+            << " film to supply the pressure and temperature, but "
+            << owner_.type() << " film model selected.  "
+            << "Use the 'useReferenceValues' flag to employ reference "
+            << "pressure and temperature" << exit(FatalError);
+    }
+
+    return refCast<const thermoSingleLayer>(owner_);
+}
+
+
+void liquidFilmThermo::initLiquid(const dictionary& dict)
+{
+    if (liquidPtr_ != NULL)
+    {
+        return;
+    }
+
+    dict.lookup("liquid") >> name_;
+
+    if (owner_.primaryMesh().foundObject<SLGThermo>("SLGThermo"))
+    {
+        // retrieve from film thermo
+        ownLiquid_ = false;
+
+        const SLGThermo& thermo =
+            owner_.primaryMesh().lookupObject<SLGThermo>("SLGThermo");
+        label id = thermo.liquidId(name_);
+        liquidPtr_ = &thermo.liquids().properties()[id];
+    }
+    else
+    {
+        // new liquid create
+        ownLiquid_ = true;
+
+        liquidPtr_ = new liquidProperties(dict.subDict(name_ + "Coeffs"));
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+liquidFilmThermo::liquidFilmThermo
+(
+    const surfaceFilmModel& owner,
+    const dictionary& dict
+)
+:
+    filmThermoModel(typeName, owner, dict),
+    name_("unknown_liquid"),
+    liquidPtr_(NULL),
+    ownLiquid_(false),
+    useReferenceValues_(readBool(coeffs_.lookup("useReferenceValues"))),
+    pRef_(0.0),
+    TRef_(0.0)
+{
+    initLiquid(coeffs_);
+
+    if (useReferenceValues_)
+    {
+        coeffs_.lookup("pRef") >> pRef_;
+        coeffs_.lookup("TRef") >> TRef_;
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+liquidFilmThermo::~liquidFilmThermo()
+{
+    if (ownLiquid_)
+    {
+        deleteDemandDrivenData(liquidPtr_);
+    }
+}
+
+
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+
+const word& liquidFilmThermo::name() const
+{
+    return name_;
+}
+
+
+scalar liquidFilmThermo::rho
+(
+    const scalar p,
+    const scalar T
+) const
+{
+    return liquidPtr_->rho(p, T);
+}
+
+
+scalar liquidFilmThermo::mu
+(
+    const scalar p,
+    const scalar T
+) const
+{
+    return liquidPtr_->mu(p, T);
+}
+
+
+scalar liquidFilmThermo::sigma
+(
+    const scalar p,
+    const scalar T
+) const
+{
+    return liquidPtr_->sigma(p, T);
+}
+
+
+scalar liquidFilmThermo::Cp
+(
+    const scalar p,
+    const scalar T
+) const
+{
+    return liquidPtr_->Cp(p, T);
+}
+
+
+scalar liquidFilmThermo::kappa
+(
+    const scalar p,
+    const scalar T
+) const
+{
+    return liquidPtr_->K(p, T);
+}
+
+
+scalar liquidFilmThermo::D
+(
+    const scalar p,
+    const scalar T
+) const
+{
+    return liquidPtr_->D(p, T);
+}
+
+
+scalar liquidFilmThermo::hl
+(
+    const scalar p,
+    const scalar T
+) const
+{
+    return liquidPtr_->hl(p, T);
+}
+
+
+scalar liquidFilmThermo::pv
+(
+    const scalar p,
+    const scalar T
+) const
+{
+    return liquidPtr_->pv(p, T);
+}
+
+
+scalar liquidFilmThermo::W() const
+{
+    return liquidPtr_->W();
+}
+
+
+scalar liquidFilmThermo::Tb(const scalar p) const
+{
+    return liquidPtr_->pvInvert(p);
+}
+
+
+tmp<volScalarField> liquidFilmThermo::rho() const
+{
+    tmp<volScalarField> trho
+    (
+        new volScalarField
+        (
+            IOobject
+            (
+                type() + ":rho",
+                owner().time().timeName(),
+                owner().regionMesh(),
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            owner().regionMesh(),
+            dimensionedScalar("0", dimDensity, 0.0),
+            zeroGradientFvPatchScalarField::typeName
+        )
+    );
+
+    scalarField& rho = trho().internalField();
+
+    if (useReferenceValues_)
+    {
+        forAll(rho, cellI)
+        {
+            rho[cellI] = this->rho(pRef_, TRef_);
+        }
+    }
+    else
+    {
+        const thermoSingleLayer& film = thermoFilm();
+
+        const volScalarField& T = film.T();
+        const volScalarField& p = film.pPrimary();
+
+        forAll(rho, cellI)
+        {
+            rho[cellI] = this->rho(p[cellI], T[cellI]);
+        }
+    }
+
+    trho().correctBoundaryConditions();
+
+    return trho;
+}
+
+
+tmp<volScalarField> liquidFilmThermo::mu() const
+{
+    tmp<volScalarField> tmu
+    (
+        new volScalarField
+        (
+            IOobject
+            (
+                type() + ":mu",
+                owner().time().timeName(),
+                owner().regionMesh(),
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            owner().regionMesh(),
+            dimensionedScalar("0", dimPressure*dimTime, 0.0),
+            zeroGradientFvPatchScalarField::typeName
+        )
+    );
+
+    scalarField& mu = tmu().internalField();
+
+    if (useReferenceValues_)
+    {
+        forAll(mu, cellI)
+        {
+            mu[cellI] = this->mu(pRef_, TRef_);
+        }
+    }
+    else
+    {
+        const thermoSingleLayer& film = thermoFilm();
+
+        const volScalarField& T = film.T();
+        const volScalarField& p = film.pPrimary();
+
+        forAll(mu, cellI)
+        {
+            mu[cellI] = this->mu(p[cellI], T[cellI]);
+        }
+    }
+
+    tmu().correctBoundaryConditions();
+
+    return tmu;
+}
+
+
+tmp<volScalarField> liquidFilmThermo::sigma() const
+{
+    tmp<volScalarField> tsigma
+    (
+        new volScalarField
+        (
+            IOobject
+            (
+                type() + ":sigma",
+                owner().time().timeName(),
+                owner().regionMesh(),
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            owner().regionMesh(),
+            dimensionedScalar("0", dimMass/sqr(dimTime), 0.0),
+            zeroGradientFvPatchScalarField::typeName
+        )
+    );
+
+    scalarField& sigma = tsigma().internalField();
+
+    if (useReferenceValues_)
+    {
+        forAll(sigma, cellI)
+        {
+            sigma[cellI] = this->sigma(pRef_, TRef_);
+        }
+    }
+    else
+    {
+        const thermoSingleLayer& film = thermoFilm();
+
+        const volScalarField& T = film.T();
+        const volScalarField& p = film.pPrimary();
+
+        forAll(sigma, cellI)
+        {
+            sigma[cellI] = this->sigma(p[cellI], T[cellI]);
+        }
+    }
+
+    tsigma().correctBoundaryConditions();
+
+    return tsigma;
+}
+
+
+tmp<volScalarField> liquidFilmThermo::Cp() const
+{
+    tmp<volScalarField> tCp
+    (
+        new volScalarField
+        (
+            IOobject
+            (
+                type() + ":Cp",
+                owner().time().timeName(),
+                owner().regionMesh(),
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            owner().regionMesh(),
+            dimensionedScalar("0", dimEnergy/dimMass/dimTemperature, 0.0),
+            zeroGradientFvPatchScalarField::typeName
+        )
+    );
+
+    scalarField& Cp = tCp().internalField();
+
+    if (useReferenceValues_)
+    {
+        forAll(Cp, cellI)
+        {
+            Cp[cellI] = this->Cp(pRef_, TRef_);
+        }
+    }
+    else
+    {
+        const thermoSingleLayer& film = thermoFilm();
+
+        const volScalarField& T = film.T();
+        const volScalarField& p = film.pPrimary();
+
+        forAll(Cp, cellI)
+        {
+            Cp[cellI] = this->Cp(p[cellI], T[cellI]);
+        }
+    }
+
+    tCp().correctBoundaryConditions();
+
+    return tCp;
+}
+
+
+tmp<volScalarField> liquidFilmThermo::kappa() const
+{
+    tmp<volScalarField> tkappa
+    (
+        new volScalarField
+        (
+            IOobject
+            (
+                type() + ":kappa",
+                owner().time().timeName(),
+                owner().regionMesh(),
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            owner().regionMesh(),
+            dimensionedScalar("0", dimPower/dimLength/dimTemperature, 0.0),
+            zeroGradientFvPatchScalarField::typeName
+        )
+    );
+
+    scalarField& kappa = tkappa().internalField();
+
+    if (useReferenceValues_)
+    {
+        forAll(kappa, cellI)
+        {
+            kappa[cellI] = this->kappa(pRef_, TRef_);
+        }
+    }
+    else
+    {
+        const thermoSingleLayer& film = thermoFilm();
+
+        const volScalarField& T = film.T();
+        const volScalarField& p = film.pPrimary();
+
+        forAll(kappa, cellI)
+        {
+            kappa[cellI] = this->kappa(p[cellI], T[cellI]);
+        }
+    }
+
+    tkappa().correctBoundaryConditions();
+
+    return tkappa;
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace surfaceFilmModels
+} // End namespace regionModels
+} // End namespace Foam
+
+// ************************************************************************* //

src/regionModels/surfaceFilmModels/submodels/kinematic/filmThermoModel/liquidFilmThermo/liquidFilmThermo.H

@@ -0,0 +1,187 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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::liquidFilmThermo
+
+Description
+    Liquid thermo model
+
+SourceFiles
+    liquidFilmThermo.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef liquidFilmThermo_H
+#define liquidFilmThermo_H
+
+#include "filmThermoModel.H"
+#include "liquidProperties.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace regionModels
+{
+namespace surfaceFilmModels
+{
+
+// Forward class declarations
+class thermoSingleLayer;
+
+/*---------------------------------------------------------------------------*\
+                      Class liquidFilmThermo Declaration
+\*---------------------------------------------------------------------------*/
+
+class liquidFilmThermo
+:
+    public filmThermoModel
+{
+protected:
+
+    // Protected data
+
+        //- Liquid name
+        word name_;
+
+        //- Pointer to the liquid properties
+        const liquidProperties* liquidPtr_;
+
+        //- Flag to indicate that model owns the liquid object
+        bool ownLiquid_;
+
+        //- Flag to indicate that reference values of p and T should be used
+        bool useReferenceValues_;
+
+        //- Reference pressure [pa]
+        scalar pRef_;
+
+        //- Reference temperature [K]
+        scalar TRef_;
+
+
+    // Protected member functions
+
+        //- Return a reference to a thermo film
+        const thermoSingleLayer& thermoFilm() const;
+
+        //- Initialise the liquid pointer
+        void initLiquid(const dictionary& dict);
+
+        //- Disallow default bitwise copy construct
+        liquidFilmThermo(const liquidFilmThermo&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const liquidFilmThermo&);
+
+
+public:
+
+    //- Runtime type information
+    TypeName("liquid");
+
+
+    // Constructors
+
+        //- Construct from surface film model and dictionary
+        liquidFilmThermo
+        (
+            const surfaceFilmModel& owner,
+            const dictionary& dict
+        );
+
+
+    //- Destructor
+    virtual ~liquidFilmThermo();
+
+
+    // Member Functions
+
+        //- Return the specie name
+        virtual const word& name() const;
+
+
+        // Elemental access
+
+            //- Return density [kg/m3]
+            virtual scalar rho(const scalar p, const scalar T) const;
+
+            //- Return dynamic viscosity [Pa.s]
+            virtual scalar mu(const scalar p, const scalar T) const;
+
+            //- Return surface tension [kg/s2]
+            virtual scalar sigma(const scalar p, const scalar T) const;
+
+            //- Return specific heat capacity [J/kg/K]
+            virtual scalar Cp(const scalar p, const scalar T) const;
+
+            //- Return thermal conductivity [W/m/K]
+            virtual scalar kappa(const scalar p, const scalar T) const;
+
+            //- Return diffusivity [m2/s]
+            virtual scalar D(const scalar p, const scalar T) const;
+
+            //- Return latent heat [J/kg]
+            virtual scalar hl(const scalar p, const scalar T) const;
+
+            //- Return vapour pressure [Pa]
+            virtual scalar pv(const scalar p, const scalar T) const;
+
+            //- Return molecular weight [kg/kmol]
+            virtual scalar W() const;
+
+            //- Return boiling temperature [K]
+            virtual scalar Tb(const scalar p) const;
+
+
+        // Field access
+
+            //- Return density [kg/m3]
+            virtual tmp<volScalarField> rho() const;
+
+            //- Return dynamic viscosity [Pa.s]
+            virtual tmp<volScalarField> mu() const;
+
+            //- Return surface tension [kg/s2]
+            virtual tmp<volScalarField> sigma() const;
+
+            //- Return specific heat capacity [J/kg/K]
+            virtual tmp<volScalarField> Cp() const;
+
+            //- Return thermal conductivity [W/m/K]
+            virtual tmp<volScalarField> kappa() const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace surfaceFilmModels
+} // End namespace regionModels
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/regionModels/surfaceFilmModels/submodels/kinematic/filmTurbulenceModel/filmTurbulenceModel/filmTurbulenceModelNew.C

@@ -44,7 +44,7 @@ autoPtr<filmTurbulenceModel> filmTurbulenceModel::New
 {
     const word modelType(dict.lookup("turbulence"));
 
-    Info<< "        " << modelType << endl;
+    Info<< "    Selecting filmTurbulenceModel " << modelType << endl;
 
     dictionaryConstructorTable::iterator cstrIter =
         dictionaryConstructorTablePtr_->find(modelType);

src/regionModels/surfaceFilmModels/submodels/kinematic/filmTurbulenceModel/laminar/laminar.C

@@ -122,6 +122,7 @@ void laminar::correct()
     // do nothing
 }
 
+
 tmp<fvVectorMatrix> laminar::Su(volVectorField& U) const
 {
     // local reference to film model
@@ -138,9 +139,7 @@ tmp<fvVectorMatrix> laminar::Su(volVectorField& U) const
 
     // employ simple coeff-based model
     volScalarField Cs("Cs", Cf_*rhop*mag(Up - Us));
-
+    volScalarField Cw("Cw", mu/(0.3333*(delta + film.deltaSmall())));
-    dimensionedScalar d0("SMALL", delta.dimensions(), SMALL);
-    volScalarField Cw("Cw", mu/(0.3333*(delta + d0)));
     Cw.min(5000.0);
 
     return

src/regionModels/surfaceFilmModels/submodels/subModelBase.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
@@ -99,11 +99,16 @@ public:
 
         // Access
 
+            virtual const word& type() const = 0;
+
             //- Return the reference to the owner surface film model
             inline const surfaceFilmModel& owner() const;
 
             //- Return the model coefficients dictionary
             inline const dictionary& coeffs() const;
+
+            template<class FilmType>
+            inline const FilmType& filmType() const;
 };
 
 
@@ -119,6 +124,12 @@ public:
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
+#ifdef NoRepository
+    #include "subModelBaseTemplates.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
 #endif
 
 // ************************************************************************* //

src/regionModels/surfaceFilmModels/submodels/subModelBaseTemplates.C

@@ -0,0 +1,56 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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/>.
+
+\*---------------------------------------------------------------------------*/
+
+namespace Foam
+{
+namespace regionModels
+{
+namespace surfaceFilmModels
+{
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+template<class FilmType>
+const FilmType& subModelBase::filmType() const
+{
+    if (!isA<FilmType>(owner_))
+    {
+        FatalErrorIn("FilmType& subModelBase::film() const")
+            << "Model " << type() << " requested film type "
+            << FilmType::typeName << " but film is type " << owner_.type()
+            << abort(FatalError);
+    }
+
+    return refCast<const FilmType>(owner_);
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace surfaceFilmModels
+} // End namespace regionModels
+} // End namespace Foam
+
+// ************************************************************************* //

src/regionModels/surfaceFilmModels/submodels/thermo/filmViscosityModel/constantViscosity/constantViscosity.C

@@ -0,0 +1,90 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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 "constantViscosity.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace regionModels
+{
+namespace surfaceFilmModels
+{
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+defineTypeNameAndDebug(constantViscosity, 0);
+
+addToRunTimeSelectionTable
+(
+    filmViscosityModel,
+    constantViscosity,
+    dictionary
+);
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+constantViscosity::constantViscosity
+(
+    const surfaceFilmModel& owner,
+    const dictionary& dict,
+    volScalarField& mu
+)
+:
+    filmViscosityModel(typeName, owner, dict, mu),
+    mu0_(readScalar(coeffs().lookup("mu0")))
+{
+    mu_.internalField() = mu0_;
+    mu_.correctBoundaryConditions();
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+constantViscosity::~constantViscosity()
+{}
+
+
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+
+void constantViscosity::correct
+(
+    const volScalarField& p,
+    const volScalarField& T
+)
+{
+    // do nothing
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace surfaceFilmModels
+} // End namespace regionModels
+} // End namespace Foam
+
+// ************************************************************************* //

src/regionModels/surfaceFilmModels/submodels/thermo/filmViscosityModel/constantViscosity/constantViscosity.H

@@ -0,0 +1,121 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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::constantViscosity
+
+Description
+    Constant viscosity model
+
+SourceFiles
+    constantViscosity.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef constantViscosity_H
+#define constantViscosity_H
+
+#include "filmViscosityModel.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace regionModels
+{
+namespace surfaceFilmModels
+{
+
+/*---------------------------------------------------------------------------*\
+                      Class constantViscosity Declaration
+\*---------------------------------------------------------------------------*/
+
+class constantViscosity
+:
+    public filmViscosityModel
+{
+private:
+
+    // Private member functions
+
+        //- Disallow default bitwise copy construct
+        constantViscosity(const constantViscosity&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const constantViscosity&);
+
+
+protected:
+
+    // Protected data
+
+        //- Constant viscosity [Pa.s]
+        scalar mu0_;
+
+
+
+public:
+
+    //- Runtime type information
+    TypeName("constant");
+
+
+    // Constructors
+
+        //- Construct from surface film model
+        constantViscosity
+        (
+            const surfaceFilmModel& owner,
+            const dictionary& dict,
+            volScalarField& mu
+        );
+
+
+    //- Destructor
+    virtual ~constantViscosity();
+
+
+    // Member Functions
+
+        // Evolution
+
+            //- Correct
+            virtual void correct
+            (
+                const volScalarField& p,
+                const volScalarField& T
+            );
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace surfaceFilmModels
+} // End namespace regionModels
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/regionModels/surfaceFilmModels/submodels/thermo/filmViscosityModel/filmViscosityModel/filmViscosityModel.C

@@ -0,0 +1,75 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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 "filmViscosityModel.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace regionModels
+{
+namespace surfaceFilmModels
+{
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+defineTypeNameAndDebug(filmViscosityModel, 0);
+defineRunTimeSelectionTable(filmViscosityModel, dictionary);
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+filmViscosityModel::filmViscosityModel
+(
+    const word& type,
+    const surfaceFilmModel& owner,
+    const dictionary& dict,
+    volScalarField& mu
+)
+:
+    subModelBase(type, owner, dict),
+    mu_(mu)
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+filmViscosityModel::~filmViscosityModel()
+{}
+
+
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+
+void filmViscosityModel::info(Ostream& os) const
+{}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // end namespace surfaceFilmModels
+} // end namespace regionModels
+} // end namespace Foam
+
+// ************************************************************************* //

src/regionModels/surfaceFilmModels/submodels/thermo/filmViscosityModel/filmViscosityModel/filmViscosityModel.H

@@ -0,0 +1,156 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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::filmViscosityModel
+
+Description
+    Base class for surface film viscosity models
+
+SourceFiles
+    filmViscosityModel.C
+    filmViscosityModelNew.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef filmViscosityModel_H
+#define filmViscosityModel_H
+
+#include "subModelBase.H"
+#include "runTimeSelectionTables.H"
+#include "scalarField.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace regionModels
+{
+namespace surfaceFilmModels
+{
+
+/*---------------------------------------------------------------------------*\
+                     Class filmViscosityModel Declaration
+\*---------------------------------------------------------------------------*/
+
+class filmViscosityModel
+:
+    public subModelBase
+{
+private:
+
+    // Private Member Functions
+
+        //- Disallow default bitwise copy construct
+        filmViscosityModel(const filmViscosityModel&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const filmViscosityModel&);
+
+
+protected:
+
+    // Protected Member Data
+
+        //- Reference to the viscosity field
+        volScalarField& mu_;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("filmViscosityModel");
+
+
+    // Declare runtime constructor selection table
+
+         declareRunTimeSelectionTable
+         (
+             autoPtr,
+             filmViscosityModel,
+             dictionary,
+             (
+                const surfaceFilmModel& owner,
+                const dictionary& dict,
+                volScalarField& mu
+             ),
+             (owner, dict, mu)
+         );
+
+    // Constructors
+
+        //- Construct from type name, dictionary and surface film model
+        filmViscosityModel
+        (
+            const word& type,
+            const surfaceFilmModel& owner,
+            const dictionary& dict,
+            volScalarField& mu
+        );
+
+
+    // Selectors
+
+        //- Return a reference to the selected phase change model
+        static autoPtr<filmViscosityModel> New
+        (
+            const surfaceFilmModel& owner,
+            const dictionary& dict,
+            volScalarField& mu
+        );
+
+
+    //- Destructor
+    virtual ~filmViscosityModel();
+
+
+    // Member Functions
+
+        // Evolution
+
+            //- Correct
+            virtual void correct
+            (
+                const volScalarField& p,
+                const volScalarField& T
+            ) = 0;
+
+
+        // I-O
+
+            //- Provide some feedback
+            virtual void info(Ostream& os) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace surfaceFilmModels
+} // End namespace regionModels
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/regionModels/surfaceFilmModels/submodels/thermo/filmViscosityModel/filmViscosityModel/filmViscosityModelNew.C

@@ -0,0 +1,79 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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 "filmViscosityModel.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace regionModels
+{
+namespace surfaceFilmModels
+{
+
+// * * * * * * * * * * * * * * * * Selectors * * * * * * * * * * * * * * * * //
+
+autoPtr<filmViscosityModel> filmViscosityModel::New
+(
+    const surfaceFilmModel& model,
+    const dictionary& dict,
+    volScalarField& mu
+)
+{
+    word modelType(dict.lookup("filmViscosityModel"));
+
+    Info<< "    Selecting filmViscosityModel " << modelType << endl;
+
+    dictionaryConstructorTable::iterator cstrIter =
+        dictionaryConstructorTablePtr_->find(modelType);
+
+    if (cstrIter == dictionaryConstructorTablePtr_->end())
+    {
+        FatalErrorIn
+        (
+            "filmViscosityModel::New"
+            "("
+                "const surfaceFilmModel&, "
+                "const dictionary&, "
+                "volScalarField&"
+            ")"
+        )   << "Unknown filmViscosityModel type " << modelType
+            << nl << nl << "Valid filmViscosityModel types are:" << nl
+            << dictionaryConstructorTablePtr_->toc()
+            << exit(FatalError);
+    }
+
+    return autoPtr<filmViscosityModel>(cstrIter()(model, dict, mu));
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // end namespace surfaceFilmModels
+} // end namespace regionModels
+} // end namespace Foam
+
+// ************************************************************************* //

src/regionModels/surfaceFilmModels/submodels/thermo/filmViscosityModel/liquidViscosity/liquidViscosity.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -23,11 +23,9 @@ License
 
 \*---------------------------------------------------------------------------*/
 
-#include "surfaceShearForce.H"
+#include "liquidViscosity.H"
+#include "thermoSingleLayer.H"
 #include "addToRunTimeSelectionTable.H"
-#include "fvmSup.H"
-#include "kinematicSingleLayer.H"
-#include "turbulenceModel.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -40,55 +38,45 @@ namespace surfaceFilmModels
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
-defineTypeNameAndDebug(surfaceShearForce, 0);
+defineTypeNameAndDebug(liquidViscosity, 0);
-addToRunTimeSelectionTable(force, surfaceShearForce, dictionary);
+
+addToRunTimeSelectionTable
+(
+    filmViscosityModel,
+    liquidViscosity,
+    dictionary
+);
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
-surfaceShearForce::surfaceShearForce
+liquidViscosity::liquidViscosity
 (
     const surfaceFilmModel& owner,
-    const dictionary& dict
+    const dictionary& dict,
+    volScalarField& mu
 )
 :
-    force(typeName, owner, dict),
+    filmViscosityModel(typeName, owner, dict, mu)
-    Cf_(readScalar(coeffs_.lookup("Cf")))
 {}
 
 
 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 
-surfaceShearForce::~surfaceShearForce()
+liquidViscosity::~liquidViscosity()
 {}
 
 
 // * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
 
-tmp<fvVectorMatrix> surfaceShearForce::correct(volVectorField& U)
+void liquidViscosity::correct
+(
+    const volScalarField& p,
+    const volScalarField& T
+)
 {
-    // local reference to film model
+    const thermoSingleLayer& film = filmType<thermoSingleLayer>();
-    const kinematicSingleLayer& film =
-        static_cast<const kinematicSingleLayer&>(owner_);
 
-    // local references to film fields
+    mu_ = film.filmThermo().mu();
-    const volScalarField& mu = film.mu();
-    const volVectorField& Uw = film.Uw();
-    const volScalarField& delta = film.delta();
-    const volVectorField& Up = film.UPrimary();
-
-    // laminar case - employ simple coeff-based model
-    const volScalarField& rhop = film.rhoPrimary();
-    volScalarField Cs("Cs", Cf_*rhop*mag(Up - U));
-
-    dimensionedScalar d0("SMALL", delta.dimensions(), SMALL);
-    volScalarField Cw("Cw", mu/(0.3333*(delta + d0)));
-    Cw.min(5000.0);
-
-    return
-    (
-       - fvm::Sp(Cs, U) + Cs*Up // surface contribution
-       - fvm::Sp(Cw, U) + Cw*Uw // wall contribution
-    );
 }
 
 

src/regionModels/surfaceFilmModels/submodels/thermo/filmViscosityModel/liquidViscosity/liquidViscosity.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) 2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -22,20 +22,20 @@ License
     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
 
 Class
-    Foam::surfaceShearForce
+    Foam::liquidViscosity
 
 Description
-    Film surface shear force
+    liquidViscosity viscosity model
 
 SourceFiles
-    surfaceShearForce.C
+    liquidViscosity.C
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef surfaceShearForce_H
+#ifndef liquidViscosity_H
-#define surfaceShearForce_H
+#define liquidViscosity_H
 
-#include "force.H"
+#include "filmViscosityModel.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -47,49 +47,49 @@ namespace surfaceFilmModels
 {
 
 /*---------------------------------------------------------------------------*\
-                      Class surfaceShearForce Declaration
+                       Class liquidViscosity Declaration
 \*---------------------------------------------------------------------------*/
 
-class surfaceShearForce
+class liquidViscosity
 :
-    public force
+    public filmViscosityModel
 {
 private:
 
-    // Private Data
-
-        //- Surface roughness coefficient
-        scalar Cf_;
-
-
-
     // Private member functions
 
         //- Disallow default bitwise copy construct
-        surfaceShearForce(const surfaceShearForce&);
+        liquidViscosity(const liquidViscosity&);
 
         //- Disallow default bitwise assignment
-        void operator=(const surfaceShearForce&);
+        void operator=(const liquidViscosity&);
+
+
+protected:
+
+    // Protected data
+
 
 
 public:
 
     //- Runtime type information
-    TypeName("surfaceShear");
+    TypeName("liquid");
 
 
     // Constructors
 
         //- Construct from surface film model
-        surfaceShearForce
+        liquidViscosity
         (
             const surfaceFilmModel& owner,
-            const dictionary& dict
+            const dictionary& dict,
+            volScalarField& mu
         );
 
 
     //- Destructor
-    virtual ~surfaceShearForce();
+    virtual ~liquidViscosity();
 
 
     // Member Functions
@@ -97,7 +97,11 @@ public:
         // Evolution
 
             //- Correct
-            virtual tmp<fvVectorMatrix> correct(volVectorField& U);
+            virtual void correct
+            (
+                const volScalarField& p,
+                const volScalarField& T
+            );
 };
 
 

src/regionModels/surfaceFilmModels/submodels/thermo/filmViscosityModel/thixotropicViscosity/thixotropicViscosity.C

@@ -0,0 +1,181 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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 "thixotropicViscosity.H"
+#include "kinematicSingleLayer.H"
+#include "addToRunTimeSelectionTable.H"
+
+#include "fvmDdt.H"
+#include "fvmDiv.H"
+#include "fvcDiv.H"
+#include "fvmSup.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace regionModels
+{
+namespace surfaceFilmModels
+{
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+defineTypeNameAndDebug(thixotropicViscosity, 0);
+
+addToRunTimeSelectionTable
+(
+    filmViscosityModel,
+    thixotropicViscosity,
+    dictionary
+);
+
+
+// * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
+
+void thixotropicViscosity::updateMu()
+{
+    const kinematicSingleLayer& film = filmType<kinematicSingleLayer>();
+
+    // blend based on mass fraction of added- to existing film mass
+    const dimensionedScalar m0("zero", dimMass, 0.0);
+    const dimensionedScalar mSMALL("SMALL", dimMass, ROOTVSMALL);
+    const volScalarField deltaMass("deltaMass", max(m0, film.deltaMass()));
+    const volScalarField filmMass("filmMass", film.netMass() + mSMALL);
+
+    // weighting field to blend new and existing mass contributions
+    const volScalarField w("w", max(0.0, min(1.0, deltaMass/filmMass)));
+
+    // evaluate thixotropic viscosity
+    volScalarField muThx("muThx", muInf_/(sqr(1.0 - K_*lambda_) + ROOTVSMALL));
+
+    // set new viscosity based on weight field
+    mu_ = w*muInf_ + (1.0 - w)*muThx;
+    mu_.correctBoundaryConditions();
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+thixotropicViscosity::thixotropicViscosity
+(
+    const surfaceFilmModel& owner,
+    const dictionary& dict,
+    volScalarField& mu
+)
+:
+    filmViscosityModel(typeName, owner, dict, mu),
+    a_(coeffs().lookup("a")),
+    b_(coeffs().lookup("b")),
+    c_(coeffs().lookup("c")),
+    d_(coeffs().lookup("d")),
+    mu0_(coeffs().lookup("mu0")),
+    muInf_(coeffs().lookup("muInf")),
+    K_(1.0 - Foam::sqrt(muInf_/mu0_)),
+    lambda_
+    (
+        IOobject
+        (
+            typeName + ":lambda",
+            owner.regionMesh().time().timeName(),
+            owner.regionMesh(),
+            IOobject::MUST_READ_IF_MODIFIED,
+            IOobject::AUTO_WRITE
+        ),
+        owner.regionMesh()
+    )
+{
+    lambda_.min(1.0);
+    lambda_.max(0.0);
+
+    // initialise viscosity to inf value
+    // - cannot call updateMu() since this calls film.netMass() which
+    //   cannot be evaluated yet (still in construction)
+    mu_ = muInf_;
+    mu_.correctBoundaryConditions();
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+thixotropicViscosity::~thixotropicViscosity()
+{}
+
+
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+
+void thixotropicViscosity::correct
+(
+    const volScalarField& p,
+    const volScalarField& T
+)
+{
+    const kinematicSingleLayer& film = filmType<kinematicSingleLayer>();
+
+    // references to film fields
+    const volVectorField& U = film.U();
+    const volVectorField& Uw = film.Uw();
+    const volScalarField& delta = film.delta();
+    const volScalarField& deltaRho = film.deltaRho();
+    const surfaceScalarField& phi = film.phi();
+
+    // gamma-dot (shear rate) raised to the power d
+    volScalarField gDotPowD
+    (
+        "gDotPowD",
+        pow(mag(U - Uw)/(delta + film.deltaSmall()), d_)
+    );
+
+    dimensionedScalar c0("SMALL", dimMass/sqr(dimLength)/dimTime, SMALL);
+    volScalarField coeff(-deltaRho*c_*gDotPowD + c0);
+
+    fvScalarMatrix lambdaEqn
+    (
+        fvm::ddt(deltaRho, lambda_)
+      + fvm::div(phi, lambda_)
+      - fvm::Sp(fvc::div(phi), lambda_)
+      ==
+        deltaRho*a_*pow((1.0 - lambda_), b_)
+      + fvm::SuSp(coeff, lambda_)
+    );
+
+    lambdaEqn.relax();
+
+    lambdaEqn.solve();
+
+    lambda_.min(1.0);
+    lambda_.max(0.0);
+
+    updateMu();
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace surfaceFilmModels
+} // End namespace regionModels
+} // End namespace Foam
+
+// ************************************************************************* //

src/regionModels/surfaceFilmModels/submodels/thermo/filmViscosityModel/thixotropicViscosity/thixotropicViscosity.H

@@ -0,0 +1,205 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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::thixotropicViscosity
+
+Description
+    Thixotropic viscosity model based on the evolution of the structural
+    parameter \f$ \lambda \f$:
+
+        \f[
+            \lambda = a(1 - \lambda)^b - c \lambda \dot{\gamma}^d
+        \f]
+
+    The viscosity is then calculated using the expression
+
+        \f[
+            \mu = \frac{\mu_{inf}}{{1 - K \lambda}^2}
+        \f]
+
+    Where the parameter K is given by:
+
+        \f[
+            K = 1 - \frac{\mu_{\inf}}{\mu_{0}}^0.5
+        \f]
+
+    Here:
+    \vartable
+        \lambda         | structural parameter
+        a               | model coefficient
+        b               | model coefficient
+        c               | model coefficient
+        d               | model coefficient
+        \dot{\gamma}    | stress rate [1/s]
+        \mu_{0}         | limiting viscosity when \f$ \lambda = 1 \f$
+        \mu_{\inf}      | limiting viscosity when \f$ \lambda = 0 \f$
+    \endvartable
+
+
+    Reference:
+    \verbatim
+        Barnes H A, 1997.  Thixotropy - a review.  J. Non-Newtonian Fluid
+        Mech 70, pp 1-33
+    \endverbatim
+
+SourceFiles
+    thixotropicViscosity.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef thixotropicViscosity_H
+#define thixotropicViscosity_H
+
+#include "filmViscosityModel.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace regionModels
+{
+namespace surfaceFilmModels
+{
+
+/*---------------------------------------------------------------------------*\
+                    Class thixotropicViscosity Declaration
+\*---------------------------------------------------------------------------*/
+
+class thixotropicViscosity
+:
+    public filmViscosityModel
+{
+private:
+
+    // Private member functions
+
+        //- Disallow default bitwise copy construct
+        thixotropicViscosity(const thixotropicViscosity&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const thixotropicViscosity&);
+
+
+protected:
+
+    // Protected data
+
+/*
+        //- Model `a' coefficient
+        scalar a_;
+
+        //- Model `b' coefficient
+        scalar b_;
+
+        //- Model `c' coefficient
+        scalar c_;
+
+        //- Model `d' coefficient
+        scalar d_;
+
+        //- Limiting viscosity when lambda = 1
+        scalar mu0_;
+
+        //- Limiting viscosity when lambda = 0
+        scalar muInf_;
+*/
+
+        //- Model `a' coefficient
+        dimensionedScalar a_;
+
+        //- Model `b' coefficient
+        dimensionedScalar b_;
+
+        //- Model `c' coefficient
+        dimensionedScalar c_;
+
+        //- Model `d' coefficient
+        dimensionedScalar d_;
+
+        //- Limiting viscosity when lambda = 1
+        dimensionedScalar mu0_;
+
+        //- Limiting viscosity when lambda = 0
+        dimensionedScalar muInf_;
+
+        //- Model coeffiicient
+        dimensionedScalar K_;
+
+        //- Structural parameter
+        //  0 = freestream value (most liquid)
+        //  1 = fully built (most solid)
+        volScalarField lambda_;
+
+
+    // Protected Member Functions
+
+        //- Update the viscosity
+        void updateMu();
+
+
+public:
+
+    //- Runtime type information
+    TypeName("thixotropic");
+
+
+    // Constructors
+
+        //- Construct from surface film model
+        thixotropicViscosity
+        (
+            const surfaceFilmModel& owner,
+            const dictionary& dict,
+            volScalarField& mu
+        );
+
+
+    //- Destructor
+    virtual ~thixotropicViscosity();
+
+
+    // Member Functions
+
+        // Evolution
+
+            //- Correct
+            virtual void correct
+            (
+                const volScalarField& p,
+                const volScalarField& T
+            );
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace surfaceFilmModels
+} // End namespace regionModels
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/regionModels/surfaceFilmModels/submodels/thermo/phaseChangeModel/solidification/solidification.C

@@ -0,0 +1,137 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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 "solidification.H"
+#include "addToRunTimeSelectionTable.H"
+#include "thermoSingleLayer.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace regionModels
+{
+namespace surfaceFilmModels
+{
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+defineTypeNameAndDebug(solidification, 0);
+
+addToRunTimeSelectionTable
+(
+    phaseChangeModel,
+    solidification,
+    dictionary
+);
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+solidification::solidification
+(
+    const surfaceFilmModel& owner,
+    const dictionary& dict
+)
+:
+    phaseChangeModel(typeName, owner, dict),
+    T0_(readScalar(coeffs_.lookup("T0"))),
+    L_(readScalar(coeffs_.lookup("L"))),
+    alpha_(readScalar(coeffs_.lookup("alpha"))),
+    mass_
+    (
+        IOobject
+        (
+            typeName + ":mass",
+            owner.regionMesh().time().timeName(),
+            owner.regionMesh(),
+            IOobject::READ_IF_PRESENT,
+            IOobject::AUTO_WRITE
+        ),
+        owner.regionMesh(),
+        dimensionedScalar("zero", dimMass, 0.0),
+        zeroGradientFvPatchScalarField::typeName
+    ),
+    thickness_
+    (
+        IOobject
+        (
+            typeName + ":thickness",
+            owner.regionMesh().time().timeName(),
+            owner.regionMesh(),
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        owner.regionMesh(),
+        dimensionedScalar("zero", dimLength, 0.0),
+        zeroGradientFvPatchScalarField::typeName
+    )
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+solidification::~solidification()
+{}
+
+
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+
+void solidification::correctModel
+(
+    const scalar dt,
+    scalarField& availableMass,
+    scalarField& dMass,
+    scalarField& dEnergy
+)
+{
+    const thermoSingleLayer& film = filmType<thermoSingleLayer>();
+
+    const scalarField& T = film.T();
+    const scalarField& alpha = film.alpha();
+
+    forAll(alpha, cellI)
+    {
+        if (alpha[cellI] > 0.5)
+        {
+            if (T[cellI] > T0_)
+            {
+                mass_[cellI] += alpha_*availableMass[cellI];
+                dMass[cellI] += alpha_*availableMass[cellI];
+                dEnergy[cellI] += alpha_*availableMass[cellI]*L_;
+            }
+        }
+    }
+
+    thickness_ = mass_/film.magSf()/film.rho();
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace surfaceFilmModels
+} // End namespace regionModels
+} // End namespace Foam
+
+// ************************************************************************* //

src/regionModels/surfaceFilmModels/submodels/thermo/phaseChangeModel/solidification/solidification.H

@@ -0,0 +1,130 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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::solidification
+
+Description
+    Solidification phase change model where all film mass is converted when
+    the local temperature > activation temperature.
+
+SourceFiles
+    solidification.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef solidification_H
+#define solidification_H
+
+#include "phaseChangeModel.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace regionModels
+{
+namespace surfaceFilmModels
+{
+
+/*---------------------------------------------------------------------------*\
+                       Class solidification Declaration
+\*---------------------------------------------------------------------------*/
+
+class solidification
+:
+    public phaseChangeModel
+{
+private:
+
+    // Private member functions
+
+        //- Disallow default bitwise copy construct
+        solidification(const solidification&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const solidification&);
+
+
+protected:
+
+    // Protected data
+
+        //- Temperature at which solidification starts
+        scalar T0_;
+
+        //- Latent heat of solidification [J/kg]
+        scalar L_;
+
+        //- Under-relaxation parameter for solidification process (0-1)
+        scalar alpha_;
+
+        //- Accumulated solid mass [kg]
+        volScalarField mass_;
+
+        //- Accumulated solid thickness [m]
+        volScalarField thickness_;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("solidification");
+
+
+    // Constructors
+
+        //- Construct from surface film model
+        solidification(const surfaceFilmModel& owner, const dictionary& dict);
+
+
+    //- Destructor
+    virtual ~solidification();
+
+
+    // Member Functions
+
+        // Evolution
+
+            //- Correct
+            virtual void correctModel
+            (
+                const scalar dt,
+                scalarField& availableMass,
+                scalarField& dMass,
+                scalarField& dEnergy
+            );
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace surfaceFilmModels
+} // End namespace regionModels
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/regionModels/surfaceFilmModels/submodels/thermo/phaseChangeModel/standardPhaseChange/standardPhaseChange.C

@@ -97,13 +97,12 @@ void standardPhaseChange::correctModel
     scalarField& dEnergy
 )
 {
-    const thermoSingleLayer& film = refCast<const thermoSingleLayer>(owner_);
+    const thermoSingleLayer& film = filmType<thermoSingleLayer>();
 
     // set local thermo properties
     const SLGThermo& thermo = film.thermo();
-    const label liqId = film.liquidId();
+    const filmThermoModel& filmThermo = film.filmThermo();
-    const liquidProperties& liq = thermo.liquids().properties()[liqId];
+    const label vapId = thermo.carrierId(filmThermo.name());
-    const label vapId = thermo.carrierId(thermo.liquids().components()[liqId]);
 
     // retrieve fields from film model
     const scalarField& delta = film.delta();
@@ -128,22 +127,22 @@ void standardPhaseChange::correctModel
             const scalar pc = pInf[cellI];
 
             // calculate the boiling temperature
-            const scalar Tb = liq.pvInvert(pc);
+            const scalar Tb = filmThermo.Tb(pc);
 
             // local temperature - impose lower limit of 200 K for stability
             const scalar Tloc = min(TbFactor_*Tb, max(200.0, T[cellI]));
 
             // saturation pressure [Pa]
-            const scalar pSat = liq.pv(pc, Tloc);
+            const scalar pSat = filmThermo.pv(pc, Tloc);
 
             // latent heat [J/kg]
-            const scalar hVap = liq.hl(pc, Tloc);
+            const scalar hVap = filmThermo.hl(pc, Tloc);
 
             // calculate mass transfer
             if (pSat >= 0.95*pc)
             {
                 // boiling
-                const scalar Cp = liq.Cp(pc, Tloc);
+                const scalar Cp = filmThermo.Cp(pc, Tloc);
                 const scalar Tcorr = max(0.0, T[cellI] - Tb);
                 const scalar qCorr = limMass[cellI]*Cp*(Tcorr);
                 dMass[cellI] = qCorr/hVap;
@@ -163,13 +162,13 @@ void standardPhaseChange::correctModel
                 const scalar Wvap = thermo.carrier().W(vapId);
 
                 // molecular weight of liquid [kg/kmol]
-                const scalar Wliq = liq.W();
+                const scalar Wliq = filmThermo.W();
 
                 // vapour mass fraction at interface
                 const scalar Ys = Wliq*pSat/(Wliq*pSat + Wvap*(pc - pSat));
 
                 // vapour diffusivity [m2/s]
-                const scalar Dab = liq.D(pc, Tloc);
+                const scalar Dab = filmThermo.D(pc, Tloc);
 
                 // Schmidt number
                 const scalar Sc = muInfc/(rhoInfc*(Dab + ROOTVSMALL));

src/regionModels/surfaceFilmModels/surfaceFilmModel/surfaceFilmModel.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
@@ -26,30 +26,6 @@ License
 #include "surfaceFilmModel.H"
 #include "fvMesh.H"
 
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-    template<>
-    const char* NamedEnum
-    <
-        regionModels::surfaceFilmModels::surfaceFilmModel::thermoModelType,
-        2
-    >::names[] =
-    {
-        "constant",
-        "singleComponent"
-    };
-}
-
-const Foam::NamedEnum
-<
-    Foam::regionModels::surfaceFilmModels::surfaceFilmModel::thermoModelType,
-    2
->
-Foam::regionModels::surfaceFilmModels::surfaceFilmModel::thermoModelTypeNames_;
-
-
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 namespace Foam
@@ -70,8 +46,6 @@ bool surfaceFilmModel::read()
 {
     if (singleLayerRegion::read())
     {
-        thermoModel_ =
-            thermoModelTypeNames_.read(coeffs_.lookup("thermoModel"));
         return true;
     }
     else
@@ -91,8 +65,7 @@ surfaceFilmModel::surfaceFilmModel
 )
 :
     singleLayerRegion(mesh, "surfaceFilm", modelType),
-    g_(g),
+    g_(g)
-    thermoModel_(tmConstant)
 {
     if (active_)
     {

src/regionModels/surfaceFilmModels/surfaceFilmModel/surfaceFilmModel.H

@@ -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
@@ -63,20 +63,6 @@ class surfaceFilmModel
 :
     public singleLayerRegion
 {
-public:
-
-    // Data types
-
-        //- Enumeration listing the possible thermo types
-        enum thermoModelType
-        {
-            tmConstant,
-            tmSingleComponent
-        };
-
-        //- Named enumeration for the thermoType
-        static const NamedEnum<thermoModelType, 2> thermoModelTypeNames_;
-
 
 private:
 
@@ -97,9 +83,6 @@ protected:
         //- Acceleration due to gravity [m/s2]
         const dimensionedVector& g_;
 
-        //- Thermo type
-        thermoModelType thermoModel_;
-
 
     // Protected member functions
 
@@ -160,9 +143,6 @@ public:
             //- Return the accleration due to gravity
             inline const dimensionedVector& g() const;
 
-            //- Return the thermo type
-            inline const thermoModelType& thermoModel() const;
-
             //- External hook to add sources to the film
             virtual void addSources
             (

src/regionModels/surfaceFilmModels/surfaceFilmModel/surfaceFilmModelI.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
@@ -42,13 +42,6 @@ inline const Foam::dimensionedVector& surfaceFilmModel::g() const
 }
 
 
-inline const surfaceFilmModel::thermoModelType&
-surfaceFilmModel::thermoModel() const
-{
-    return thermoModel_;
-}
-
-
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 } // End namespace surfaceFilmModels

src/regionModels/surfaceFilmModels/thermoSingleLayer/thermoSingleLayer.C

@@ -33,6 +33,8 @@ License
 #include "mapDistribute.H"
 
 // Sub-models
+#include "filmThermoModel.H"
+#include "filmViscosityModel.H"
 #include "heatTransferModel.H"
 #include "phaseChangeModel.H"
 #include "filmRadiationModel.H"
@@ -59,7 +61,11 @@ wordList thermoSingleLayer::hsBoundaryTypes()
     wordList bTypes(T_.boundaryField().types());
     forAll(bTypes, patchI)
     {
-        if (bTypes[patchI] == mappedFieldFvPatchField<scalar>::typeName)
+        if
+        (
+            T_.boundaryField()[patchI].fixesValue()
+         || bTypes[patchI] == mappedFieldFvPatchField<scalar>::typeName
+        )
         {
             bTypes[patchI] = fixedValueFvPatchField<scalar>::typeName;
         }
@@ -93,53 +99,10 @@ void thermoSingleLayer::resetPrimaryRegionSourceTerms()
 
 void thermoSingleLayer::correctThermoFields()
 {
-    switch (thermoModel_)
+    rho_ == filmThermo_->rho();
-    {
+    sigma_ == filmThermo_->sigma();
-        case tmConstant:
+    Cp_ == filmThermo_->Cp();
-        {
+    kappa_ == filmThermo_->kappa();
-            const dictionary&
-                constDict(coeffs_.subDict("constantThermoCoeffs"));
-            rho_ == dimensionedScalar(constDict.lookup("rho0"));
-            mu_ == dimensionedScalar(constDict.lookup("mu0"));
-            sigma_ == dimensionedScalar(constDict.lookup("sigma0"));
-            Cp_ == dimensionedScalar(constDict.lookup("Cp0"));
-            kappa_ == dimensionedScalar(constDict.lookup("kappa0"));
-
-            break;
-        }
-        case tmSingleComponent:
-        {
-            const liquidProperties& liq =
-                thermo_.liquids().properties()[liquidId_];
-
-            forAll(rho_, cellI)
-            {
-                const scalar T = T_[cellI];
-                const scalar p = pPrimary_[cellI];
-                rho_[cellI] = liq.rho(p, T);
-                mu_[cellI] = liq.mu(p, T);
-                sigma_[cellI] = liq.sigma(p, T);
-                Cp_[cellI] = liq.Cp(p, T);
-                kappa_[cellI] = liq.K(p, T);
-            }
-
-            rho_.correctBoundaryConditions();
-            mu_.correctBoundaryConditions();
-            sigma_.correctBoundaryConditions();
-            Cp_.correctBoundaryConditions();
-            kappa_.correctBoundaryConditions();
-
-            break;
-        }
-        default:
-        {
-            FatalErrorIn
-            (
-                "void thermoSingleLayer::"
-                "correctThermoFields()"
-            )   << "Unknown thermoType enumeration" << abort(FatalError);
-        }
-    }
 }
 
 
@@ -286,7 +249,7 @@ void thermoSingleLayer::updateSubmodels()
     rhoSp_ += primaryMassPCTrans_/magSf()/time().deltaT();
 
     // Vapour recoil pressure
-    pSp_ -= sqr(primaryMassPCTrans_/magSf()/time_.deltaT())/2.0/rhoPrimary_;
+    pSp_ -= sqr(primaryMassPCTrans_/magSf()/time().deltaT())/2.0/rhoPrimary_;
 }
 
 
@@ -294,20 +257,22 @@ tmp<fvScalarMatrix> thermoSingleLayer::q(volScalarField& hs) const
 {
     dimensionedScalar Tstd("Tstd", dimTemperature, 298.15);
 
-    volScalarField htcst = htcs_->h()();
+    volScalarField htcst(htcs_->h());
-    volScalarField htcwt = htcw_->h()();
+    volScalarField htcwt(htcw_->h());
+
     forAll(alpha_, i)
     {
         htcst[i] *= max(alpha_[i], ROOTVSMALL);
         htcwt[i] *= max(alpha_[i], ROOTVSMALL);
     }
+
     htcst.correctBoundaryConditions();
     htcwt.correctBoundaryConditions();
 
     return
     (
       - fvm::Sp(htcst/Cp_, hs) - htcst*(Tstd - TPrimary_)
-      - fvm::Sp(htcwt/Cp_, hs) -htcwt*(Tstd - Tw_)
+      - fvm::Sp(htcwt/Cp_, hs) - htcwt*(Tstd - Tw_)
     );
 }
 
@@ -334,6 +299,9 @@ void thermoSingleLayer::solveEnergy()
     );
 
     correctThermoFields();
+
+    // evaluate viscosity from user-model
+    viscosity_->correct(pPrimary_, T_);
 }
 
 
@@ -349,7 +317,6 @@ thermoSingleLayer::thermoSingleLayer
 :
     kinematicSingleLayer(modelType, mesh, g, false),
     thermo_(mesh.lookupObject<SLGThermo>("SLGThermo")),
-    liquidId_(thermo_.liquidId(coeffs_.lookup("liquid"))),
     Cp_
     (
         IOobject
@@ -517,6 +484,7 @@ thermoSingleLayer::thermoSingleLayer
 
     YPrimary_(),
 
+    viscosity_(filmViscosityModel::New(*this, coeffs(), mu_)),
     htcs_
     (
         heatTransferModel::New(*this, coeffs().subDict("upperSurfaceModels"))
@@ -576,12 +544,18 @@ thermoSingleLayer::thermoSingleLayer
     if (readFields)
     {
         transferPrimaryRegionThermoFields();
+
+        correctAlpha();
+
         correctThermoFields();
 
         // Update derived fields
         hs_ == hs(T_);
         deltaRho_ == delta_*rho_;
         phi_ = fvc::interpolate(deltaRho_*U_) & regionMesh().Sf();
+
+        // evaluate viscosity from user-model
+        viscosity_->correct(pPrimary_, T_);
     }
 }
 
@@ -802,8 +776,7 @@ tmp<DimensionedField<scalar, volMesh> > thermoSingleLayer::Srho
     const label i
 ) const
 {
-    const label vapId =
+    const label vapId = thermo_.carrierId(filmThermo_->name());
-        thermo_.carrierId(thermo_.liquids().components()[liquidId_]);
 
     tmp<DimensionedField<scalar, volMesh> > tSrho
     (

src/regionModels/surfaceFilmModels/thermoSingleLayer/thermoSingleLayer.H

@@ -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
@@ -52,6 +52,7 @@ namespace surfaceFilmModels
 {
 
 // Forward declaration of classes
+class filmViscosityModel;
 class heatTransferModel;
 class phaseChangeModel;
 class filmRadiationModel;
@@ -87,11 +88,6 @@ protected:
             //- Reference to the SLGThermo
             const SLGThermo& thermo_;
 
-            // Single component
-
-                //- Id of component in thermo database
-                label liquidId_;
-
 
             // Fields
 
@@ -170,6 +166,9 @@ protected:
 
         // Sub-models
 
+            //- Viscosity model
+            autoPtr<filmViscosityModel> viscosity_;
+
             //- Heat transfer coefficient bewteen film surface and primary
             //  region [W/m2/K]
             autoPtr<heatTransferModel> htcs_;
@@ -261,11 +260,6 @@ public:
             //- Return const reference to the SLGThermo object
             inline const SLGThermo& thermo() const;
 
-            // Single component
-
-                //- Return the Id of component in thermo database
-                inline label liquidId() const;
-
 
             // Fields
 

src/regionModels/surfaceFilmModels/thermoSingleLayer/thermoSingleLayerI.H

@@ -24,8 +24,8 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "thermoSingleLayer.H"
-#include "heatTransferModel.H"
 #include "filmRadiationModel.H"
+#include "heatTransferModel.H"
 #include "phaseChangeModel.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -45,12 +45,6 @@ inline const SLGThermo& thermoSingleLayer::thermo() const
 }
 
 
-inline label thermoSingleLayer::liquidId() const
-{
-    return liquidId_;
-}
-
-
 inline tmp<scalarField> thermoSingleLayer::hs
 (
     const scalarField& T,

src/thermophysicalModels/basic/basicThermo/basicThermo.H

@@ -405,32 +405,35 @@ public:
             //- Thermal diffusivity for temperature of mixture [J/m/s/K]
             virtual tmp<volScalarField> kappa() const = 0;
 
-            //- Thermal diffusivity of mixture for patch [J/m/s/K]
+            //- Thermal diffusivity for temperature
+            //  of mixture for patch [J/m/s/K]
             virtual tmp<scalarField> kappa
             (
                 const label patchi
             ) const = 0;
 
-            //- Effective thermal diffusivity of mixture [J/m/s/K]
+            //- Effective thermal diffusivity for temperature
+            //  of mixture [J/m/s/K]
             virtual tmp<volScalarField> kappaEff
             (
                 const volScalarField&
             ) const = 0;
 
-            //- Effective thermal diffusivity of mixture for patch [J/m/s/K]
+            //- Effective thermal diffusivity for temperature
+            //  of mixture for patch [J/m/s/K]
             virtual tmp<scalarField> kappaEff
             (
                 const scalarField& alphat,
                 const label patchi
             ) const = 0;
 
-            //- Effective thermal diffusivity of mixture [J/m/s/K]
+            //- Effective thermal diffusivity of mixture [kg/m/s]
             virtual tmp<volScalarField> alphaEff
             (
                 const volScalarField& alphat
             ) const = 0;
 
-            //- Effective thermal diffusivity of mixture for patch [J/m/s/K]
+            //- Effective thermal diffusivity of mixture for patch [kg/m/s]
             virtual tmp<scalarField> alphaEff
             (
                 const scalarField& alphat,

src/thermophysicalModels/basic/heThermo/heThermo.H

@@ -259,30 +259,32 @@ public:
             //- Thermal diffusivity for temperature of mixture [J/m/s/K]
             virtual tmp<volScalarField> kappa() const;
 
-            //- Thermal diffusivity of mixture for patch [J/m/s/K]
+            //- Thermal diffusivity for temperature
+            //  of mixture for patch [J/m/s/K]
             virtual tmp<scalarField> kappa
             (
                 const label patchi
             ) const;
 
-            //- Effective thermal diffusivity of mixture [J/m/s/K]
+            //- Effective thermal diffusivity for temperature
+            //  of mixture [J/m/s/K]
             virtual tmp<volScalarField> kappaEff(const volScalarField&) const;
 
-            //- Effective thermal diffusivity of mixture for patch [J/m/s/K]
+            //- Effective thermal diffusivity for temperature
+            //  of mixture for patch [J/m/s/K]
             virtual tmp<scalarField> kappaEff
             (
                 const scalarField& alphat,
                 const label patchi
             ) const;
 
-
+            //- Effective thermal diffusivity of mixture [kg/m/s]
-            //- Effective thermal diffusivity of mixture [J/m/s/K]
             virtual tmp<volScalarField> alphaEff
             (
                 const volScalarField& alphat
             ) const;
 
-            //- Effective thermal diffusivity of mixture for patch [J/m/s/K]
+            //- Effective thermal diffusivity of mixture for patch [kg/m/s]
             virtual tmp<scalarField> alphaEff
             (
                 const scalarField& alphat,

src/thermophysicalModels/basic/psiThermo/hePsiThermo.H

@@ -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
@@ -91,7 +91,7 @@ public:
 
 } // End namespace Foam
 
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 #ifdef NoRepository
 #   include "hePsiThermo.C"

src/thermophysicalModels/basic/rhoThermo/heRhoThermo.H

@@ -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
@@ -92,7 +92,7 @@ public:
 
 } // End namespace Foam
 
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 #ifdef NoRepository
 #   include "heRhoThermo.C"

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;
     }
 

tutorials/incompressible/pimpleDyMFoam/propeller/Allclean

@@ -5,10 +5,10 @@ cd ${0%/*} || exit 1    # run from this directory
 . $WM_PROJECT_DIR/bin/tools/CleanFunctions
 
 # remove surface
-\rm -f constant/triSurface/propellerTip.obj.gz
+\rm -f constant/triSurface/propellerTip.obj.gz > /dev/null 2>&1
 
 cleanCase
 
-\rm -rf 0
+\rm -rf 0 > /dev/null 2>&1
-\rm -rf constant/extendedFeatureEdgeMesh/
+\rm -rf constant/extendedFeatureEdgeMesh > /dev/null 2>&1
-\rm -f constant/triSurface/*.eMesh*
+\rm -f constant/triSurface/*.eMesh* > /dev/null 2>&1

tutorials/incompressible/porousSimpleFoam/straightDuctImplicit/Allclean

@@ -4,10 +4,10 @@ cd ${0%/*} || exit 1    # run from this directory
 # Source tutorial clean functions
 . $WM_PROJECT_DIR/bin/tools/CleanFunctions
 
-rm *.obj
+rm *.obj > /dev/null 2>&1
-rm -r constant/extendedFeatureEdgeMesh
+rm -r constant/extendedFeatureEdgeMesh > /dev/null 2>&1
-rm constant/triSurface/boundaryAndFaceZones.eMesh
+rm constant/triSurface/boundaryAndFaceZones.eMesh > /dev/null 2>&1
-rm constant/polyMesh/boundary
+rm constant/polyMesh/boundary > /dev/null 2>&1
 
 cleanCase
 

tutorials/incompressible/simpleFoam/motorBike/Allclean

@@ -4,9 +4,9 @@
 . $WM_PROJECT_DIR/bin/tools/CleanFunctions
 
 # remove surface and features
-\rm -f constant/triSurface/motorBike.obj.gz
+\rm -f constant/triSurface/motorBike.obj.gz > /dev/null 2>&1
-\rm -rf constant/extendedFeatureEdgeMesh/
+\rm -rf constant/extendedFeatureEdgeMesh > /dev/null 2>&1
-\rm -f constant/triSurface/motorBike.eMesh
+\rm -f constant/triSurface/motorBike.eMesh > /dev/null 2>&1
 
 rm -rf 0 > /dev/null 2>&1
 

tutorials/lagrangian/reactingParcelFilmFoam/cylinder/constant/surfaceFilmProperties

@@ -23,10 +23,18 @@ active          true;
 
 thermoSingleLayerCoeffs
 {
-    thermoModel singleComponent;
+    filmThermoModel liquid;
-    liquid      H2O;
+
+    liquidCoeffs
+    {
+        useReferenceValues  no;
+        liquid      H2O;
+    }
+
+    filmViscosityModel liquid;
 
     deltaWet    1e-4;
+
     hydrophilic no;
 
     turbulence  laminar;

tutorials/lagrangian/reactingParcelFilmFoam/hotBoxes/constant/surfaceFilmProperties

@@ -24,11 +24,18 @@ active          true;
 
 thermoSingleLayerCoeffs
 {
-    thermoModel singleComponent; // constant
+    filmThermoModel liquid;
 
-    liquid      H2O;
+    liquidCoeffs
+    {
+        useReferenceValues no;
+        liquid      H2O;
+    }
+
+    filmViscosityModel liquid;
 
     deltaWet    1e-4;
+
     hydrophilic no;
 
     radiationModel  none;

tutorials/lagrangian/reactingParcelFilmFoam/rivuletPanel/constant/surfaceFilmProperties

@@ -23,12 +23,13 @@ active          true;
 
 kinematicSingleLayerCoeffs
 {
-    thermoModel constant;
+    filmThermoModel constant;
-    constantThermoCoeffs
+    constantCoeffs
     {
-        rho0        rho0 [1 -3 0 0 0] 1000;
+        specieName  water;
-        mu0         mu0 [1 -1 -1 0 0] 1e-3;
+        rho0        1000;
-        sigma0      sigma0 [1 0 -2 0 0] 0.07;
+        mu0         1e-3;
+        sigma0      0.07;
     }
 
     injectionModels ();

tutorials/lagrangian/reactingParcelFilmFoam/splashPanel/constant/surfaceFilmProperties

@@ -23,9 +23,15 @@ active          true;
 
 thermoSingleLayerCoeffs
 {
-    thermoModel singleComponent;
+    filmThermoModel liquid;
 
-    liquid      H2O;
+    liquidCoeffs
+    {
+        useReferenceValues  no;
+        liquid      H2O;
+    }
+
+    filmViscosityModel liquid;
 
     deltaWet    1e-4;
     hydrophilic no;

tutorials/multiphase/interDyMFoam/ras/mixerVesselAMI/Allclean

@@ -6,8 +6,8 @@ cd ${0%/*} || exit 1    # run from this directory
 
 cleanCase
 
-\rm -rf 0
+\rm -rf 0 > /dev/null 2>&1
-\rm -rf constant/extendedFeatureEdgeMesh/
+\rm -rf constant/extendedFeatureEdgeMesh > /dev/null 2>&1
-\rm -f constant/triSurface/*.eMesh*
+\rm -f constant/triSurface/*.eMesh* > /dev/null 2>&1
 
 # ----------------------------------------------------------------- end-of-file

tutorials/multiphase/interPhaseChangeDyMFoam/propeller/Allclean

@@ -9,6 +9,6 @@ cd ${0%/*} || exit 1    # run from this directory
 
 cleanCase
 
-\rm -rf 0
+\rm -rf 0 > /dev/null 2>&1
-\rm -rf constant/extendedFeatureEdgeMesh/
+\rm -rf constant/extendedFeatureEdgeMesh > /dev/null 2>&1
-\rm -f constant/triSurface/*.eMesh*
+\rm -f constant/triSurface/*.eMesh* > /dev/null 2>&1

tutorials/multiphase/twoPhaseEulerFoam/LES/bubbleColumn/constant/phaseProperties

@@ -36,6 +36,9 @@ water
     }
 }
 
+// Surface tension coefficient
+sigma           0.07;
+
 drag
 {
     air     SchillerNaumann;

tutorials/multiphase/twoPhaseEulerFoam/RAS/bubbleColumn/constant/phaseProperties

@@ -36,6 +36,9 @@ water
     }
 }
 
+// Surface tension coefficient
+sigma           0.07;
+
 drag
 {
     air     SchillerNaumann;

tutorials/multiphase/twoPhaseEulerFoam/RAS/fluidisedBed/constant/phaseProperties

@@ -35,6 +35,9 @@ air
     }
 }
 
+// Surface tension coefficient
+sigma           0;
+
 drag
 {
     particles   GidaspowErgunWenYu;

tutorials/multiphase/twoPhaseEulerFoam/laminar/bubbleColumn/constant/phaseProperties

@@ -36,6 +36,9 @@ water
     }
 }
 
+// Surface tension coefficient
+sigma           0.07;
+
 drag
 {
     air     SchillerNaumann;

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;
+    }
+}
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

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;
+    }
+}
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

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;
+    }
+}
+
+// ************************************************************************* //

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);
+    }
+}
+
+// ************************************************************************* //

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);
+    }
+}
+
+// ************************************************************************* //

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;
+    }
+}
+
+// ************************************************************************* //

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;
+    }
+}
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

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;
+    }
+}
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

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);
+
+
+// ************************************************************************* //