Merge branch 'master' of /home/noisy3/OpenFOAM/OpenFOAM-dev

applications/solvers/lagrangian/reactingParcelFilmFoam/UEqn.H

@@ -4,13 +4,23 @@
       + fvm::div(phi, U)
       + turbulence->divDevRhoReff(U)
      ==
-        rho.dimensionedInternalField()*g
+        parcels.SU(U)
-      + parcels.SU(U)
     );
 
     UEqn.relax();
 
     if (momentumPredictor)
     {
-        solve(UEqn == -fvc::grad(p));
+        solve
+        (
+            UEqn
+          ==
+            fvc::reconstruct
+            (
+                (
+                  - ghf*fvc::snGrad(rho)
+                  - fvc::snGrad(p_rgh)
+                )*mesh.magSf()
+            )
+        );
     }

applications/solvers/lagrangian/reactingParcelFilmFoam/YEqn.H

@@ -5,7 +5,7 @@ tmp<fv::convectionScheme<scalar> > mvConvection
         mesh,
         fields,
         phi,
-        mesh.divScheme("div(phi,Yi_h)")
+        mesh.divScheme("div(phi,Yi_hs)")
     )
 );
 
@@ -23,7 +23,7 @@ tmp<fv::convectionScheme<scalar> > mvConvection
             (
                 fvm::ddt(rho, Yi)
               + mvConvection->fvmDiv(phi, Yi)
-              - fvm::laplacian(turbulence->muEff(), Yi)
+              - fvm::laplacian(turbulence->alphaEff(), Yi)
               ==
                 parcels.SYi(i, Yi)
               + surfaceFilm.Srho(i)

applications/solvers/lagrangian/reactingParcelFilmFoam/chemistry.H

@@ -1,3 +1,4 @@
+if (chemistry.chemistry())
 {
     Info << "Solving chemistry" << endl;
 

applications/solvers/lagrangian/reactingParcelFilmFoam/createFields.H

@@ -15,11 +15,6 @@
 
     const word inertSpecie(thermo.lookup("inertSpecie"));
 
-    volScalarField& p = thermo.p();
-    volScalarField& hs = thermo.hs();
-    const volScalarField& T = thermo.T();
-    const volScalarField& psi = thermo.psi();
-
     Info<< "Creating field rho\n" << endl;
     volScalarField rho
     (
@@ -34,6 +29,11 @@
         thermo.rho()
     );
 
+    volScalarField& p = thermo.p();
+    volScalarField& hs = thermo.hs();
+    const volScalarField& T = thermo.T();
+    const volScalarField& psi = thermo.psi();
+
     Info<< "\nReading field U\n" << endl;
     volVectorField U
     (
@@ -84,6 +84,28 @@
         fvc::DDt(surfaceScalarField("phiU", phi/fvc::interpolate(rho)), p)
     );
 
+
+    Info<< "Calculating field g.h\n" << endl;
+    volScalarField gh("gh", g & mesh.C());
+
+    surfaceScalarField ghf("gh", g & mesh.Cf());
+
+    volScalarField p_rgh
+    (
+        IOobject
+        (
+            "p_rgh",
+            runTime.timeName(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh
+    );
+
+    // Force p_rgh to be consistent with p
+    p_rgh = p - rho*gh;
+
     multivariateSurfaceInterpolationScheme<scalar>::fieldTable fields;
 
     forAll(Y, i)

applications/solvers/lagrangian/reactingParcelFilmFoam/hsEqn.H

@@ -19,4 +19,6 @@
     thermo.correct();
 
     radiation->correct();
+
+    Info<< "min/max(T) = " << min(T).value() << ", " << max(T).value() << endl;
 }

applications/solvers/lagrangian/reactingParcelFilmFoam/pEqn.H

@@ -1,74 +1,58 @@
 rho = thermo.rho();
 
 volScalarField rAU(1.0/UEqn.A());
+surfaceScalarField rhorAUf(rAU.name(), fvc::interpolate(rho*rAU));
 U = rAU*UEqn.H();
 
-if (transonic)
+surfaceScalarField phiU
-{
-    surfaceScalarField phid
 (
-        "phid",
+    fvc::interpolate(rho)
-        fvc::interpolate(psi)
    *(
         (fvc::interpolate(U) & mesh.Sf())
       + fvc::ddtPhiCorr(rAU, rho, U, phi)
     )
 );
 
-    for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
+phi = phiU - rhorAUf*ghf*fvc::snGrad(rho)*mesh.magSf();
-    {
-        fvScalarMatrix pEqn
-        (
-            fvm::ddt(psi, p)
-          + fvm::div(phid, p)
-          - fvm::laplacian(rho*rAU, p)
-         ==
-            parcels.Srho()
-          + surfaceFilm.Srho()
-        );
-
-        pEqn.solve();
-
-        if (nonOrth == nNonOrthCorr)
-        {
-            phi == pEqn.flux();
-        }
-    }
-}
-else
-{
-    phi =
-        fvc::interpolate(rho)
-       *(
-            (fvc::interpolate(U) & mesh.Sf())
-          + fvc::ddtPhiCorr(rAU, rho, U, phi)
-        );
 
 for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
 {
-        fvScalarMatrix pEqn
+    fvScalarMatrix p_rghEqn
     (
-            fvm::ddt(psi, p)
+        fvc::ddt(psi, rho)*gh
       + fvc::div(phi)
-          - fvm::laplacian(rho*rAU, p)
+      + fvm::ddt(psi, p_rgh)
+      - fvm::laplacian(rhorAUf, p_rgh)
      ==
         parcels.Srho()
       + surfaceFilm.Srho()
     );
 
-        pEqn.solve();
+    p_rghEqn.solve
+    (
+        mesh.solver
+        (
+            p_rgh.select
+            (
+                pimpleCorr.finalIter()
+             && corr == nCorr-1
+             && nonOrth == nNonOrthCorr
+            )
+        )
+    );
 
     if (nonOrth == nNonOrthCorr)
     {
-            phi += pEqn.flux();
+        phi += p_rghEqn.flux();
-        }
     }
 }
 
+p = p_rgh + rho*gh;
+
 #include "rhoEqn.H"
 #include "compressibleContinuityErrs.H"
 
-U -= rAU*fvc::grad(p);
+U += rAU*fvc::reconstruct((phi - phiU)/rhorAUf);
 U.correctBoundaryConditions();
 
 DpDt = fvc::DDt(surfaceScalarField("phiU", phi/fvc::interpolate(rho)), p);

applications/solvers/lagrangian/reactingParcelFilmFoam/reactingParcelFilmFoam.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2009-2010 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2009-2011 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -39,6 +39,7 @@ Description
 #include "chemistrySolver.H"
 #include "radiationModel.H"
 #include "SLGThermo.H"
+#include "pimpleLoop.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -66,8 +67,9 @@ int main(int argc, char *argv[])
     while (runTime.run())
     {
         #include "readTimeControls.H"
-        #include "readPISOControls.H"
+        #include "readPIMPLEControls.H"
         #include "compressibleCourantNo.H"
+        #include "setMultiRegionDeltaT.H"
         #include "setDeltaT.H"
 
         runTime++;
@@ -84,20 +86,22 @@ int main(int argc, char *argv[])
             #include "rhoEqn.H"
 
             // --- PIMPLE loop
-            for (int ocorr=1; ocorr<=nOuterCorr; ocorr++)
+            for
+            (
+                pimpleLoop pimpleCorr(mesh, nOuterCorr);
+                pimpleCorr.loop();
+                pimpleCorr++
+            )
             {
                 #include "UEqn.H"
                 #include "YEqn.H"
+                #include "hsEqn.H"
 
                 // --- PISO loop
                 for (int corr=1; corr<=nCorr; corr++)
                 {
-                    #include "hsEqn.H"
                     #include "pEqn.H"
                 }
-
-                Info<< "T gas min/max   = " << min(T).value() << ", "
-                    << max(T).value() << endl;
             }
 
             turbulence->correct();

applications/solvers/lagrangian/reactingParcelFilmFoam/setMultiRegionDeltaT.H

@@ -0,0 +1,57 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2010-2011 OpenCFD Ltd.
+     \\/     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/>.
+
+Global
+    setMultiRegionDeltaT
+
+Description
+    Reset the timestep to maintain a constant maximum Courant numbers.
+    Reduction of time-step is immediate, but increase is damped to avoid
+    unstable oscillations.
+
+\*---------------------------------------------------------------------------*/
+
+if (adjustTimeStep)
+{
+    if (CoNum == -GREAT)
+    {
+        CoNum = SMALL;
+    }
+
+    const scalar TFactorFluid = maxCo/(CoNum + SMALL);
+    const scalar TFactorFilm = maxCo/(surfaceFilm.CourantNumber() + SMALL);
+
+    const scalar dt0 = runTime.deltaTValue();
+
+    runTime.setDeltaT
+    (
+        min
+        (
+            dt0*min(min(TFactorFluid, TFactorFilm), 1.2),
+            maxDeltaT
+        )
+    );
+}
+
+
+// ************************************************************************* //

applications/solvers/multiphase/bubbleFoam/wallFunctions.H

@@ -34,7 +34,7 @@
         {
             const scalarField& nuw = nutb.boundaryField()[patchi];
 
-            scalarField magFaceGradU(mag(U.boundaryField()[patchi].snGrad()));
+            scalarField magFaceGradU(mag(Ub.boundaryField()[patchi].snGrad()));
 
             forAll(currPatch, facei)
             {

applications/solvers/multiphase/interMixingFoam/incompressibleThreePhaseMixture/threePhaseMixture.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2004-2010 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2004-2011 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -36,6 +36,10 @@ Class
 //- Calculate and return the laminar viscosity
 void Foam::threePhaseMixture::calcNu()
 {
+    nuModel1_->correct();
+    nuModel2_->correct();
+    nuModel3_->correct();
+
     // Average kinematic viscosity calculated from dynamic viscosity
     nu_ = mu()/(alpha1_*rho1_ + alpha2_*rho2_ + alpha3_*rho3_);
 }

etc/settings.csh

@@ -52,7 +52,7 @@ case Linux:
     case x86_64:
         switch ($WM_ARCH_OPTION)
         case 32:
-            setenv WM_COMPILER_ARCH '-64'
+            setenv WM_COMPILER_ARCH 64
             setenv WM_CC 'gcc'
             setenv WM_CXX 'g++'
             setenv WM_CFLAGS '-m32 -fPIC'

etc/settings.sh

@@ -76,7 +76,7 @@ Linux)
     x86_64)
         case "$WM_ARCH_OPTION" in
         32)
-            export WM_COMPILER_ARCH='-64'
+            export WM_COMPILER_ARCH=64
             export WM_CC='gcc'
             export WM_CXX='g++'
             export WM_CFLAGS='-m32 -fPIC'

src/OSspecific/POSIX/regExp.H

@@ -62,6 +62,7 @@ class regExp
         //- Precompiled regular expression
         mutable regex_t* preg_;
 
+
     // Private Member Functions
 
         //- Disallow default bitwise copy construct
@@ -72,6 +73,8 @@ class regExp
 
 public:
 
+    // Static Member Functions
+
         //- Is character a regular expression meta-character?
         //  any character: '.' \n
         //  quantifiers: '*', '+', '?' \n
@@ -102,13 +105,14 @@ public:
         //- Construct from std::string (or string), optionally ignoring case
         regExp(const std::string&, const bool ignoreCase=false);
 
+
     //- Destructor
     ~regExp();
 
 
     // Member functions
 
-    //- Access
+        // Access
 
             //- Return true if a precompiled expression does not exist
             inline bool empty() const
@@ -129,12 +133,14 @@ public:
             }
 
 
-    //- Editing
+        // Editing
 
-        //- Compile pattern into a regular expression, optionally ignoring case
+            //- Compile pattern into a regular expression, optionally ignoring
+            //  case
             void set(const char*, const bool ignoreCase=false) const;
 
-        //- Compile pattern into a regular expression, optionally ignoring case
+            //- Compile pattern into a regular expression, optionally ignoring
+            //  case
             void set(const std::string&, const bool ignoreCase=false) const;
 
 
@@ -143,7 +149,7 @@ public:
             bool clear() const;
 
 
-    //- Searching
+        // Searching
 
             //- Find position within string.
             //  Returns the index where it begins or string::npos if not found
@@ -173,7 +179,6 @@ public:
         //- Assign and compile pattern from string
         //  Always case sensitive
         void operator=(const std::string&);
-
 };
 
 

src/OpenFOAM/primitives/strings/wordRe/wordRe.H

@@ -110,6 +110,7 @@ public:
         //- Test string for regular expression meta characters
         static inline bool isPattern(const string&);
 
+
     // Constructors
 
         //- Construct null
@@ -140,14 +141,16 @@ public:
         //  Words are treated as literals, strings with an auto-test
         wordRe(Istream&);
 
+
     // Member functions
 
-    //- Access
+        // Access
 
             //- Should be treated as a match rather than a literal string?
             inline bool isPattern() const;
 
-    //- Infrastructure
+
+        // Infrastructure
 
             //- Compile the regular expression
             inline bool compile() const;
@@ -162,7 +165,8 @@ public:
             //  Optionally strips invalid word characters
             inline void uncompile(const bool doStripInvalid = false) const;
 
-    //- Editing
+
+        // Editing
 
             //- Copy string, auto-test for regular expression or other options
             inline void set(const std::string&, const compOption = DETECT);
@@ -173,13 +177,19 @@ public:
             //- Clear string and precompiled regular expression
             inline void clear();
 
-    //- Searching
+
+        // Searching
 
             //- Smart match as regular expression or as a string
             //  Optionally force a literal match only
-        inline bool match(const std::string&, bool literalMatch=false) const;
+            inline bool match
+            (
+                const std::string&,
+                bool literalMatch = false
+            ) const;
 
-    //- Miscellaneous
+
+        // Miscellaneous
 
             //- Return a string with quoted meta-characters
             inline string quotemeta() const;

src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloud.C

@@ -100,13 +100,7 @@ template<class CloudType>
 template<class TrackData>
 void Foam::KinematicCloud<CloudType>::solve(TrackData& td)
 {
-    if (solution_.transient())
+    if (solution_.steadyState())
-    {
-        td.cloud().preEvolve();
-
-        evolveCloud(td);
-    }
-    else
     {
         td.cloud().storeState();
 
@@ -114,8 +108,22 @@ void Foam::KinematicCloud<CloudType>::solve(TrackData& td)
 
         evolveCloud(td);
 
+        if (solution_.coupled())
+        {
             td.cloud().relaxSources(td.cloud().cloudCopy());
         }
+    }
+    else
+    {
+        td.cloud().preEvolve();
+
+        evolveCloud(td);
+
+        if (solution_.coupled())
+        {
+            td.cloud().scaleSources();
+        }
+    }
 
     td.cloud().info();
 
@@ -258,6 +266,7 @@ void Foam::KinematicCloud<CloudType>::cloudReset(KinematicCloud<CloudType>& c)
     injectionModel_.reset(c.injectionModel_.ptr());
     patchInteractionModel_.reset(c.patchInteractionModel_.ptr());
     postProcessingModel_.reset(c.postProcessingModel_.ptr());
+    surfaceFilmModel_.reset(c.surfaceFilmModel_.ptr());
 
     UIntegrator_.reset(c.UIntegrator_.ptr());
 }
@@ -556,11 +565,23 @@ void Foam::KinematicCloud<CloudType>::relax
 ) const
 {
     const scalar coeff = solution_.relaxCoeff(name);
-
     field = field0 + coeff*(field - field0);
 }
 
 
+template<class CloudType>
+template<class Type>
+void Foam::KinematicCloud<CloudType>::scale
+(
+    DimensionedField<Type, volMesh>& field,
+    const word& name
+) const
+{
+    const scalar coeff = solution_.relaxCoeff(name);
+    field *= coeff;
+}
+
+
 template<class CloudType>
 void Foam::KinematicCloud<CloudType>::relaxSources
 (
@@ -568,6 +589,15 @@ void Foam::KinematicCloud<CloudType>::relaxSources
 )
 {
     this->relax(UTrans_(), cloudOldTime.UTrans(), "U");
+    this->relax(UCoeff_(), cloudOldTime.UCoeff(), "U");
+}
+
+
+template<class CloudType>
+void Foam::KinematicCloud<CloudType>::scaleSources()
+{
+    this->scale(UTrans_(), "U");
+    this->scale(UCoeff_(), "U");
 }
 
 

src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloud.H

@@ -498,9 +498,20 @@ public:
                 const word& name
             ) const;
 
+            //- Scale field
+            template<class Type>
+            void scale
+            (
+                DimensionedField<Type, volMesh>& field,
+                const word& name
+            ) const;
+
             //- Apply relaxation to (steady state) cloud sources
             void relaxSources(const KinematicCloud<CloudType>& cloudOldTime);
 
+            //- Apply scaling to (transient) cloud sources
+            void scaleSources();
+
             //- Evolve the cloud
             void evolve();
 

src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloudI.H

@@ -370,42 +370,22 @@ Foam::KinematicCloud<CloudType>::theta() const
                 false
             ),
             mesh_,
-            dimensionedScalar("zero", dimless, 0.0)
+            dimensionedScalar("zero", dimless, 0.0),
+            zeroGradientFvPatchScalarField::typeName
         )
     );
 
     volScalarField& theta = ttheta();
-    theta.boundaryField() == 0;
-
     forAllConstIter(typename KinematicCloud<CloudType>, *this, iter)
     {
         const parcelType& p = iter();
         const label cellI = p.cell();
 
-        if ((p.face() != -1))
-        {
-            const label patchI = p.patch(p.face());
-            if (patchI != -1)
-            {
-                scalarField& thetap = theta.boundaryField()[patchI];
-                const label faceI = p.patchFace(patchI, p.face());
-                thetap[faceI] += p.nParticle()*p.areaP();
-            }
-        }
-
         theta[cellI] += p.nParticle()*p.volume();
     }
 
     theta.internalField() /= mesh_.V();
-
+    theta.correctBoundaryConditions();
-    forAll(theta.boundaryField(), patchI)
-    {
-        scalarField& thetap = theta.boundaryField()[patchI];
-        if (thetap.size() > 0)
-        {
-            thetap /= mesh_.magSf().boundaryField()[patchI];
-        }
-    }
 
     return ttheta;
 }

src/lagrangian/intermediate/clouds/Templates/KinematicCloud/cloudSolution/cloudSolution.C

@@ -116,10 +116,16 @@ void Foam::cloudSolution::read()
         dict_.lookup("calcFrequency") >> calcFrequency_;
         dict_.lookup("maxCo") >> maxCo_;
         dict_.lookup("maxTrackTime") >> maxTrackTime_;
+
+        if (coupled_)
+        {
             dict_.subDict("sourceTerms").lookup("resetOnStartup")
                 >> resetSourcesOnStartup_;
         }
+    }
 
+    if (coupled_)
+    {
         const dictionary&
             schemesDict(dict_.subDict("sourceTerms").subDict("schemes"));
 
@@ -152,6 +158,7 @@ void Foam::cloudSolution::read()
             is  >> schemes_[i].second().second();
         }
     }
+}
 
 
 Foam::scalar Foam::cloudSolution::relaxCoeff(const word& fieldName) const

src/lagrangian/intermediate/clouds/Templates/ReactingCloud/ReactingCloud.C

@@ -295,10 +295,10 @@ void Foam::ReactingCloud<CloudType>::relaxSources
     const ReactingCloud<CloudType>& cloudOldTime
 )
 {
-    typedef DimensionedField<scalar, volMesh> dsfType;
-
     CloudType::relaxSources(cloudOldTime);
 
+    typedef DimensionedField<scalar, volMesh> dsfType;
+
     forAll(rhoTrans_, fieldI)
     {
         dsfType& rhoT = rhoTrans_[fieldI];
@@ -308,6 +308,21 @@ void Foam::ReactingCloud<CloudType>::relaxSources
 }
 
 
+template<class CloudType>
+void Foam::ReactingCloud<CloudType>::scaleSources()
+{
+    CloudType::scaleSources();
+
+    typedef DimensionedField<scalar, volMesh> dsfType;
+
+    forAll(rhoTrans_, fieldI)
+    {
+        dsfType& rhoT = rhoTrans_[fieldI];
+        this->scale(rhoT, "rho");
+    }
+}
+
+
 template<class CloudType>
 void Foam::ReactingCloud<CloudType>::evolve()
 {

src/lagrangian/intermediate/clouds/Templates/ReactingCloud/ReactingCloud.H

@@ -287,6 +287,9 @@ public:
             //- Apply relaxation to (steady state) cloud sources
             void relaxSources(const ReactingCloud<CloudType>& cloudOldTime);
 
+            //- Apply scaling to (transient) cloud sources
+            void scaleSources();
+
             //- Evolve the cloud
             void evolve();
 

src/lagrangian/intermediate/clouds/Templates/ThermoCloud/ThermoCloud.C

@@ -290,6 +290,17 @@ void Foam::ThermoCloud<CloudType>::relaxSources
     CloudType::relaxSources(cloudOldTime);
 
     this->relax(hsTrans_(), cloudOldTime.hsTrans(), "hs");
+    this->relax(hsCoeff_(), cloudOldTime.hsCoeff(), "hs");
+}
+
+
+template<class CloudType>
+void Foam::ThermoCloud<CloudType>::scaleSources()
+{
+    CloudType::scaleSources();
+
+    this->scale(hsTrans_(), "hs");
+    this->scale(hsCoeff_(), "hs");
 }
 
 

src/lagrangian/intermediate/clouds/Templates/ThermoCloud/ThermoCloud.H

@@ -301,6 +301,9 @@ public:
             //- Apply relaxation to (steady state) cloud sources
             void relaxSources(const ThermoCloud<CloudType>& cloudOldTime);
 
+            //- Apply scaling to (transient) cloud sources
+            void scaleSources();
+
             //- Evolve the cloud
             void evolve();
 

src/lagrangian/intermediate/parcels/Templates/ReactingMultiphaseParcel/ReactingMultiphaseParcel.C

@@ -327,7 +327,6 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calc
 
 
 
-
     // Heat transfer
     // ~~~~~~~~~~~~~
 
@@ -369,7 +368,7 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calc
             d0,
             U0,
             rho0,
-            mass0,
+            0.5*(mass0 + mass1),
             Su,
             dUTrans,
             Spu
@@ -386,15 +385,11 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calc
         {
             label gid = composition.localToGlobalCarrierId(GAS, i);
             td.cloud().rhoTrans(gid)[cellI] += np0*dMassGas[i];
-//            td.cloud().hsTrans()[cellI] +=
-//                np0*dMassGas[i]*composition.carrier().Hs(gid, T0);
         }
         forAll(YLiquid_, i)
         {
             label gid = composition.localToGlobalCarrierId(LIQ, i);
             td.cloud().rhoTrans(gid)[cellI] += np0*dMassLiquid[i];
-//            td.cloud().hsTrans()[cellI] +=
-//                np0*dMassLiquid[i]*composition.carrier().Hs(gid, T0);
         }
 /*
         // No mapping between solid components and carrier phase
@@ -402,15 +397,11 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calc
         {
             label gid = composition.localToGlobalCarrierId(SLD, i);
             td.cloud().rhoTrans(gid)[cellI] += np0*dMassSolid[i];
-//            td.cloud().hsTrans()[cellI] +=
-//                np0*dMassSolid[i]*composition.carrier().Hs(gid, T0);
         }
 */
         forAll(dMassSRCarrier, i)
         {
             td.cloud().rhoTrans(i)[cellI] += np0*dMassSRCarrier[i];
-//            td.cloud().hsTrans()[cellI] +=
-//                np0*dMassSRCarrier[i]*composition.carrier().Hs(i, T0);
         }
 
         // Update momentum transfer
@@ -458,8 +449,8 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calc
             }
 */
             td.cloud().UTrans()[cellI] += np0*mass1*U1;
-            td.cloud().hsTrans()[cellI] +=
+
-                np0*mass1*HEff(td, pc, T1, idG, idL, idS); // using total h
+            // enthalpy transfer accounted for via change in mass fractions
         }
     }
 

src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcel.C

@@ -379,13 +379,12 @@ void Foam::ReactingParcel<ParcelType>::calc
             d0,
             U0,
             rho0,
-            mass0,
+            0.5*(mass0 + mass1),
             Su,
             dUTrans,
             Spu
         );
 
-    dUTrans += 0.5*(mass0 - mass1)*(U0 + U1);
 
     // Accumulate carrier phase source terms
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
@@ -396,8 +395,6 @@ void Foam::ReactingParcel<ParcelType>::calc
         {
             label gid = composition.localToGlobalCarrierId(0, i);
             td.cloud().rhoTrans(gid)[cellI] += np0*dMassPC[i];
-//            td.cloud().hsTrans()[cellI] +=
-//                np0*dMassPC[i]*composition.carrier().Hs(gid, T0);
         }
 
         // Update momentum transfer
@@ -429,8 +426,8 @@ void Foam::ReactingParcel<ParcelType>::calc
                 td.cloud().rhoTrans(gid)[cellI] += np0*mass1*Y_[i];
             }
             td.cloud().UTrans()[cellI] += np0*mass1*U1;
-            td.cloud().hsTrans()[cellI] +=
+
-                np0*mass1*composition.H(0, Y_, pc_, T1);
+            // enthalpy transfer accounted for via change in mass fractions
         }
     }
 

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

@@ -148,7 +148,8 @@ void Foam::SurfaceFilmModel<CloudType>::inject(TrackData& td)
     const labelList& filmPatches = filmModel.intCoupledPatchIDs();
     const labelList& primaryPatches = filmModel.primaryPatchIDs();
 
-    const polyBoundaryMesh& pbm = this->owner().mesh().boundaryMesh();
+    const fvMesh& mesh = this->owner().mesh();
+    const polyBoundaryMesh& pbm = mesh.boundaryMesh();
 
     forAll(filmPatches, i)
     {
@@ -163,6 +164,10 @@ void Foam::SurfaceFilmModel<CloudType>::inject(TrackData& td)
 
         cacheFilmFields(filmPatchI, primaryPatchI, distMap, filmModel);
 
+        const vectorField& Cf = mesh.C().boundaryField()[primaryPatchI];
+        const vectorField& Sf = mesh.Sf().boundaryField()[primaryPatchI];
+        const scalarField& magSf = mesh.magSf().boundaryField()[primaryPatchI];
+
         forAll(injectorCellsPatch, j)
         {
             if (diameterParcelPatch_[j] > 0)
@@ -179,7 +184,15 @@ void Foam::SurfaceFilmModel<CloudType>::inject(TrackData& td)
                 const label tetFaceI = this->owner().mesh().cells()[cellI][0];
                 const label tetPtI = 1;
 
-                const point& pos = this->owner().mesh().C()[cellI];
+//                const point& pos = this->owner().mesh().C()[cellI];
+
+                const scalar offset =
+                    max
+                    (
+                        diameterParcelPatch_[j],
+                        deltaFilmPatch_[primaryPatchI][j]
+                    );
+                const point pos = Cf[j] - 1.1*offset*Sf[j]/magSf[j];
 
                 // Create a new parcel
                 parcelType* pPtr =
@@ -217,11 +230,11 @@ void Foam::SurfaceFilmModel<CloudType>::cacheFilmFields
     const regionModels::surfaceFilmModels::surfaceFilmModel& filmModel
 )
 {
-    massParcelPatch_ = filmModel.massForPrimary().boundaryField()[filmPatchI];
+    massParcelPatch_ = filmModel.cloudMassTrans().boundaryField()[filmPatchI];
     distMap.distribute(massParcelPatch_);
 
     diameterParcelPatch_ =
-        filmModel.diametersForPrimary().boundaryField()[filmPatchI];
+        filmModel.cloudDiameterTrans().boundaryField()[filmPatchI];
     distMap.distribute(diameterParcelPatch_);
 
     UFilmPatch_ = filmModel.Us().boundaryField()[filmPatchI];

src/lagrangian/intermediate/submodels/Thermodynamic/SurfaceFilmModel/ThermoSurfaceFilm/ThermoSurfaceFilm.C

@@ -391,8 +391,8 @@ void Foam::ThermoSurfaceFilm<CloudType>::splashInteraction
     const scalar dBarSplash = 1/cbrt(6.0)*cbrt(mRatio/Ns)*d + ROOTVSMALL;
 
     // cumulative diameter splash distribution
-    const scalar dMax = cbrt(mRatio)*d;
+    const scalar dMax = 0.9*cbrt(mRatio)*d;
-    const scalar dMin = 0.001*dMax;
+    const scalar dMin = 0.1*dMax;
     const scalar K = exp(-dMin/dBarSplash) - exp(-dMax/dBarSplash);
 
     // surface energy of secondary parcels [J]
@@ -437,7 +437,7 @@ void Foam::ThermoSurfaceFilm<CloudType>::splashInteraction
 
     // magnitude of the normal velocity of the first splashed parcel
     const scalar magUns0 =
-        sqrt(2.0*parcelsPerSplash_*EKs/mSplash/(1 + coeff1/sqr(coeff2)));
+        sqrt(2.0*parcelsPerSplash_*EKs/mSplash/(1.0 + coeff1/sqr(coeff2)));
 
     // Set splashed parcel properties
     forAll(dNew, i)
@@ -467,7 +467,7 @@ void Foam::ThermoSurfaceFilm<CloudType>::splashInteraction
 
         // Apply correction to velocity for 2-D cases
         meshTools::constrainDirection(mesh, mesh.solutionD(), pPtr->U());
-
+Info<< "NEW PARTICLE: " << *pPtr << endl;
         // Add the new parcel
         this->owner().addParticle(pPtr);
 

src/regionModels/regionModel/regionModel1D/regionModel1D.C

@@ -114,8 +114,6 @@ void Foam::regionModels::regionModel1D::initialise()
             boundaryFaceCells_[localPyrolysisFaceI].transfer(cellIDs);
 
             localPyrolysisFaceI++;
-
-            nLayers_ = nCells;
         }
     }
 
@@ -268,7 +266,6 @@ Foam::regionModels::regionModel1D::regionModel1D(const fvMesh& mesh)
     boundaryFaceFaces_(),
     boundaryFaceCells_(),
     boundaryFaceOppositeFace_(),
-    nLayers_(0),
     nMagSfPtr_(NULL),
     moveMesh_(false)
 {}
@@ -286,7 +283,6 @@ Foam::regionModels::regionModel1D::regionModel1D
     boundaryFaceFaces_(regionMesh().nCells()),
     boundaryFaceCells_(regionMesh().nCells()),
     boundaryFaceOppositeFace_(regionMesh().nCells()),
-    nLayers_(0),
     nMagSfPtr_(NULL),
     moveMesh_(true)
 {

src/regionModels/regionModel/regionModel1D/regionModel1D.H

@@ -88,9 +88,6 @@ protected:
             //- Global boundary face IDs oppossite coupled patch
             labelList boundaryFaceOppositeFace_;
 
-            //- Number of layers in the region
-            label nLayers_;
-
 
         // Geometry
 
@@ -155,9 +152,6 @@ public:
                 //- Return the global boundary face IDs oppossite coupled patch
                 inline const labelList& boundaryFaceOppositeFace() const;
 
-                //- Return the number of layers in the region
-                inline label nLayers() const;
-
 
             // Geometry
 

src/regionModels/regionModel/regionModel1D/regionModel1DI.H

@@ -49,12 +49,6 @@ Foam::regionModels::regionModel1D::boundaryFaceOppositeFace() const
 }
 
 
-inline Foam::label Foam::regionModels::regionModel1D::nLayers() const
-{
-    return nLayers_;
-}
-
-
 inline const Foam::surfaceScalarField&
 Foam::regionModels::regionModel1D::nMagSf() const
 {

src/regionModels/regionModel/singleLayerRegion/singleLayerRegion.C

@@ -111,14 +111,12 @@ void Foam::regionModels::singleLayerRegion::initialise()
 
     if (nBoundaryFaces != regionMesh().nCells())
     {
-        /*
         FatalErrorIn("singleLayerRegion::initialise()")
             << "Number of primary region coupled boundary faces not equal to "
             << "the number of cells in the local region" << nl << nl
             << "Number of cells = " << regionMesh().nCells() << nl
             << "Boundary faces  = " << nBoundaryFaces << nl
             << abort(FatalError);
-        */
     }
 
     scalarField passiveMagSf(magSf.size(), 0.0);
@@ -178,12 +176,11 @@ Foam::regionModels::singleLayerRegion::singleLayerRegion
     bool readFields
 )
 :
-    regionModel(mesh, regionType, modelName, readFields),
+    regionModel(mesh, regionType, modelName, false),
     nHatPtr_(NULL),
     magSfPtr_(NULL),
     passivePatchIDs_()
 {
-    Info << "singleLayerRegion" << endl;
     if (active_)
     {
         constructMeshObjects();

src/regionModels/surfaceFilmModels/Make/files

@@ -12,9 +12,10 @@ submodels/subModelBase.C
 KINEMATICMODELS=submodels/kinematic
 $(KINEMATICMODELS)/injectionModel/injectionModel/injectionModel.C
 $(KINEMATICMODELS)/injectionModel/injectionModel/injectionModelNew.C
-$(KINEMATICMODELS)/injectionModel/noInjection/noInjection.C
+$(KINEMATICMODELS)/injectionModel/injectionModelList/injectionModelList.C
-$(KINEMATICMODELS)/injectionModel/cloudInjection/cloudInjection.C
+$(KINEMATICMODELS)/injectionModel/drippingInjection/drippingInjection.C
 $(KINEMATICMODELS)/injectionModel/removeInjection/removeInjection.C
+$(KINEMATICMODELS)/injectionModel/curvatureSeparation/curvatureSeparation.C
 
 THERMOMODELS=submodels/thermo
 $(THERMOMODELS)/phaseChangeModel/phaseChangeModel/phaseChangeModel.C

src/regionModels/surfaceFilmModels/derivedFvPatchFields/filmHeightInletVelocity/filmHeightInletVelocityFvPatchVectorField.C

@@ -123,10 +123,10 @@ void Foam::filmHeightInletVelocityFvPatchVectorField::updateCoeffs()
     const fvPatchField<scalar>& deltafp =
         patch().lookupPatchField<volScalarField, scalar>(deltafName_);
 
-    const vectorField& n = patch().nf();
+    vectorField n(patch().nf());
     const scalarField& magSf = patch().magSf();
 
-    operator==(deltafp*n*phip/(rhop*magSf*sqr(deltafp) + ROOTVSMALL));
+    operator==(n*phip/(rhop*magSf*deltafp + ROOTVSMALL));
 
     fixedValueFvPatchVectorField::updateCoeffs();
 }

src/regionModels/surfaceFilmModels/derivedFvPatchFields/wallFunctions/alphatFilmWallFunction/alphatFilmWallFunctionFvPatchScalarField.C

@@ -156,9 +156,9 @@ void alphatFilmWallFunctionFvPatchScalarField::updateCoeffs()
 
     const mapDistribute& distMap = filmModel.mappedPatches()[filmPatchI].map();
 
-    scalarField mDotFilm =
+    tmp<volScalarField> mDotFilm(filmModel.primaryMassTrans());
-        filmModel.massPhaseChangeForPrimary().boundaryField()[filmPatchI];
+    scalarField mDotFilmp = mDotFilm().boundaryField()[filmPatchI];
-    distMap.distribute(mDotFilm);
+    distMap.distribute(mDotFilmp);
 
     // Retrieve RAS turbulence model
     const RASModel& rasModel = db().lookupObject<RASModel>("RASProperties");
@@ -185,7 +185,7 @@ void alphatFilmWallFunctionFvPatchScalarField::updateCoeffs()
         scalar Pr = muw[faceI]/alphaw[faceI];
 
         scalar factor = 0.0;
-        scalar mStar = mDotFilm[faceI]/(y[faceI]*uTau);
+        scalar mStar = mDotFilmp[faceI]/(y[faceI]*uTau);
         if (yPlus > yPlusCrit_)
         {
             scalar expTerm = exp(min(50.0, yPlusCrit_*mStar*Pr));
@@ -209,6 +209,7 @@ void alphatFilmWallFunctionFvPatchScalarField::updateCoeffs()
 }
 
 
+
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 void alphatFilmWallFunctionFvPatchScalarField::write(Ostream& os) const

src/regionModels/surfaceFilmModels/derivedFvPatchFields/wallFunctions/mutkFilmWallFunction/mutkFilmWallFunctionFvPatchScalarField.C

@@ -72,9 +72,9 @@ tmp<scalarField> mutkFilmWallFunctionFvPatchScalarField::calcUTau
 
     const mapDistribute& distMap = filmModel.mappedPatches()[filmPatchI].map();
 
-    scalarField mDotFilm =
+    tmp<volScalarField> mDotFilm(filmModel.primaryMassTrans());
-        filmModel.massPhaseChangeForPrimary().boundaryField()[filmPatchI];
+    scalarField mDotFilmp = mDotFilm().boundaryField()[filmPatchI];
-    distMap.distribute(mDotFilm);
+    distMap.distribute(mDotFilmp);
 
 
     // Retrieve RAS turbulence model
@@ -95,7 +95,7 @@ tmp<scalarField> mutkFilmWallFunctionFvPatchScalarField::calcUTau
 
         scalar yPlus = y[faceI]*ut/(muw[faceI]/rhow[faceI]);
 
-        scalar mStar = mDotFilm[faceI]/(y[faceI]*ut);
+        scalar mStar = mDotFilmp[faceI]/(y[faceI]*ut);
 
         scalar factor = 0.0;
         if (yPlus > yPlusCrit_)

src/regionModels/surfaceFilmModels/kinematicSingleLayer/kinematicSingleLayer.C

@@ -35,9 +35,6 @@ License
 #include "directMappedWallPolyPatch.H"
 #include "mapDistribute.H"
 
-// Sub-models
-#include "injectionModel.H"
-
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 namespace Foam
@@ -66,7 +63,6 @@ bool kinematicSingleLayer::read()
         solution.lookup("nNonOrthCorr") >> nNonOrthCorr_;
 
         coeffs_.lookup("Cf") >> Cf_;
-        coeffs_.lookup("deltaStable") >> deltaStable_;
 
         return true;
     }
@@ -99,6 +95,11 @@ void kinematicSingleLayer::correctThermoFields()
 
 void kinematicSingleLayer::resetPrimaryRegionSourceTerms()
 {
+    if (debug)
+    {
+        Info<< "kinematicSingleLayer::resetPrimaryRegionSourceTerms()" << endl;
+    }
+
     rhoSpPrimary_ == dimensionedScalar("zero", rhoSp_.dimensions(), 0.0);
     USpPrimary_ == dimensionedVector("zero", USp_.dimensions(), vector::zero);
     pSpPrimary_ == dimensionedScalar("zero", pSp_.dimensions(), 0.0);
@@ -107,6 +108,11 @@ void kinematicSingleLayer::resetPrimaryRegionSourceTerms()
 
 void kinematicSingleLayer::transferPrimaryRegionThermoFields()
 {
+    if (debug)
+    {
+        Info<< "kinematicSingleLayer::"
+            << "transferPrimaryRegionThermoFields()" << endl;
+    }
     // Update fields from primary region via direct mapped
     // (coupled) boundary conditions
     UPrimary_.correctBoundaryConditions();
@@ -118,6 +124,12 @@ void kinematicSingleLayer::transferPrimaryRegionThermoFields()
 
 void kinematicSingleLayer::transferPrimaryRegionSourceFields()
 {
+    if (debug)
+    {
+        Info<< "kinematicSingleLayer::"
+            << "transferPrimaryRegionSourceFields()" << endl;
+    }
+
     // Retrieve the source fields from the primary region via direct mapped
     // (coupled) boundary conditions
     // - fields require transfer of values for both patch AND to push the
@@ -132,10 +144,6 @@ void kinematicSingleLayer::transferPrimaryRegionSourceFields()
     rhoSp_.field() /= magSf()*deltaT;
     USp_.field() /= magSf()*deltaT;
     pSp_.field() /= magSf()*deltaT;
-
-    // reset transfer to primary fields
-    massForPrimary_ == dimensionedScalar("zero", dimMass, 0.0);
-    diametersForPrimary_ == dimensionedScalar("zero", dimLength, -1.0);
 }
 
 
@@ -154,7 +162,7 @@ tmp<volScalarField> kinematicSingleLayer::pu()
                 IOobject::NO_WRITE
             ),
             pPrimary_                  // pressure (mapped from primary region)
-          + pSp_                           // accumulated particle impingement
+          - pSp_                           // accumulated particle impingement
           - fvc::laplacian(sigma_, delta_) // surface tension
         )
     );
@@ -181,51 +189,25 @@ tmp<volScalarField> kinematicSingleLayer::pp()
 }
 
 
-void kinematicSingleLayer::correctDetachedFilm()
-{
-    tmp<volScalarField> tgNorm(this->gNorm());
-    const scalarField& gNorm = tgNorm();
-    const scalarField& magSf = this->magSf();
-
-    forAll(gNorm, i)
-    {
-        if (gNorm[i] > SMALL)
-        {
-            const scalar ddelta = max(0.0, delta_[i] - deltaStable_.value());
-            massForPrimary_[i] =
-                max
-                (
-                    0.0,
-                    ddelta*rho_[i]*magSf[i] - massPhaseChangeForPrimary_[i]
-                );
-        }
-    }
-}
-
-
 void kinematicSingleLayer::updateSubmodels()
 {
-    correctDetachedFilm();
+    if (debug)
+    {
+        Info<< "kinematicSingleLayer::updateSubmodels()" << endl;
+    }
 
-    // Update injection model - mass returned is actual mass injected
+    // Update injection model - mass returned is mass available for injection
-    injection_->correct(massForPrimary_, diametersForPrimary_);
+    injection_.correct(availableMass_, cloudMassTrans_, cloudDiameterTrans_);
-
-    // Update cumulative detached mass counter
-    injectedMassTotal_ += sum(massForPrimary_.field());
-
-    // Push values to boundaries ready for transfer to the primary region
-    massForPrimary_.correctBoundaryConditions();
-    diametersForPrimary_.correctBoundaryConditions();
 
     // Update source fields
     const dimensionedScalar deltaT = time().deltaT();
-    rhoSp_ -= (massForPrimary_ + massPhaseChangeForPrimary_)/magSf()/deltaT;
+    rhoSp_ += cloudMassTrans_/magSf()/deltaT;
 }
 
 
 void kinematicSingleLayer::continuityCheck()
 {
-    const volScalarField deltaRho0 = deltaRho_;
+    const volScalarField deltaRho0(deltaRho_);
 
     solveContinuity();
 
@@ -268,7 +250,7 @@ void kinematicSingleLayer::solveContinuity()
         fvm::ddt(deltaRho_)
       + fvc::div(phi_)
      ==
-        rhoSp_
+      - rhoSp_
     );
 }
 
@@ -305,8 +287,8 @@ tmp<fvVectorMatrix> kinematicSingleLayer::tau(volVectorField& U) const
 
     return
     (
-       - fvm::Sp(Cs, U) + Cs*Us_
+       - fvm::Sp(Cs, U) + Cs*Us_ // surface contribution
-       - fvm::Sp(Cw, U) + Cw*Uw_
+       - fvm::Sp(Cw, U) + Cw*Uw_ // wall contribution
     );
 }
 
@@ -330,15 +312,10 @@ tmp<Foam::fvVectorMatrix> kinematicSingleLayer::solveMomentum
         fvm::ddt(deltaRho_, U_)
       + fvm::div(phi_, U_)
      ==
-        USp_
+      - USp_
       + tau(U_)
       + fvc::grad(sigma_)
-      - fvm::Sp
+      - fvm::SuSp(rhoSp_, U_)
-        (
-            (massForPrimary_ + massPhaseChangeForPrimary_)
-           /magSf()/time().deltaT(),
-            U_
-        )
     );
 
     fvVectorMatrix& UEqn = tUEqn();
@@ -415,6 +392,7 @@ void kinematicSingleLayer::solveThickness
 
     surfaceScalarField ddrhorUAppf
     (
+        "deltaCoeff",
         fvc::interpolate(delta_)*deltarUAf*rhof*fvc::interpolate(pp)
     );
 //    constrainFilmField(ddrhorUAppf, 0.0);
@@ -428,7 +406,7 @@ void kinematicSingleLayer::solveThickness
           + fvm::div(phid, delta_)
           - fvm::laplacian(ddrhorUAppf, delta_)
          ==
-            rhoSp_
+          - rhoSp_
         );
 
         deltaEqn.solve();
@@ -483,7 +461,6 @@ kinematicSingleLayer::kinematicSingleLayer
     cumulativeContErr_(0.0),
 
     Cf_(readScalar(coeffs().lookup("Cf"))),
-    deltaStable_(coeffs().lookup("deltaStable")),
 
     rho_
     (
@@ -607,11 +584,11 @@ kinematicSingleLayer::kinematicSingleLayer
         dimLength*dimMass/dimTime
     ),
 
-    massForPrimary_
+    primaryMassTrans_
     (
         IOobject
         (
-            "massForPrimary",
+            "primaryMassTrans",
             time().timeName(),
             regionMesh(),
             IOobject::NO_READ,
@@ -621,11 +598,25 @@ kinematicSingleLayer::kinematicSingleLayer
         dimensionedScalar("zero", dimMass, 0.0),
         zeroGradientFvPatchScalarField::typeName
     ),
-    diametersForPrimary_
+    cloudMassTrans_
     (
         IOobject
         (
-            "diametersForPrimary",
+            "cloudMassTrans",
+            time().timeName(),
+            regionMesh(),
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        regionMesh(),
+        dimensionedScalar("zero", dimMass, 0.0),
+        zeroGradientFvPatchScalarField::typeName
+    ),
+    cloudDiameterTrans_
+    (
+        IOobject
+        (
+            "cloudDiameterTrans",
             time().timeName(),
             regionMesh(),
             IOobject::NO_READ,
@@ -635,20 +626,6 @@ kinematicSingleLayer::kinematicSingleLayer
         dimensionedScalar("zero", dimLength, -1.0),
         zeroGradientFvPatchScalarField::typeName
     ),
-    massPhaseChangeForPrimary_
-    (
-        IOobject
-        (
-            "massPhaseChangeForPrimary",
-            time().timeName(),
-            regionMesh(),
-            IOobject::NO_READ,
-            IOobject::NO_WRITE
-        ),
-        regionMesh(),
-        dimensionedScalar("zero", dimMass, 0),
-        zeroGradientFvPatchScalarField::typeName
-    ),
 
     USp_
     (
@@ -793,10 +770,11 @@ kinematicSingleLayer::kinematicSingleLayer
         this->mappedFieldAndInternalPatchTypes<scalar>()
     ),
 
-    injection_(injectionModel::New(*this, coeffs_)),
+    availableMass_(regionMesh().nCells(), 0.0),
 
-    addedMassTotal_(0.0),
+    injection_(*this, coeffs_),
-    injectedMassTotal_(0.0)
+
+    addedMassTotal_(0.0)
 {
     if (readFields)
     {
@@ -836,9 +814,9 @@ void kinematicSingleLayer::addSources
             << "    pressure = " << pressureSource << endl;
     }
 
-    rhoSpPrimary_.boundaryField()[patchI][faceI] += massSource;
+    rhoSpPrimary_.boundaryField()[patchI][faceI] -= massSource;
-    USpPrimary_.boundaryField()[patchI][faceI] += momentumSource;
+    USpPrimary_.boundaryField()[patchI][faceI] -= momentumSource;
-    pSpPrimary_.boundaryField()[patchI][faceI] += pressureSource;
+    pSpPrimary_.boundaryField()[patchI][faceI] -= pressureSource;
 
     addedMassTotal_ += massSource;
 }
@@ -846,22 +824,38 @@ void kinematicSingleLayer::addSources
 
 void kinematicSingleLayer::preEvolveRegion()
 {
+    if (debug)
+    {
+        Info<< "kinematicSingleLayer::preEvolveRegion()" << endl;
+    }
+
     transferPrimaryRegionThermoFields();
 
     correctThermoFields();
 
     transferPrimaryRegionSourceFields();
+
+    // Reset transfer fields
+//    availableMass_ = mass();
+    availableMass_ = netMass();
+    cloudMassTrans_ == dimensionedScalar("zero", dimMass, 0.0);
+    cloudDiameterTrans_ == dimensionedScalar("zero", dimLength, -1.0);
 }
 
 
 void kinematicSingleLayer::evolveRegion()
 {
+    if (debug)
+    {
+        Info<< "kinematicSingleLayer::evolveRegion()" << endl;
+    }
+
     updateSubmodels();
 
     // Solve continuity for deltaRho_
     solveContinuity();
 
-    // Implicit pressure source coefficient
+    // Implicit pressure source coefficient - constant
     tmp<volScalarField> tpp(this->pp());
 
     for (int oCorr=0; oCorr<nOuterCorr_; oCorr++)
@@ -941,6 +935,12 @@ const volVectorField& kinematicSingleLayer::Uw() const
 }
 
 
+const surfaceScalarField& kinematicSingleLayer::phi() const
+{
+    return phi_;
+}
+
+
 const volScalarField& kinematicSingleLayer::rho() const
 {
     return rho_;
@@ -1002,21 +1002,37 @@ const volScalarField& kinematicSingleLayer::kappa() const
 }
 
 
-const volScalarField& kinematicSingleLayer::massForPrimary() const
+tmp<volScalarField> kinematicSingleLayer::primaryMassTrans() const
 {
-    return massForPrimary_;
+    return tmp<volScalarField>
+    (
+        new volScalarField
+        (
+            IOobject
+            (
+                "kinematicSingleLayer::primaryMassTrans",
+                time().timeName(),
+                primaryMesh(),
+                IOobject::NO_READ,
+                IOobject::NO_WRITE,
+                false
+            ),
+            primaryMesh(),
+            dimensionedScalar("zero", dimMass/dimVolume/dimTime, 0.0)
+        )
+    );
 }
 
 
-const volScalarField& kinematicSingleLayer::diametersForPrimary() const
+const volScalarField& kinematicSingleLayer::cloudMassTrans() const
 {
-    return diametersForPrimary_;
+    return cloudMassTrans_;
 }
 
 
-const volScalarField& kinematicSingleLayer::massPhaseChangeForPrimary() const
+const volScalarField& kinematicSingleLayer::cloudDiameterTrans() const
 {
-    return massPhaseChangeForPrimary_;
+    return cloudDiameterTrans_;
 }
 
 
@@ -1028,18 +1044,18 @@ void kinematicSingleLayer::info() const
         << returnReduce<scalar>(addedMassTotal_, sumOp<scalar>()) << nl
         << indent << "current mass       = "
         << gSum((deltaRho_*magSf())()) << nl
-        << indent << "injected mass      = "
-        << returnReduce<scalar>(injectedMassTotal_, sumOp<scalar>()) << nl
         << indent << "min/max(mag(U))    = " << min(mag(U_)).value() << ", "
         << max(mag(U_)).value() << nl
         << indent << "min/max(delta)     = " << min(delta_).value() << ", "
         << max(delta_).value() << nl;
+
+    injection_.info(Info);
 }
 
 
 tmp<DimensionedField<scalar, volMesh> > kinematicSingleLayer::Srho() const
 {
-    tmp<DimensionedField<scalar, volMesh> > tSrho
+    return tmp<DimensionedField<scalar, volMesh> >
     (
         new DimensionedField<scalar, volMesh>
         (
@@ -1056,37 +1072,12 @@ tmp<DimensionedField<scalar, volMesh> >  kinematicSingleLayer::Srho() const
             dimensionedScalar("zero", dimMass/dimVolume/dimTime, 0.0)
         )
     );
-
-    scalarField& Srho = tSrho();
-    const scalarField& V = primaryMesh().V();
-    const scalar dt = time_.deltaTValue();
-
-    forAll(intCoupledPatchIDs(), i)
-    {
-        const label filmPatchI = intCoupledPatchIDs()[i];
-        const mapDistribute& distMap = mappedPatches_[filmPatchI].map();
-
-        scalarField patchMass =
-            massPhaseChangeForPrimary_.boundaryField()[filmPatchI];
-        distMap.distribute(patchMass);
-
-        const label primaryPatchI = primaryPatchIDs()[i];
-        const unallocLabelList& cells =
-            primaryMesh().boundaryMesh()[primaryPatchI].faceCells();
-
-        forAll(patchMass, j)
-        {
-            Srho[cells[j]] = patchMass[j]/(V[cells[j]]*dt);
-        }
-    }
-
-    return tSrho;
 }
 
 
 tmp<DimensionedField<scalar, volMesh> > kinematicSingleLayer::Srho
 (
-    const label
+    const label i
 ) const
 {
     return tmp<DimensionedField<scalar, volMesh> >
@@ -1095,7 +1086,7 @@ tmp<DimensionedField<scalar, volMesh> > kinematicSingleLayer::Srho
         (
             IOobject
             (
-                "kinematicSingleLayer::Srho(i)",
+                "kinematicSingleLayer::Srho(" + Foam::name(i) + ")",
                 time().timeName(),
                 primaryMesh(),
                 IOobject::NO_READ,

src/regionModels/surfaceFilmModels/kinematicSingleLayer/kinematicSingleLayer.H

@@ -42,6 +42,8 @@ SourceFiles
 #include "surfaceFields.H"
 #include "fvMatrices.H"
 
+#include "injectionModelList.H"
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 namespace Foam
@@ -51,9 +53,6 @@ namespace regionModels
 namespace surfaceFilmModels
 {
 
-// Forward declaration of classes
-class injectionModel;
-
 /*---------------------------------------------------------------------------*\
                    Class kinematicSingleLayer Declaration
 \*---------------------------------------------------------------------------*/
@@ -100,9 +99,6 @@ protected:
             //- Skin frition coefficient for film/primary region interface
             scalar Cf_;
 
-            //- Stable film thickness
-            dimensionedScalar deltaStable_;
-
 
         // Thermo properties
 
@@ -139,16 +135,16 @@ protected:
             surfaceScalarField phi_;
 
 
-            // Transfer fields - to the primary region
+            // Transfer fields
 
-                //- Film mass available for transfer
+                //- Film mass available for transfer to the primary region
-                volScalarField massForPrimary_;
+                volScalarField primaryMassTrans_;
 
-                //- Parcel diameters originating from film
+                //- Film mass available for transfer to cloud
-                volScalarField diametersForPrimary_;
+                volScalarField cloudMassTrans_;
 
-                //- Film mass evolved via phase change
+                //- Parcel diameters originating from film to cloud
-                volScalarField massPhaseChangeForPrimary_;
+                volScalarField cloudDiameterTrans_;
 
 
         // Source term fields
@@ -198,8 +194,11 @@ protected:
 
         // Sub-models
 
-            //- Injection
+            //- Available mass for transfer via sub-models
-            autoPtr<injectionModel> injection_;
+            scalarField availableMass_;
+
+            //- Cloud injection
+            injectionModelList injection_;
 
 
        // Checks
@@ -208,12 +207,6 @@ protected:
            scalar addedMassTotal_;
 
 
-       // Detached surface properties
-
-           //- Cumulative mass detached [kg]
-           scalar injectedMassTotal_;
-
-
     // Protected member functions
 
         //- Read control parameters from dictionary
@@ -231,9 +224,6 @@ protected:
         //- Transfer source fields from the primary region to the film region
         virtual void transferPrimaryRegionSourceFields();
 
-        //- Correct the source terms for film that detaches from film region
-        virtual void correctDetachedFilm();
-
         // Explicit pressure source contribution
         virtual tmp<volScalarField> pu();
 
@@ -354,6 +344,9 @@ public:
             //- Return the film wall velocity [m/s]
             virtual const volVectorField& Uw() const;
 
+            //- Return the film flux [kg.m/s]
+            virtual const surfaceScalarField& phi() const;
+
             //- Return the film density [kg/m3]
             virtual const volScalarField& rho() const;
 
@@ -375,14 +368,14 @@ public:
 
             // Transfer fields - to the primary region
 
-               //- Return the film mass available for transfer
+                //- Return mass transfer source - Eulerian phase only
-               virtual const volScalarField& massForPrimary() const;
+                virtual tmp<volScalarField> primaryMassTrans() const;
 
-               //- Return the parcel diameters originating from film
+                //- Return the film mass available for transfer to cloud
-               virtual const volScalarField& diametersForPrimary() const;
+                virtual const volScalarField& cloudMassTrans() const;
 
-               //- Return the film mass evolved via phase change
+                //- Return the parcel diameters originating from film to cloud
-               virtual const volScalarField& massPhaseChangeForPrimary() const;
+                virtual const volScalarField& cloudDiameterTrans() const;
 
 
         // External helper functions
@@ -452,13 +445,16 @@ public:
         // Sub-models
 
             //- Injection
-            inline injectionModel& injection();
+            inline injectionModelList& injection();
 
 
         // Helper functions
 
-            //- Return the gravity tangential component contributions
+            //- Return the current film mass
-            inline tmp<volVectorField> gTan() const;
+            inline tmp<volScalarField> mass() const;
+
+            //- Return the net film mass available over the next integration
+            inline tmp<volScalarField> netMass() const;
 
             //- Return the gravity normal-to-patch component contribution
             inline tmp<volScalarField> gNorm() const;
@@ -467,6 +463,9 @@ public:
             //  Clipped so that only non-zero if g & nHat_ < 0
             inline tmp<volScalarField> gNormClipped() const;
 
+            //- Return the gravity tangential component contributions
+            inline tmp<volVectorField> gTan() const;
+
 
         // Evolution
 

src/regionModels/surfaceFilmModels/kinematicSingleLayer/kinematicSingleLayerI.H

@@ -24,7 +24,8 @@ License
 
 \*---------------------------------------------------------------------------*/
 
-#include "kinematicSingleLayer.H"
+#include "surfaceInterpolate.H"
+#include "fvcSurfaceIntegrate.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -163,9 +164,24 @@ inline const volScalarField& kinematicSingleLayer::muPrimary() const
 }
 
 
-inline injectionModel& kinematicSingleLayer::injection()
+inline injectionModelList& kinematicSingleLayer::injection()
 {
-    return injection_();
+    return injection_;
+}
+
+
+inline tmp<volScalarField> kinematicSingleLayer::mass() const
+{
+    return rho_*delta_*magSf();
+}
+
+
+inline tmp<volScalarField> kinematicSingleLayer::netMass() const
+{
+    dimensionedScalar d0("SMALL", dimLength, ROOTVSMALL);
+    return
+        fvc::surfaceSum(phi_/(fvc::interpolate(delta_) + d0))*time().deltaT()
+      + rho_*delta_*magSf();
 }
 
 

src/regionModels/surfaceFilmModels/noFilm/noFilm.C

@@ -103,6 +103,15 @@ const volScalarField& noFilm::delta() const
 }
 
 
+const volScalarField& noFilm::sigma() const
+{
+    FatalErrorIn("const volScalarField& noFilm::sigma() const")
+        << "sigma field not available for " << type() << abort(FatalError);
+
+    return volScalarField::null();
+}
+
+
 const volVectorField& noFilm::U() const
 {
     FatalErrorIn("const volVectorField& noFilm::U() const")
@@ -184,32 +193,42 @@ const volScalarField& noFilm::kappa() const
 }
 
 
-const volScalarField& noFilm::massForPrimary() const
+tmp<volScalarField> noFilm::primaryMassTrans() const
 {
-    FatalErrorIn("const volScalarField& noFilm::massForPrimary() const")
+    return tmp<volScalarField>
-        << "massForPrimary field not available for " << type()
-        << abort(FatalError);
-
-    return volScalarField::null();
-}
-
-
-const volScalarField& noFilm::diametersForPrimary() const
-{
-    FatalErrorIn("const volScalarField& noFilm::diametersForPrimary() const")
-        << "diametersForPrimary field not available for " << type()
-        << abort(FatalError);
-
-    return volScalarField::null();
-}
-
-
-const volScalarField& noFilm::massPhaseChangeForPrimary() const
-{
-    FatalErrorIn
     (
-        "const volScalarField& noFilm::massPhaseChangeForPrimary() const"
+        new volScalarField
-    )   << "massPhaseChange field not available for " << type()
+        (
+            IOobject
+            (
+                "noFilm::primaryMassTrans",
+                time().timeName(),
+                primaryMesh(),
+                IOobject::NO_READ,
+                IOobject::NO_WRITE,
+                false
+            ),
+            primaryMesh(),
+            dimensionedScalar("zero", dimMass/dimVolume/dimTime, 0.0)
+        )
+    );
+}
+
+
+const volScalarField& noFilm::cloudMassTrans() const
+{
+    FatalErrorIn("const volScalarField& noFilm::cloudMassTrans() const")
+        << "cloudMassTrans field not available for " << type()
+        << abort(FatalError);
+
+    return volScalarField::null();
+}
+
+
+const volScalarField& noFilm::cloudDiameterTrans() const
+{
+    FatalErrorIn("const volScalarField& noFilm::cloudDiameterTrans() const")
+        << "cloudDiameterTrans field not available for " << type()
         << abort(FatalError);
 
     return volScalarField::null();
@@ -238,7 +257,7 @@ tmp<DimensionedField<scalar, volMesh> > noFilm::Srho() const
 }
 
 
-tmp<DimensionedField<scalar, volMesh> > noFilm::Srho(const label) const
+tmp<DimensionedField<scalar, volMesh> > noFilm::Srho(const label i) const
 {
     return tmp<DimensionedField<scalar, volMesh> >
     (
@@ -246,7 +265,7 @@ tmp<DimensionedField<scalar, volMesh> > noFilm::Srho(const label) const
         (
             IOobject
             (
-                "noFilm::Srho(i)",
+                "noFilm::Srho(" + Foam::name(i) + ")",
                 time().timeName(),
                 primaryMesh(),
                 IOobject::NO_READ,

src/regionModels/surfaceFilmModels/noFilm/noFilm.H

@@ -116,6 +116,9 @@ public:
             //- Return the film thickness [m]
             virtual const volScalarField& delta() const;
 
+            //- Return const access to the surface tension / [m/s2]
+            inline const volScalarField& sigma() const;
+
             //- Return the film velocity [m/s]
             virtual const volVectorField& U() const;
 
@@ -146,14 +149,14 @@ public:
 
             // Transfer fields - to the primary region
 
+                //- Return mass transfer source - Eulerian phase only
+                virtual tmp<volScalarField> primaryMassTrans() const;
+
                 //- Return the film mass available for transfer
-               virtual const volScalarField& massForPrimary() const;
+                virtual const volScalarField& cloudMassTrans() const;
 
                 //- Return the parcel diameters originating from film
-               virtual const volScalarField& diametersForPrimary() const;
+                virtual const volScalarField& cloudDiameterTrans() const;
-
-               //- Return the film mass evolved via phase change
-               virtual const volScalarField& massPhaseChangeForPrimary() const;
 
 
         // Source fields

src/regionModels/surfaceFilmModels/submodels/kinematic/injectionModel/curvatureSeparation/curvatureSeparation.C

@@ -0,0 +1,355 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2011 OpenCFD Ltd.
+     \\/     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 2 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, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+\*---------------------------------------------------------------------------*/
+
+#include "curvatureSeparation.H"
+#include "addToRunTimeSelectionTable.H"
+#include "fvMesh.H"
+#include "Time.H"
+#include "volFields.H"
+#include "kinematicSingleLayer.H"
+#include "surfaceInterpolate.H"
+#include "fvcDiv.H"
+#include "fvcGrad.H"
+#include "stringListOps.H"
+#include "cyclicPolyPatch.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace regionModels
+{
+namespace surfaceFilmModels
+{
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+defineTypeNameAndDebug(curvatureSeparation, 0);
+addToRunTimeSelectionTable
+(
+    injectionModel,
+    curvatureSeparation,
+    dictionary
+);
+
+// * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
+
+tmp<volScalarField> curvatureSeparation::calcInvR1
+(
+    const volVectorField& U
+) const
+{
+    // method 1
+/*
+    tmp<volScalarField> tinvR1
+    (
+        new volScalarField("invR1", fvc::div(owner().nHat()))
+    );
+*/
+
+    // method 2
+    dimensionedScalar smallU("smallU", dimVelocity, ROOTVSMALL);
+    volVectorField UHat(U/(mag(U) + smallU));
+    tmp<volScalarField> tinvR1
+    (
+        new volScalarField("invR1", UHat & (UHat & gradNHat_))
+    );
+
+
+    scalarField& invR1 = tinvR1().internalField();
+
+    // apply defined patch radii
+    const scalar rMin = 1e-6;
+    const fvMesh& mesh = owner().regionMesh();
+    const polyBoundaryMesh& pbm = mesh.boundaryMesh();
+    forAll(definedPatchRadii_, i)
+    {
+        label patchI = definedPatchRadii_[i].first();
+        scalar definedInvR1 = 1.0/max(rMin, definedPatchRadii_[i].second());
+        UIndirectList<scalar>(invR1, pbm[patchI].faceCells()) = definedInvR1;
+    }
+
+    // filter out large radii
+    const scalar rMax = 1e6;
+    forAll(invR1, i)
+    {
+        if (mag(invR1[i]) < 1/rMax)
+        {
+            invR1[i] = -1.0;
+        }
+    }
+
+    if (debug && mesh.time().outputTime())
+    {
+        tinvR1().write();
+    }
+
+    return tinvR1;
+}
+
+
+tmp<scalarField> curvatureSeparation::calcCosAngle
+(
+    const surfaceScalarField& phi
+) const
+{
+    const fvMesh& mesh = owner().regionMesh();
+    const vectorField nf(mesh.Sf()/mesh.magSf());
+    const unallocLabelList& own = mesh.owner();
+    const unallocLabelList& nbr = mesh.neighbour();
+
+    scalarField phiMax(mesh.nCells(), -GREAT);
+    scalarField cosAngle(mesh.nCells(), 0.0);
+    forAll(nbr, faceI)
+    {
+        label cellO = own[faceI];
+        label cellN = nbr[faceI];
+
+        if (phi[faceI] > phiMax[cellO])
+        {
+            phiMax[cellO] = phi[faceI];
+            cosAngle[cellO] = -gHat_ & nf[faceI];
+        }
+        if (-phi[faceI] > phiMax[cellN])
+        {
+            phiMax[cellN] = -phi[faceI];
+            cosAngle[cellN] = -gHat_ & -nf[faceI];
+        }
+    }
+
+    forAll(phi.boundaryField(), patchI)
+    {
+        const fvsPatchScalarField& phip = phi.boundaryField()[patchI];
+        const fvPatch& pp = phip.patch();
+        const labelList& faceCells = pp.faceCells();
+        const vectorField nf(pp.nf());
+        forAll(phip, i)
+        {
+            label cellI = faceCells[i];
+            if (phip[i] > phiMax[cellI])
+            {
+                phiMax[cellI] = phip[i];
+                cosAngle[cellI] = -gHat_ & nf[i];
+            }
+        }
+    }
+/*
+    // correction for cyclics - use cyclic pairs' face normal instead of
+    // local face normal
+    const fvBoundaryMesh& pbm = mesh.boundary();
+    forAll(phi.boundaryField(), patchI)
+    {
+        if (isA<cyclicPolyPatch>(pbm[patchI]))
+        {
+            const scalarField& phip = phi.boundaryField()[patchI];
+            const vectorField nf(pbm[patchI].nf());
+            const labelList& faceCells = pbm[patchI].faceCells();
+            const label sizeBy2 = pbm[patchI].size()/2;
+
+            for (label face0=0; face0<sizeBy2; face0++)
+            {
+                label face1 = face0 + sizeBy2;
+                label cell0 = faceCells[face0];
+                label cell1 = faceCells[face1];
+
+                // flux leaving half 0, entering half 1
+                if (phip[face0] > phiMax[cell0])
+                {
+                    phiMax[cell0] = phip[face0];
+                    cosAngle[cell0] = -gHat_ & -nf[face1];
+                }
+
+                // flux leaving half 1, entering half 0
+                if (-phip[face1] > phiMax[cell1])
+                {
+                    phiMax[cell1] = -phip[face1];
+                    cosAngle[cell1] = -gHat_ & nf[face0];
+                }
+            }
+        }
+    }
+*/
+    // checks
+    if (debug && mesh.time().outputTime())
+    {
+        volScalarField volCosAngle
+        (
+            IOobject
+            (
+                "cosAngle",
+                mesh.time().timeName(),
+                mesh,
+                IOobject::NO_READ
+            ),
+            mesh,
+            dimensionedScalar("zero", dimless, 0.0),
+            zeroGradientFvPatchScalarField::typeName
+        );
+        volCosAngle.internalField() = cosAngle;
+        volCosAngle.correctBoundaryConditions();
+        volCosAngle.write();
+    }
+
+    return max(min(cosAngle, 1.0), -1.0);
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+curvatureSeparation::curvatureSeparation
+(
+    const surfaceFilmModel& owner,
+    const dictionary& dict
+)
+:
+    injectionModel(type(), owner, dict),
+    gradNHat_(fvc::grad(owner.nHat())),
+    deltaByR1Min_(coeffs().lookupOrDefault<scalar>("deltaByR1Min", 0.0)),
+    definedPatchRadii_(),
+    magG_(mag(owner.g().value())),
+    gHat_(owner.g().value()/magG_)
+{
+    List<Tuple2<word, scalar> > prIn(coeffs().lookup("definedPatchRadii"));
+    const wordList& allPatchNames = owner.regionMesh().boundaryMesh().names();
+
+    DynamicList<Tuple2<label, scalar> > prData(allPatchNames.size());
+
+    labelHashSet uniquePatchIDs;
+
+    forAllReverse(prIn, i)
+    {
+        labelList patchIDs = findStrings(prIn[i].first(), allPatchNames);
+        forAll(patchIDs, j)
+        {
+            const label patchI = patchIDs[j];
+
+            if (!uniquePatchIDs.found(patchI))
+            {
+                const scalar radius = prIn[i].second();
+                prData.append(Tuple2<label, scalar>(patchI, radius));
+
+                uniquePatchIDs.insert(patchI);
+            }
+        }
+    }
+
+    definedPatchRadii_.transfer(prData);
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+curvatureSeparation::~curvatureSeparation()
+{}
+
+
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+
+void curvatureSeparation::correct
+(
+    scalarField& availableMass,
+    scalarField& massToInject,
+    scalarField& diameterToInject
+)
+{
+    const kinematicSingleLayer& film =
+        refCast<const kinematicSingleLayer>(this->owner());
+    const fvMesh& mesh = film.regionMesh();
+
+    const volScalarField& delta = film.delta();
+    const volVectorField& U = film.U();
+    const surfaceScalarField& phi = film.phi();
+    const volScalarField& rho = film.rho();
+    const scalarField magSqrU(magSqr(film.U()));
+    const volScalarField& sigma = film.sigma();
+
+    const scalarField invR1(calcInvR1(U));
+    const scalarField cosAngle(calcCosAngle(phi));
+
+    // calculate force balance
+    const scalar Fthreshold = 1e-10;
+    scalarField Fnet(mesh.nCells(), 0.0);
+    scalarField separated(mesh.nCells(), 0.0);
+    forAll(invR1, i)
+    {
+        if ((invR1[i] > 0) && (delta[i]*invR1[i] > deltaByR1Min_))
+        {
+            scalar R1 = 1.0/(invR1[i] + ROOTVSMALL);
+            scalar R2 = R1 + delta[i];
+
+            // inertial force
+            scalar Fi = -delta[i]*rho[i]*magSqrU[i]*72.0/60.0*invR1[i];
+
+            // body force
+            scalar Fb =
+              - 0.5*rho[i]*magG_*invR1[i]*(sqr(R1) - sqr(R2))*cosAngle[i];
+
+            // surface force
+            scalar Fs = sigma[i]/R2;
+
+            Fnet[i] = Fi + Fb + Fs;
+
+            if (Fnet[i] + Fthreshold < 0)
+            {
+                separated[i] = 1.0;
+            }
+        }
+    }
+
+    // inject all available mass
+    massToInject = separated*availableMass;
+    diameterToInject = separated*delta;
+    availableMass -= separated*availableMass;
+
+    if (debug && mesh.time().outputTime())
+    {
+        volScalarField volFnet
+        (
+            IOobject
+            (
+                "Fnet",
+                mesh.time().timeName(),
+                mesh,
+                IOobject::NO_READ
+            ),
+            mesh,
+            dimensionedScalar("zero", dimForce, 0.0),
+            zeroGradientFvPatchScalarField::typeName
+        );
+        volFnet.internalField() = Fnet;
+        volFnet.correctBoundaryConditions();
+        volFnet.write();
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace surfaceFilmModels
+} // End namespace regionModels
+} // End namespace Foam
+
+// ************************************************************************* //

src/regionModels/surfaceFilmModels/submodels/kinematic/injectionModel/curvatureSeparation/curvatureSeparation.H

@@ -0,0 +1,158 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2011 OpenCFD Ltd.
+     \\/     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 2 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, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Class
+    Foam::curvatureSeparation
+
+Description
+    Curvature film separation model
+
+    Assesses film curvature via the mesh geometry and calculates a force
+    balance of the form:
+
+        F_sum = F_inertial + F_body + F_surface
+
+    If F_sum < 0, the film separates. Similarly, if F_sum > 0 the film will
+    remain attached.
+
+    Based on description given by
+        Owen and D. J. Ryley. The flow of thin liquid films around corners.
+        International Journal of Multiphase Flow, 11(1):51-62, 1985.
+
+
+SourceFiles
+    curvatureSeparation.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef curvatureSeparation_H
+#define curvatureSeparation_H
+
+#include "injectionModel.H"
+#include "surfaceFields.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace regionModels
+{
+namespace surfaceFilmModels
+{
+
+/*---------------------------------------------------------------------------*\
+                   Class curvatureSeparation Declaration
+\*---------------------------------------------------------------------------*/
+
+class curvatureSeparation
+:
+    public injectionModel
+{
+private:
+
+    // Private member functions
+
+        //- Disallow default bitwise copy construct
+        curvatureSeparation(const curvatureSeparation&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const curvatureSeparation&);
+
+
+protected:
+
+    // Protected data
+
+        //- Gradient of surface normals
+        volTensorField gradNHat_;
+
+        //- Minimum gravity driven film thickness (non-dimensionalised delta/R1)
+        scalar deltaByR1Min_;
+
+        //- List of radii for patches - if patch not defined, radius
+        // calculated based on mesh geometry
+        List<Tuple2<label, scalar> > definedPatchRadii_;
+
+        //- Magnitude of gravity vector
+        scalar magG_;
+
+        //- Direction of gravity vector
+        vector gHat_;
+
+
+    // Protected Member Functions
+
+        //- Calculate local (inverse) radius of curvature
+        tmp<volScalarField> calcInvR1(const volVectorField& U) const;
+
+        //- Calculate the cosine of the angle between gravity vector and
+        //  cell out flow direction 
+        tmp<scalarField> calcCosAngle(const surfaceScalarField& phi) const;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("curvatureSeparation");
+
+
+    // Constructors
+
+        //- Construct from surface film model
+        curvatureSeparation
+        (
+            const surfaceFilmModel& owner,
+            const dictionary& dict
+        );
+
+
+    //- Destructor
+    virtual ~curvatureSeparation();
+
+
+    // Member Functions
+
+        // Evolution
+
+            //- Correct
+            virtual void correct
+            (
+                scalarField& availableMass,
+                scalarField& massToInject,
+                scalarField& diameterToInject
+            );
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace surfaceFilmModels
+} // End namespace regionModels
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/regionModels/surfaceFilmModels/submodels/kinematic/injectionModel/drippingInjection/drippingInjection.C

@@ -24,13 +24,14 @@ License
 
 \*---------------------------------------------------------------------------*/
 
-#include "cloudInjection.H"
+#include "drippingInjection.H"
 #include "addToRunTimeSelectionTable.H"
 #include "fvMesh.H"
 #include "Time.H"
 #include "mathematicalConstants.H"
 #include "Random.H"
 #include "volFields.H"
+#include "kinematicSingleLayer.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -43,18 +44,19 @@ namespace surfaceFilmModels
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
-defineTypeNameAndDebug(cloudInjection, 0);
+defineTypeNameAndDebug(drippingInjection, 0);
-addToRunTimeSelectionTable(injectionModel, cloudInjection, dictionary);
+addToRunTimeSelectionTable(injectionModel, drippingInjection, dictionary);
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
-cloudInjection::cloudInjection
+drippingInjection::drippingInjection
 (
     const surfaceFilmModel& owner,
     const dictionary& dict
 )
 :
     injectionModel(type(), owner, dict),
+    deltaStable_(readScalar(coeffs_.lookup("deltaStable"))),
     particlesPerParcel_(readScalar(coeffs_.lookup("particlesPerParcel"))),
     rndGen_(label(0), -1),
     parcelDistribution_
@@ -76,31 +78,59 @@ cloudInjection::cloudInjection
 
 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 
-cloudInjection::~cloudInjection()
+drippingInjection::~drippingInjection()
 {}
 
 
 // * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
 
-void cloudInjection::correct
+void drippingInjection::correct
 (
+    scalarField& availableMass,
     scalarField& massToInject,
     scalarField& diameterToInject
 )
 {
+    const kinematicSingleLayer& film =
+        refCast<const kinematicSingleLayer>(this->owner());
+
     const scalar pi = constant::mathematical::pi;
-    const scalarField& rhoFilm = owner().rho();
+
+    // calculate available dripping mass
+    tmp<volScalarField> tgNorm(film.gNorm());
+    const scalarField& gNorm = tgNorm();
+    const scalarField& magSf = film.magSf();
+
+    const scalarField& delta = film.delta();
+    const scalarField& rho = film.rho();
+
+    scalarField massDrip(film.regionMesh().nCells(), 0.0);
+
+    forAll(gNorm, i)
+    {
+        if (gNorm[i] > SMALL)
+        {
+            const scalar ddelta = max(0.0, delta[i] - deltaStable_);
+            massDrip[i] +=
+                min(availableMass[i], max(0.0, ddelta*rho[i]*magSf[i]));
+        }
+    }
+
 
     // Collect the data to be transferred
     forAll(massToInject, cellI)
     {
-        scalar rho = rhoFilm[cellI];
+        scalar rhoc = rho[cellI];
         scalar diam = diameter_[cellI];
-        scalar minMass = particlesPerParcel_*rho*pi/6*pow3(diam);
+        scalar minMass = particlesPerParcel_*rhoc*pi/6*pow3(diam);
 
-        if (massToInject[cellI] > minMass)
+        if (massDrip[cellI] > minMass)
         {
-            // All mass can be injected - set particle diameter
+            // All drip mass can be injected
+            massToInject[cellI] += massDrip[cellI];
+            availableMass[cellI] -= massDrip[cellI];
+
+            // Set particle diameter
             diameterToInject[cellI] = diameter_[cellI];
 
             // Retrieve new particle diameter sample

src/regionModels/surfaceFilmModels/submodels/kinematic/injectionModel/drippingInjection/drippingInjection.H

@@ -23,18 +23,23 @@ License
     Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
 
 Class
-    Foam::cloudInjection
+    Foam::drippingInjection
 
 Description
-    Cloud injection model
+    Film Dripping mass transfer model.
+
+    If the film mass exceeds that needed to generate a valid parcel, the
+    equivalent mass is removed from the film.
+
+    New parcel diameters are sampled from a PDF.
 
 SourceFiles
-    cloudInjection.C
+    drippingInjection.C
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef cloudInjection_H
+#ifndef drippingInjection_H
-#define cloudInjection_H
+#define drippingInjection_H
 
 #include "injectionModel.H"
 #include "distributionModel.H"
@@ -50,10 +55,10 @@ namespace surfaceFilmModels
 {
 
 /*---------------------------------------------------------------------------*\
-                        Class cloudInjection Declaration
+                    Class drippingInjection Declaration
 \*---------------------------------------------------------------------------*/
 
-class cloudInjection
+class drippingInjection
 :
     public injectionModel
 {
@@ -62,16 +67,20 @@ private:
     // Private member functions
 
         //- Disallow default bitwise copy construct
-        cloudInjection(const cloudInjection&);
+        drippingInjection(const drippingInjection&);
 
         //- Disallow default bitwise assignment
-        void operator=(const cloudInjection&);
+        void operator=(const drippingInjection&);
 
 
 protected:
 
     // Protected data
 
+        //- Stable film thickness - drips only formed if thickness
+        //  execeeds this threhold value
+        scalar deltaStable_;
+
         //- Number of particles per parcel
         scalar particlesPerParcel_;
 
@@ -82,24 +91,28 @@ protected:
         const autoPtr<distributionModels::distributionModel>
             parcelDistribution_;
 
-        //- Diameters of particles to inject into the cloud
+        //- Diameters of particles to inject into the dripping
         scalarList diameter_;
 
 
 public:
 
     //- Runtime type information
-    TypeName("cloudInjection");
+    TypeName("drippingInjection");
 
 
     // Constructors
 
         //- Construct from surface film model
-        cloudInjection(const surfaceFilmModel& owner, const dictionary& dict);
+        drippingInjection
+        (
+            const surfaceFilmModel& owner,
+            const dictionary& dict
+        );
 
 
     //- Destructor
-    virtual ~cloudInjection();
+    virtual ~drippingInjection();
 
 
     // Member Functions
@@ -109,6 +122,7 @@ public:
             //- Correct
             virtual void correct
             (
+                scalarField& availableMass,
                 scalarField& massToInject,
                 scalarField& diameterToInject
             );

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

@@ -26,12 +26,12 @@ Class
     Foam::injectionModel
 
 Description
-    Base class for film injection models
+    Base class for film injection models, handling mass transfer from the
+    film.
 
 SourceFiles
     injectionModel.C
     injectionModelNew.C
-
 \*---------------------------------------------------------------------------*/
 
 #ifndef injectionModel_H
@@ -109,7 +109,8 @@ public:
         static autoPtr<injectionModel> New
         (
             const surfaceFilmModel& owner,
-            const dictionary& dict
+            const dictionary& dict,
+            const word& mdoelType
         );
 
 
@@ -124,6 +125,7 @@ public:
             //- Correct
             virtual void correct
             (
+                scalarField& availableMass,
                 scalarField& massToInject,
                 scalarField& diameterToInject
             ) = 0;

src/regionModels/surfaceFilmModels/submodels/kinematic/injectionModel/injectionModel/injectionModelNew.C

@@ -40,12 +40,11 @@ namespace surfaceFilmModels
 autoPtr<injectionModel> injectionModel::New
 (
     const surfaceFilmModel& model,
-    const dictionary& dict
+    const dictionary& dict,
+    const word& modelType
 )
 {
-    word modelType(dict.lookup("injectionModel"));
+    Info<< "        " << modelType << endl;
-
-    Info<< "    Selecting injectionModel " << modelType << endl;
 
     dictionaryConstructorTable::iterator cstrIter =
         dictionaryConstructorTablePtr_->find(modelType);

src/regionModels/surfaceFilmModels/submodels/kinematic/injectionModel/injectionModelList/injectionModelList.C

@@ -0,0 +1,135 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2011 OpenCFD Ltd.
+     \\/     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 2 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, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+\*---------------------------------------------------------------------------*/
+
+#include "injectionModelList.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace regionModels
+{
+namespace surfaceFilmModels
+{
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+injectionModelList::injectionModelList(const surfaceFilmModel& owner)
+:
+    PtrList<injectionModel>(),
+    owner_(owner),
+    dict_(dictionary::null),
+    injectedMassTotal_(0.0)
+{}
+
+
+injectionModelList::injectionModelList
+(
+    const surfaceFilmModel& owner,
+    const dictionary& dict
+)
+:
+    PtrList<injectionModel>(),
+    owner_(owner),
+    dict_(dict),
+    injectedMassTotal_(0.0)
+{
+    const wordList activeModels(dict.lookup("injectionModels"));
+
+    wordHashSet models;
+    forAll(activeModels, i)
+    {
+        models.insert(activeModels[i]);
+    }
+
+    Info<< "    Selecting film injection models" << endl;
+    if (models.size() > 0)
+    {
+        this->setSize(models.size());
+
+        label i = 0;
+        forAllConstIter(wordHashSet, models, iter)
+        {
+            const word& model = iter.key();
+            set
+            (
+                i,
+                injectionModel::New(owner, dict, model)
+            );
+            i++;
+        }
+    }
+    else
+    {
+        Info<< "        none" << endl;
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+injectionModelList::~injectionModelList()
+{}
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void injectionModelList::correct
+(
+    scalarField& availableMass,
+    volScalarField& massToInject,
+    volScalarField& diameterToInject
+)
+{
+    // Correct models that accumulate mass and diameter transfers
+    forAll(*this, i)
+    {
+        injectionModel& im = operator[](i);
+        im.correct(availableMass, massToInject, diameterToInject);
+    }
+
+    injectedMassTotal_ += sum(massToInject.internalField());
+
+
+    // Push values to boundaries ready for transfer to the primary region
+    massToInject.correctBoundaryConditions();
+    diameterToInject.correctBoundaryConditions();
+}
+
+
+void injectionModelList::info(Ostream& os) const
+{
+    os  << indent << "injected mass      = "
+        << returnReduce<scalar>(injectedMassTotal_, sumOp<scalar>()) << nl;
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace surfaceFilmModels
+} // End namespace regionModels
+} // End namespace Foam
+
+// ************************************************************************* //

src/regionModels/surfaceFilmModels/submodels/kinematic/injectionModel/injectionModelList/injectionModelList.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2009-2011 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2011-2011 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -23,19 +23,20 @@ License
     Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
 
 Class
-    Foam::noInjection
+    Foam::injectionModelList
 
 Description
-    Dummy injection model for 'none'
+    List container for film injection models
 
 SourceFiles
-    noInjection.C
+    injectionModelList.C
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef noInjection_H
+#ifndef injectionModelList_H
-#define noInjection_H
+#define injectionModelList_H
 
+#include "PtrList.H"
 #include "injectionModel.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -48,38 +49,53 @@ namespace surfaceFilmModels
 {
 
 /*---------------------------------------------------------------------------*\
-                        Class noInjection Declaration
+                    Class injectionModelList Declaration
 \*---------------------------------------------------------------------------*/
 
-class noInjection
+class injectionModelList
 :
-    public injectionModel
+    public PtrList<injectionModel>
 {
 private:
 
-    // Private member functions
+    // Private data
+
+        //- Reference to the owner surface film model
+        const surfaceFilmModel& owner_;
+
+        //- Dictionary
+        dictionary dict_;
+
+        //- Cumulative mass injected total
+        scalar injectedMassTotal_;
+
+
+    // Private Member Functions
 
         //- Disallow default bitwise copy construct
-        noInjection(const noInjection&);
+        injectionModelList(const injectionModelList&);
 
         //- Disallow default bitwise assignment
-        void operator=(const noInjection&);
+        void operator=(const injectionModelList&);
 
 
 public:
 
-    //- Runtime type information
-    TypeName("none");
-
-
     // Constructors
 
-        //- Construct from surface film model
+        //- Construct null
-        noInjection(const surfaceFilmModel& owner, const dictionary& dict);
+        injectionModelList(const surfaceFilmModel& owner);
+
+        //- Construct from type name, dictionary and surface film model
+        injectionModelList
+        (
+            const surfaceFilmModel& owner,
+            const dictionary& dict
+        );
 
 
     //- Destructor
-    virtual ~noInjection();
+    virtual ~injectionModelList();
 
 
     // Member Functions
@@ -89,9 +105,16 @@ public:
             //- Correct
             virtual void correct
             (
-                scalarField& massToInject,
+                scalarField& availableMass,
-                scalarField& diameterToInject
+                volScalarField& massToInject,
+                volScalarField& diameterToInject
             );
+
+
+        // I-O
+
+            //- Provide some info
+            void info(Ostream& os) const;
 };
 
 

src/regionModels/surfaceFilmModels/submodels/kinematic/injectionModel/noInjection/noInjection.C

@@ -1,81 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2009-2011 OpenCFD Ltd.
-     \\/     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 2 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, write to the Free Software Foundation,
-    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
-
-\*---------------------------------------------------------------------------*/
-
-#include "noInjection.H"
-#include "addToRunTimeSelectionTable.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-namespace regionModels
-{
-namespace surfaceFilmModels
-{
-
-// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
-
-defineTypeNameAndDebug(noInjection, 0);
-addToRunTimeSelectionTable(injectionModel, noInjection, dictionary);
-
-// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
-
-noInjection::noInjection
-(
-    const surfaceFilmModel& owner,
-    const dictionary&
-)
-:
-    injectionModel(owner)
-{}
-
-
-// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
-
-noInjection::~noInjection()
-{}
-
-
-// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
-
-void noInjection::correct
-(
-    scalarField& massToInject,
-    scalarField& diameterToInject
-)
-{
-    // no mass injected
-    massToInject = 0.0;
-}
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace surfaceFilmModels
-} // End namespace regionModels
-} // End namespace Foam
-
-// ************************************************************************* //

src/regionModels/surfaceFilmModels/submodels/kinematic/injectionModel/removeInjection/removeInjection.C

@@ -63,11 +63,13 @@ removeInjection::~removeInjection()
 
 void removeInjection::correct
 (
-    scalarField&,
+    scalarField& availableMass,
+    scalarField& massToInject,
     scalarField&
 )
 {
-    // do nothing - all mass available to be removed
+    massToInject = availableMass;
+    availableMass = 0.0;
 }
 
 

src/regionModels/surfaceFilmModels/submodels/kinematic/injectionModel/removeInjection/removeInjection.H

@@ -89,6 +89,7 @@ public:
             //- Correct
             virtual void correct
             (
+                scalarField& availableMass,
                 scalarField& massToInject,
                 scalarField& diameterToInject
             );

src/regionModels/surfaceFilmModels/submodels/subModelBaseI.H

@@ -24,8 +24,6 @@ License
 
 \*---------------------------------------------------------------------------*/
 
-#include "subModelBase.H"
-
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 namespace Foam

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

@@ -61,10 +61,11 @@ noPhaseChange::~noPhaseChange()
 
 // * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
 
-void noPhaseChange::correct
+void noPhaseChange::correctModel
 (
     const scalar,
     scalarField&,
+    scalarField&,
     scalarField&
 )
 {

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

@@ -87,9 +87,10 @@ public:
         // Evolution
 
             //- Correct
-            virtual void correct
+            virtual void correctModel
             (
                 const scalar dt,
+                scalarField& availableMass,
                 scalarField& dMass,
                 scalarField& dEnergy
             );

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

@@ -47,7 +47,9 @@ phaseChangeModel::phaseChangeModel
     const surfaceFilmModel& owner
 )
 :
-    subModelBase(owner)
+    subModelBase(owner),
+    latestMassPC_(0.0),
+    totalMassPC_(0.0)
 {}
 
 
@@ -58,7 +60,9 @@ phaseChangeModel::phaseChangeModel
     const dictionary& dict
 )
 :
-    subModelBase(type, owner, dict)
+    subModelBase(type, owner, dict),
+    latestMassPC_(0.0),
+    totalMassPC_(0.0)
 {}
 
 
@@ -68,6 +72,44 @@ phaseChangeModel::~phaseChangeModel()
 {}
 
 
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+
+void phaseChangeModel::correct
+(
+    const scalar dt,
+    scalarField& availableMass,
+    volScalarField& dMass,
+    volScalarField& dEnergy
+)
+{
+    correctModel
+    (
+        dt,
+        availableMass,
+        dMass,
+        dEnergy
+    );
+
+    latestMassPC_ = sum(dMass.internalField());
+    totalMassPC_ += latestMassPC_; 
+
+    availableMass -= dMass;
+    dMass.correctBoundaryConditions();
+}
+
+
+void phaseChangeModel::info(Ostream& os) const
+{
+    const scalar massPCRate =
+        returnReduce(latestMassPC_, sumOp<scalar>())
+       /owner_.time().deltaTValue();
+
+    os  << indent << "mass phase change  = "
+        << returnReduce(totalMassPC_, sumOp<scalar>()) << nl
+        << indent << "vapourisation rate = " << massPCRate << nl;
+}
+
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 } // end namespace surfaceFilmModels

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

@@ -69,6 +69,17 @@ private:
         void operator=(const phaseChangeModel&);
 
 
+protected:
+
+    // Protected Member Functions
+
+        //- Latest mass transfer due to phase change
+        scalar latestMassPC_;
+
+        //- Total mass transfer due to phase change
+        scalar totalMassPC_;
+
+
 public:
 
     //- Runtime type information
@@ -125,9 +136,25 @@ public:
             virtual void correct
             (
                 const scalar dt,
+                scalarField& availableMass,
+                volScalarField& dMass,
+                volScalarField& dEnergy
+            );
+
+            //- Correct
+            virtual void correctModel
+            (
+                const scalar dt,
+                scalarField& availableMass,
                 scalarField& dMass,
                 scalarField& dEnergy
             ) = 0;
+
+
+        // I-O
+
+            //- Provide some feedback
+            virtual void info(Ostream& os) const;
 };
 
 

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

@@ -81,9 +81,7 @@ standardPhaseChange::standardPhaseChange
     Tb_(readScalar(coeffs_.lookup("Tb"))),
     deltaMin_(readScalar(coeffs_.lookup("deltaMin"))),
     L_(readScalar(coeffs_.lookup("L"))),
-    TbFactor_(coeffs_.lookupOrDefault<scalar>("TbFactor", 1.1)),
+    TbFactor_(coeffs_.lookupOrDefault<scalar>("TbFactor", 1.1))
-    totalMass_(0.0),
-    vapourRate_(0.0)
 {}
 
 
@@ -95,9 +93,10 @@ standardPhaseChange::~standardPhaseChange()
 
 // * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
 
-void standardPhaseChange::correct
+void standardPhaseChange::correctModel
 (
     const scalar dt,
+    scalarField& availableMass,
     scalarField& dMass,
     scalarField& dEnergy
 )
@@ -124,8 +123,7 @@ void standardPhaseChange::correct
     const scalarField hInf(film.htcs().h());
     const scalarField hFilm(film.htcw().h());
     const vectorField dU(film.UPrimary() - film.Us());
-    const scalarField availableMass((delta - deltaMin_)*rho*magSf);
+    const scalarField limMass(max(0.0, availableMass - deltaMin_*rho*magSf));
-
 
     forAll(dMass, cellI)
     {
@@ -152,8 +150,7 @@ void standardPhaseChange::correct
 
                 const scalar Cp = liq.Cp(pc, Tloc);
                 const scalar Tcorr = max(0.0, T[cellI] - Tb_);
-                const scalar qCorr = availableMass[cellI]*Cp*(Tcorr);
+                const scalar qCorr = limMass[cellI]*Cp*(Tcorr);
-
                 dMass[cellI] =
                     dt*magSf[cellI]/hVap*(qDotInf + qDotFilm)
                   + qCorr/hVap;
@@ -195,23 +192,10 @@ void standardPhaseChange::correct
                     dt*magSf[cellI]*rhoInfc*hm*(Ys - YInf[cellI])/(1.0 - Ys);
             }
 
-            dMass[cellI] = min(availableMass[cellI], max(0.0, dMass[cellI]));
+            dMass[cellI] = min(limMass[cellI], max(0.0, dMass[cellI]));
             dEnergy[cellI] = dMass[cellI]*hVap;
         }
     }
-
-    const scalar sumdMass = sum(dMass);
-    totalMass_ += sumdMass;
-    vapourRate_ = sumdMass/owner().time().deltaTValue();
-}
-
-
-void standardPhaseChange::info() const
-{
-    Info<< indent << "mass phase change  = "
-        << returnReduce(totalMass_, sumOp<scalar>()) << nl
-        << indent << "vapourisation rate = "
-        << returnReduce(vapourRate_, sumOp<scalar>()) << nl;
 }
 
 

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

@@ -83,12 +83,6 @@ protected:
         //  Used to set max limit on temperature to Tb*TbFactor
         const scalar TbFactor_;
 
-        //- Total mass evolved / [kg]
-        scalar totalMass_;
-
-        //- Vapouristaion rate / kg/s
-        scalar vapourRate_;
-
 
     // Protected member functions
 
@@ -121,18 +115,13 @@ public:
         // Evolution
 
             //- Correct
-            virtual void correct
+            virtual void correctModel
             (
                 const scalar dt,
+                scalarField& availableMass,
                 scalarField& dMass,
                 scalarField& dEnergy
             );
-
-
-        // Input/output
-
-            //- Output model statistics
-            virtual void info() const;
 };
 
 

src/regionModels/surfaceFilmModels/surfaceFilmModel/surfaceFilmModel.H

@@ -161,6 +161,9 @@ public:
 
         // Access
 
+            //- Return the accleration due to gravity
+            inline const dimensionedVector& g() const;
+
             //- Return the thermo type
             inline const thermoModelType& thermoModel() const;
 
@@ -214,18 +217,20 @@ public:
             //- Return the film thermal conductivity [W/m/K]
             virtual const volScalarField& kappa() const = 0;
 
+            //- Return the film surface tension [N/m]
+            virtual const volScalarField& sigma() const = 0;
+
 
             // Transfer fields - to the primary region
 
+                //- Return mass transfer source - Eulerian phase only
+                virtual tmp<volScalarField> primaryMassTrans() const = 0;
+
                 //- Return the film mass available for transfer
-               virtual const volScalarField& massForPrimary() const = 0;
+                virtual const volScalarField& cloudMassTrans() const = 0;
 
                 //- Return the parcel diameters originating from film
-               virtual const volScalarField& diametersForPrimary() const = 0;
+                virtual const volScalarField& cloudDiameterTrans() const = 0;
-
-               //- Return the film mass evolved via phase change
-               virtual const volScalarField& massPhaseChangeForPrimary()
-                   const = 0;
 
 
         // Source fields

src/regionModels/surfaceFilmModels/surfaceFilmModel/surfaceFilmModelI.H

@@ -37,6 +37,12 @@ namespace surfaceFilmModels
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
+inline const Foam::dimensionedVector& surfaceFilmModel::g() const
+{
+    return g_;
+}
+
+
 inline const surfaceFilmModel::thermoModelType&
 surfaceFilmModel::thermoModel() const
 {

src/regionModels/surfaceFilmModels/thermoSingleLayer/thermoSingleLayer.C

@@ -81,6 +81,11 @@ bool thermoSingleLayer::read()
 
 void thermoSingleLayer::resetPrimaryRegionSourceTerms()
 {
+    if (debug)
+    {
+        Info<< "thermoSingleLayer::resetPrimaryRegionSourceTerms()" << endl;
+    }
+
     kinematicSingleLayer::resetPrimaryRegionSourceTerms();
 
     hsSpPrimary_ == dimensionedScalar("zero", hsSp_.dimensions(), 0.0);
@@ -105,6 +110,7 @@ void thermoSingleLayer::correctThermoFields()
         {
             const liquidProperties& liq =
                 thermo_.liquids().properties()[liquidId_];
+
             forAll(rho_, cellI)
             {
                 const scalar T = T_[cellI];
@@ -173,6 +179,11 @@ void thermoSingleLayer::updateSurfaceTemperatures()
 
 void thermoSingleLayer::transferPrimaryRegionThermoFields()
 {
+    if (debug)
+    {
+        Info<< "thermoSingleLayer::transferPrimaryRegionThermoFields()" << endl;
+    }
+
     kinematicSingleLayer::transferPrimaryRegionThermoFields();
 
     // Update primary region fields on local region via direct mapped (coupled)
@@ -187,6 +198,11 @@ void thermoSingleLayer::transferPrimaryRegionThermoFields()
 
 void thermoSingleLayer::transferPrimaryRegionSourceFields()
 {
+    if (debug)
+    {
+        Info<< "thermoSingleLayer::transferPrimaryRegionSourceFields()" << endl;
+    }
+
     kinematicSingleLayer::transferPrimaryRegionSourceFields();
 
     // Retrieve the source fields from the primary region via direct mapped
@@ -199,27 +215,30 @@ void thermoSingleLayer::transferPrimaryRegionSourceFields()
     // Note: boundary values will still have original (neat) values
     const scalar deltaT = time_.deltaTValue();
     hsSp_.field() /= magSf()*deltaT;
+
+    // Apply enthalpy source as difference between incoming and actual states
+    hsSp_ -= rhoSp_*hs_;
 }
 
 
 void thermoSingleLayer::updateSubmodels()
 {
+    if (debug)
+    {
+        Info<< "thermoSingleLayer::updateSubmodels()" << endl;
+    }
+
     // Update heat transfer coefficient sub-models
     htcs_->correct();
     htcw_->correct();
 
-    // Update phase change
-    massPhaseChangeForPrimary_.internalField() = 0.0;
-    energyPhaseChangeForPrimary_.internalField() = 0.0;
-
     phaseChange_->correct
     (
         time_.deltaTValue(),
-        massPhaseChangeForPrimary_,
+        availableMass_,
-        energyPhaseChangeForPrimary_
+        primaryMassPCTrans_,
+        primaryEnergyPCTrans_
     );
-    massPhaseChangeForPrimary_.correctBoundaryConditions();
-    totalMassPhaseChange_ += sum(massPhaseChangeForPrimary_).value();
 
     // Update radiation
     radiation_->correct();
@@ -228,14 +247,12 @@ void thermoSingleLayer::updateSubmodels()
     kinematicSingleLayer::updateSubmodels();
 
     // Update source fields
-    hsSp_ -= energyPhaseChangeForPrimary_/magSf()/time().deltaT();
+    hsSp_ += primaryEnergyPCTrans_/magSf()/time().deltaT();
+    rhoSp_ += primaryMassPCTrans_/magSf()/time().deltaT();
 }
 
 
-tmp<fvScalarMatrix> thermoSingleLayer::q
+tmp<fvScalarMatrix> thermoSingleLayer::q(volScalarField& hs) const
-(
-    volScalarField& hs
-) const
 {
     dimensionedScalar Tstd("Tstd", dimTemperature, 298.15);
 
@@ -263,10 +280,11 @@ void thermoSingleLayer::solveEnergy()
         fvm::ddt(deltaRho_, hs_)
       + fvm::div(phi_, hs_)
      ==
-        fvm::Sp(hsSp_/(hs_ + hs0), hs_)
+//      - hsSp_
+      - fvm::Sp(hsSp_/(hs_ + hs0), hs_)
       + q(hs_)
       + radiation_->Shs()
-      - fvm::Sp(massForPrimary_/magSf()/time().deltaT(), hs_)
+      - fvm::SuSp(rhoSp_, hs_)
     );
 
     correctThermoFields();
@@ -370,10 +388,38 @@ thermoSingleLayer::thermoSingleLayer
         ),
         regionMesh(),
         dimensionedScalar("zero", dimEnergy/dimMass, 0.0),
-//        T_.boundaryField().types()
         hsBoundaryTypes()
     ),
 
+    primaryMassPCTrans_
+    (
+        IOobject
+        (
+            "primaryMassPCTrans",
+            time().timeName(),
+            regionMesh(),
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        regionMesh(),
+        dimensionedScalar("zero", dimMass, 0),
+        zeroGradientFvPatchScalarField::typeName
+    ),
+    primaryEnergyPCTrans_
+    (
+        IOobject
+        (
+            "primaryEnergyPCTrans",
+            time().timeName(),
+            regionMesh(),
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        regionMesh(),
+        dimensionedScalar("zero", dimEnergy, 0),
+        zeroGradientFvPatchScalarField::typeName
+    ),
+
     hsSp_
     (
         IOobject
@@ -429,22 +475,7 @@ thermoSingleLayer::thermoSingleLayer
         heatTransferModel::New(*this, coeffs().subDict("lowerSurfaceModels"))
     ),
     phaseChange_(phaseChangeModel::New(*this, coeffs())),
-    radiation_(filmRadiationModel::New(*this, coeffs())),
+    radiation_(filmRadiationModel::New(*this, coeffs()))
-    totalMassPhaseChange_(0.0),
-    energyPhaseChangeForPrimary_
-    (
-        IOobject
-        (
-            "energyPhaseChangeForPrimary",
-            time().timeName(),
-            regionMesh(),
-            IOobject::NO_READ,
-            IOobject::NO_WRITE
-        ),
-        regionMesh(),
-        dimensionedScalar("zero", dimEnergy, 0),
-        zeroGradientFvPatchScalarField::typeName
-    )
 {
     if (thermo_.hasMultiComponentCarrier())
     {
@@ -519,24 +550,34 @@ void thermoSingleLayer::addSources
         Info<< "    energy   = " << energySource << nl << endl;
     }
 
-    hsSpPrimary_.boundaryField()[patchI][faceI] += energySource;
+    hsSpPrimary_.boundaryField()[patchI][faceI] -= energySource;
 }
 
 
 void thermoSingleLayer::preEvolveRegion()
 {
-    transferPrimaryRegionThermoFields();
+    if (debug)
+    {
+        Info<< "thermoSingleLayer::preEvolveRegion()" << endl;
+    }
 
 //    correctHsForMappedT();
 
-    correctThermoFields();
+    kinematicSingleLayer::preEvolveRegion();
 
-    transferPrimaryRegionSourceFields();
+    // Update phase change
+    primaryMassPCTrans_ == dimensionedScalar("zero", dimMass, 0.0);
+    primaryEnergyPCTrans_ == dimensionedScalar("zero", dimEnergy, 0.0);
 }
 
 
 void thermoSingleLayer::evolveRegion()
 {
+    if (debug)
+    {
+        Info<< "thermoSingleLayer::evolveRegion()" << endl;
+    }
+
     updateSubmodels();
 
     // Solve continuity for deltaRho_
@@ -617,16 +658,67 @@ const volScalarField& thermoSingleLayer::hs() const
 }
 
 
+tmp<volScalarField> thermoSingleLayer::primaryMassTrans() const
+{
+    return primaryMassPCTrans_;
+}
+
+
 void thermoSingleLayer::info() const
 {
     kinematicSingleLayer::info();
 
     Info<< indent << "min/max(T)         = " << min(T_).value() << ", "
-        << max(T_).value() << nl
+        << max(T_).value() << nl;
-        << indent << "mass phase change  = "
+
-        << returnReduce(totalMassPhaseChange_, sumOp<scalar>()) << nl
+    phaseChange_->info(Info);
-        << indent << "vapourisation rate = "
+}
-        << sum(massPhaseChangeForPrimary_).value()/time_.deltaTValue() << nl;
+
+
+tmp<DimensionedField<scalar, volMesh> > thermoSingleLayer::Srho() const
+{
+    tmp<DimensionedField<scalar, volMesh> > tSrho
+    (
+        new DimensionedField<scalar, volMesh>
+        (
+            IOobject
+            (
+                "thermoSingleLayer::Srho",
+                time().timeName(),
+                primaryMesh(),
+                IOobject::NO_READ,
+                IOobject::NO_WRITE,
+                false
+            ),
+            primaryMesh(),
+            dimensionedScalar("zero", dimMass/dimVolume/dimTime, 0.0)
+        )
+    );
+
+    scalarField& Srho = tSrho();
+    const scalarField& V = primaryMesh().V();
+    const scalar dt = time_.deltaTValue();
+
+    forAll(intCoupledPatchIDs(), i)
+    {
+        const label filmPatchI = intCoupledPatchIDs()[i];
+        const mapDistribute& distMap = mappedPatches_[filmPatchI].map();
+
+        scalarField patchMass =
+            primaryMassPCTrans_.boundaryField()[filmPatchI];
+        distMap.distribute(patchMass);
+
+        const label primaryPatchI = primaryPatchIDs()[i];
+        const unallocLabelList& cells =
+            primaryMesh().boundaryMesh()[primaryPatchI].faceCells();
+
+        forAll(patchMass, j)
+        {
+            Srho[cells[j]] = patchMass[j]/(V[cells[j]]*dt);
+        }
+    }
+
+    return tSrho;
 }
 
 
@@ -644,7 +736,7 @@ tmp<DimensionedField<scalar, volMesh> > thermoSingleLayer::Srho
         (
             IOobject
             (
-                "thermoSingleLayer::Srho(i)",
+                "thermoSingleLayer::Srho(" + Foam::name(i) + ")",
                 time_.timeName(),
                 primaryMesh(),
                 IOobject::NO_READ,
@@ -668,7 +760,7 @@ tmp<DimensionedField<scalar, volMesh> > thermoSingleLayer::Srho
             const mapDistribute& distMap = mappedPatches_[filmPatchI].map();
 
             scalarField patchMass =
-                massPhaseChangeForPrimary_.boundaryField()[filmPatchI];
+                primaryMassPCTrans_.boundaryField()[filmPatchI];
             distMap.distribute(patchMass);
 
             const label primaryPatchI = primaryPatchIDs()[i];
@@ -706,8 +798,10 @@ tmp<DimensionedField<scalar, volMesh> > thermoSingleLayer::Sh() const
         )
     );
 /*
+    phase change energy fed back into the film...
+
     scalarField& Sh = tSh();
-    const scalarField& V = mesh_.V();
+    const scalarField& V = primaryMesh().V();
     const scalar dt = time_.deltaTValue();
 
     forAll(intCoupledPatchIDs_, i)
@@ -716,14 +810,14 @@ tmp<DimensionedField<scalar, volMesh> > thermoSingleLayer::Sh() const
         const mapDistribute& distMap = mappedPatches_[filmPatchI].map();
 
         scalarField patchEnergy =
-            energyPhaseChangeForPrimary_.boundaryField()[filmPatchI];
+            primaryEnergyPCTrans_.boundaryField()[filmPatchI];
         distMap.distribute(patchEnergy);
 
         const label primaryPatchI = primaryPatchIDs()[i];
         const unallocLabelList& cells =
             primaryMesh().boundaryMesh()[primaryPatchI].faceCells();
 
-        forAll(patchMass, j)
+        forAll(patchEnergy, j)
         {
             Sh[cells[j]] += patchEnergy[j]/(V[cells[j]]*dt);
         }

src/regionModels/surfaceFilmModels/thermoSingleLayer/thermoSingleLayer.H

@@ -115,6 +115,15 @@ protected:
                 volScalarField hs_;
 
 
+            // Transfer fields - to the primary region
+
+                //- Film mass evolved via phase change
+                volScalarField primaryMassPCTrans_;
+
+                //- Film energy evolved via phase change
+                volScalarField primaryEnergyPCTrans_;
+
+
         // Source term fields
 
             // Film region - registered to the film region mesh
@@ -157,11 +166,6 @@ protected:
             //- Radiation
             autoPtr<filmRadiationModel> radiation_;
 
-            //- Total mass transferred to primary region [kg]
-            scalar totalMassPhaseChange_;
-
-            //- Film energy evolved via phase change
-            volScalarField energyPhaseChangeForPrimary_;
 
 
     // Protected member functions
@@ -256,6 +260,13 @@ public:
                 virtual const volScalarField& hs() const;
 
 
+
+            // Transfer fields - to the primary region
+
+                //- Return mass transfer source - Eulerian phase only
+                virtual tmp<volScalarField> primaryMassTrans() const;
+
+
             // Helper functions
 
                 //- Return sensible enthalpy as a function of temperature
@@ -345,6 +356,9 @@ public:
 
             // Mapped into primary region
 
+                //- Return total mass source - Eulerian phase only
+                virtual tmp<DimensionedField<scalar, volMesh> > Srho() const;
+
                 //- Return mass source for specie i - Eulerian phase only
                 virtual tmp<DimensionedField<scalar, volMesh> > Srho
                 (

src/thermophysicalModels/basic/mixtures/basicMixture/basicMixtures.C

@@ -113,6 +113,12 @@ makeBasicPolyMixture
     3
 );
 
+makeBasicPolyMixture
+(
+    pureMixture,
+    8
+);
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 } // End namespace Foam

src/thermophysicalModels/basic/rhoThermo/hRhoThermo/hRhoThermos.C

@@ -92,6 +92,13 @@ makeBasicRhoPolyThermo
     3
 );
 
+makeBasicRhoPolyThermo
+(
+    hRhoThermo,
+    pureMixture,
+    8
+);
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 } // End namespace Foam

tutorials/combustion/PDRFoam/flamePropagationWithObstacles/constant/polyMesh/blockMeshDict

@@ -167,9 +167,12 @@ blocks
     hex (29 13 12 28 73 61 60 72) (5 5 7) simpleGrading (1 1 2.985984)
 );
 
-patches
+boundary
 (
-    patch outer
+    outer
+    {
+        type    patch;
+        faces
         (
             (91 90 86 87)
             (90 89 85 86)
@@ -204,8 +207,13 @@ patches
             (74 62 61 73)
             (89 73 72 88)
             (73 61 60 72)
-    )
+        );
-    wall ground
+    }
+
+    ground
+    {
+        type wall;
+        faces
         (
             (0 4 5 1)
             (1 5 6 2)
@@ -228,15 +236,32 @@ patches
             (15 14 62 63)
             (14 13 61 62)
             (13 12 60 61)
-    )
+        );
-    wall blockedFaces
+    }
-    ()
+
-    wall baffleWall
+    blockedFaces
-    ()
+    {
-    cycic baffleCyclic_half0
+        type wall;
-    ()
+        faces ();
-    cycic baffleCyclic_half1
+    }
-    ()
+
-
+    baffleWall
-
+    {
+        type wall;
+        faces ();
+    }
+
+    baffleCyclic_half0
+    {
+        type cyclic;
+        neighbourPatch  baffleCyclic_half1;
+        faces ();
+    }
+
+    baffleCyclic_half1
+    {
+        type cyclic;
+        neighbourPatch  baffleCyclic_half0;
+        faces ();
+    }
 );

tutorials/incompressible/simpleFoam/motorBike/system/controlDict

@@ -16,7 +16,7 @@ FoamFile
 
 application     simpleFoam;
 
-startFrom       startTime;
+startFrom       latestTime;
 
 startTime       0;
 
@@ -36,7 +36,7 @@ writeFormat     ascii;
 
 writePrecision  6;
 
-writeCompression off;
+writeCompression compressed;
 
 timeFormat      general;
 
@@ -52,59 +52,10 @@ libs
 
 functions
 {
-    // Make sure all fields for functionObjects are loaded. Prevents any
+    #include "readFields"
-    // problems running with execFlowFunctionObjects.
+    #include "streamLines"
-    readFields
+    #include "cuttingPlane"
-    {
+    #include "forceCoeffs"
-        // Where to load it from (if not already in solver)
-        functionObjectLibs ("libfieldFunctionObjects.so");
-
-        type            readFields;
-        fields          (p U k);
-    }
-
-    streamLines
-    {
-        type            streamLine;
-
-        // Output every
-        outputControl   outputTime;
-        // outputInterval 10;
-
-        setFormat       vtk; //gnuplot; //xmgr; //raw; //jplot;
-
-        // Velocity field to use for tracking.
-        U U;
-
-        // Tracked forwards (+U) or backwards (-U)
-        trackForward    true;
-
-        // Names of fields to sample. Should contain above velocity field!
-        fields (p U k);
-
-        // Steps particles can travel before being removed
-        lifeTime        10000;
-
-        // Number of steps per cell (estimate). Set to 1 to disable subcycling.
-        nSubCycle 5;
-
-        // Cloud name to use
-        cloudName       particleTracks;
-
-        // Seeding method. See the sampleSets in sampleDict.
-        seedSampleSet   uniform;  //cloud;//triSurfaceMeshPointSet;
-
-        uniformCoeffs
-        {
-            type        uniform;
-            axis        x;  //distance;
-
-            // Note: tracks slightly offset so as not to be on a face
-            start       (-1.001 1e-7 0.0011);
-            end         (-1.001 1e-7 1.0011);
-            nPoints     20;
-        }
-    }
 }
 
 

tutorials/incompressible/simpleFoam/motorBike/system/cuttingPlane

@@ -0,0 +1,37 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+
+cuttingPlane
+{
+    type            surfaces;
+    functionObjectLibs ("libsampling.so");
+    outputControl   outputTime;
+
+    surfaceFormat   vtk;
+    fields          ( p U );
+
+    interpolationScheme cellPoint;
+
+    surfaces
+    (
+        yNormal
+        {
+            type            cuttingPlane;
+            planeType       pointAndNormal;
+            pointAndNormalDict
+            {
+                basePoint       (0 0 0);
+                normalVector    (0 1 0);
+            }
+            interpolate     true;
+        }
+    );
+}
+
+
+// ************************************************************************* //

tutorials/incompressible/simpleFoam/motorBike/system/forceCoeffs

@@ -5,35 +5,28 @@
 |   \\  /    A nd           | Web:      www.OpenFOAM.com                      |
 |    \\/     M anipulation  |                                                 |
 \*---------------------------------------------------------------------------*/
-FoamFile
-{
-    version     2.0;
-    format      ascii;
-    class       volScalarField;
-    location    "0";
-    object      htcConv;
-}
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-dimensions      [1 0 -1 1 0 0 0];
+forces
+{
+    type        forceCoeffs;
+    functionObjectLibs ( "libforces.so" );
+    outputControl timeStep;
+    outputInterval 1;
 
-internalField   uniform 0;
+    patches     ( "motorBike.*" );
-
+    pName       p;
-boundaryField
+    UName       U;
-{
+    rhoName     rhoInf;      // Indicates incompressible
-    "(.*)"
+    log         true;
-    {
+    rhoInf      1;           // Redundant for incompressible
-        type            calculated;
+    liftDir     (0 0 1);
-        value           uniform 0;
+    dragDir     (1 0 0);
-    }
+    CofR        (0.72 0 0);  // Axle midpoint on ground
-    region0_to_wallFilmRegion_wallFilmFaces
+    pitchAxis   (0 1 0);
-    {
+    magUInf     20;
-        type            htcConvection;
+    lRef        1.42;        // Wheelbase length
-        L               1.0;
+    Aref        0.75;        // Estimated
-        value           uniform 0;
-    }
 }
 
 
 // ************************************************************************* //
-/* vim: set filetype=cpp : */

tutorials/incompressible/simpleFoam/motorBike/system/readFields

@@ -5,19 +5,17 @@
 |   \\  /    A nd           | Web:      www.OpenFOAM.com                      |
 |    \\/     M anipulation  |                                                 |
 \*---------------------------------------------------------------------------*/
-FoamFile
+
+// Make sure all fields for functionObjects are loaded. Prevents any
+// problems running with execFlowFunctionObjects.
+readFields
 {
-    version     2.0;
+    // Where to load it from (if not already in solver)
-    format      ascii;
+    functionObjectLibs ("libfieldFunctionObjects.so");
-    class       uniformDimensionedVectorField;
+
-    location    "constant";
+    type            readFields;
-    object      g;
+    fields          (p U k);
 }
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-dimensions      [0 1 -2 0 0 0 0];
-
-value           (0 0 -9.81); 
 
 
 // ************************************************************************* //

tutorials/incompressible/simpleFoam/motorBike/system/streamLines

@@ -0,0 +1,53 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+
+streamLines
+{
+    type            streamLine;
+
+    // Output every
+    outputControl   outputTime;
+    // outputInterval 10;
+
+    setFormat       vtk; //gnuplot; //xmgr; //raw; //jplot;
+
+    // Velocity field to use for tracking.
+    U U;
+
+    // Tracked forwards (+U) or backwards (-U)
+    trackForward    true;
+
+    // Names of fields to sample. Should contain above velocity field!
+    fields (p U k);
+
+    // Steps particles can travel before being removed
+    lifeTime        10000;
+
+    // Number of steps per cell (estimate). Set to 1 to disable subcycling.
+    nSubCycle 5;
+
+    // Cloud name to use
+    cloudName       particleTracks;
+
+    // Seeding method. See the sampleSets in sampleDict.
+    seedSampleSet   uniform;  //cloud;//triSurfaceMeshPointSet;
+
+    uniformCoeffs
+    {
+        type        uniform;
+        axis        x;  //distance;
+
+        // Note: tracks slightly offset so as not to be on a face
+        start       (-1.001 1e-7 0.0011);
+        end         (-1.001 1e-7 1.0011);
+        nPoints     20;
+    }
+}
+
+
+// ************************************************************************* //

tutorials/lagrangian/icoUncoupledKinematicParcelFoam/hopperInitialState/constant/kinematicCloudProperties

@@ -22,13 +22,6 @@ solution
     transient       yes;
     cellValueSourceCorrection off;
 
-    sourceTerms
-    {
-        schemes
-        {
-        }
-    }
-
     interpolationSchemes
     {
         rho             cell;

tutorials/lagrangian/reactingParcelFilmFoam/cylinder/0.org/0/H2O

@@ -22,7 +22,7 @@ internalField   uniform 0.0;
 
 boundaryField
 {
-    "(.*)"
+    "(sides|frontAndBack)"
     {
         type            zeroGradient;
     }

tutorials/lagrangian/reactingParcelFilmFoam/cylinder/0.org/0/N2

@@ -22,7 +22,7 @@ internalField   uniform 0.79;
 
 boundaryField
 {
-    "(.*)"
+    "(sides|frontAndBack)"
     {
         type            zeroGradient;
     }

tutorials/lagrangian/reactingParcelFilmFoam/cylinder/0.org/0/O2

@@ -22,7 +22,7 @@ internalField   uniform 0.21;
 
 boundaryField
 {
-    "(.*)"
+    "(sides|frontAndBack)"
     {
         type            zeroGradient;
     }

tutorials/lagrangian/reactingParcelFilmFoam/cylinder/0.org/0/T

@@ -9,8 +9,8 @@ FoamFile
 {
     version     2.0;
     format      ascii;
-    class       volScalarField;
     location    "0";
+    class       volScalarField;
     object      T;
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -21,16 +21,16 @@ internalField   uniform 300;
 
 boundaryField
 {
-    walls
+    "(sides|frontAndBack)"
     {
         type            fixedValue;
         value           uniform 300;
     }
-    wallFilm
+    region0_to_wallFilmRegion_wallFilmFaces
     {
-        type            zeroGradient;
+        type            fixedValue;
+        value           uniform 300;
     }
 }
 
-
 // ************************************************************************* //

tutorials/lagrangian/reactingParcelFilmFoam/cylinder/0.org/0/U

@@ -21,17 +21,14 @@ internalField   uniform (0 0 0);
 
 boundaryField
 {
-    "(.*)"
+    "(sides|frontAndBack)"
     {
-        type            fixedValue;
+        type            pressureInletOutletVelocity;
         value           uniform (0 0 0);
     }
     region0_to_wallFilmRegion_wallFilmFaces
     {
-        type            directMapped;
+        type            fixedValue;
-        fieldName       Usf;
-        average         ( 0 0 0 );
-        setAverage      no;
         value           uniform (0 0 0);
     }
 }

tutorials/lagrangian/reactingParcelFilmFoam/cylinder/0.org/0/p

@@ -0,0 +1,36 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    location    "0";
+    object      p;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [1 -1 -2 0 0 0 0];
+
+internalField   uniform 100000;
+
+boundaryField
+{
+    "(sides|frontAndBack)"
+    {
+        type            calculated;
+        value           $internalField;
+    }
+    region0_to_wallFilmRegion_wallFilmFaces
+    {
+        type            calculated;
+        value           $internalField;
+    }
+}
+
+// ************************************************************************* //

tutorials/lagrangian/reactingParcelFilmFoam/cylinder/0.org/0/p_rgh

@@ -11,7 +11,7 @@ FoamFile
     format      ascii;
     class       volScalarField;
     location    "0";
-    object      p;
+    object      p_rgh;
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -21,13 +21,13 @@ internalField   uniform 100000;
 
 boundaryField
 {
-    "(.*)"
+    "(sides|frontAndBack)"
     {
-        type            zeroGradient; // buoyantPressure;
+        type            buoyantPressure;
     }
     region0_to_wallFilmRegion_wallFilmFaces
     {
-        type            zeroGradient; // buoyantPressure;
+        type            buoyantPressure;
     }
 }
 

tutorials/lagrangian/reactingParcelFilmFoam/cylinder/0.org/0/wallFilmRegion/Tf

@@ -0,0 +1,41 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    location    "0/wallFilmRegion";
+    object      Tf;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 0 0 1 0 0 0];
+
+internalField   uniform  300;
+
+boundaryField
+{
+    wallFilmFaces_top
+    {
+        type            zeroGradient;
+    }
+
+    region0_to_wallFilmRegion_wallFilmFaces
+    {
+        type            zeroGradient;
+    }
+
+    filmWalls
+    {
+        type            zeroGradient;
+    }
+}
+
+
+// ************************************************************************* //

tutorials/lagrangian/reactingParcelFilmFoam/cylinder/0.org/0/wallFilmRegion/Uf

@@ -10,8 +10,8 @@ FoamFile
     version     2.0;
     format      ascii;
     class       volVectorField;
-    location    "0";
+    location    "0/wallFilmRegion";
-    object      U;
+    object      Uf;
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -21,12 +21,18 @@ internalField   uniform (0 0 0);
 
 boundaryField
 {
-    walls
+    wallFilmFaces_top
+    {
+        type            slip;
+    }
+
+    region0_to_wallFilmRegion_wallFilmFaces
     {
         type            fixedValue;
         value           uniform (0 0 0);
     }
-    wallFilm
+
+    filmWalls
     {
         type            fixedValue;
         value           uniform (0 0 0);

tutorials/lagrangian/reactingParcelFilmFoam/cylinder/0.org/0/wallFilmRegion/deltaf

@@ -0,0 +1,41 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    location    "0/wallFilmRegion";
+    object      deltaf;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 1 0 0 0 0 0];
+
+internalField   uniform 0;
+
+boundaryField
+{
+    wallFilmFaces_top
+    {
+        type            zeroGradient;
+    }
+
+    region0_to_wallFilmRegion_wallFilmFaces
+    {
+        type            zeroGradient;
+    }
+
+    filmWalls
+    {
+        type            zeroGradient;
+    }
+}
+
+
+// ************************************************************************* //

tutorials/lagrangian/reactingParcelFilmFoam/cylinder/0.org/H2O

@@ -22,7 +22,7 @@ internalField   uniform 0.0;
 
 boundaryField
 {
-    "(.*)"
+    "(sides|frontAndBack)"
     {
         type            zeroGradient;
     }

tutorials/lagrangian/reactingParcelFilmFoam/cylinder/0.org/N2

@@ -22,7 +22,7 @@ internalField   uniform 0.79;
 
 boundaryField
 {
-    "(.*)"
+    "(sides|frontAndBack)"
     {
         type            zeroGradient;
     }

tutorials/lagrangian/reactingParcelFilmFoam/cylinder/0.org/O2

@@ -22,7 +22,7 @@ internalField   uniform 0.21;
 
 boundaryField
 {
-    "(.*)"
+    "(sides|frontAndBack)"
     {
         type            zeroGradient;
     }

tutorials/lagrangian/reactingParcelFilmFoam/cylinder/0.org/T

@@ -9,8 +9,8 @@ FoamFile
 {
     version     2.0;
     format      ascii;
-    class       volScalarField;
     location    "0";
+    class       volScalarField;
     object      T;
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -21,16 +21,16 @@ internalField   uniform 300;
 
 boundaryField
 {
-    walls
+    "(sides|frontAndBack)"
     {
         type            fixedValue;
         value           uniform 300;
     }
-    wallFilm
+    region0_to_wallFilmRegion_wallFilmFaces
     {
-        type            zeroGradient;
+        type            fixedValue;
+        value           uniform 300;
     }
 }
 
-
 // ************************************************************************* //

tutorials/lagrangian/reactingParcelFilmFoam/cylinder/0.org/U

@@ -21,17 +21,14 @@ internalField   uniform (0 0 0);
 
 boundaryField
 {
-    "(.*)"
+    "(sides|frontAndBack)"
     {
-        type            fixedValue;
+        type            pressureInletOutletVelocity;
         value           uniform (0 0 0);
     }
     region0_to_wallFilmRegion_wallFilmFaces
     {
-        type            directMapped;
+        type            fixedValue;
-        fieldName       Usf;
-        average         ( 0 0 0 );
-        setAverage      no;
         value           uniform (0 0 0);
     }
 }

tutorials/lagrangian/reactingParcelFilmFoam/cylinder/0.org/p

@@ -0,0 +1,36 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    location    "0";
+    object      p;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [1 -1 -2 0 0 0 0];
+
+internalField   uniform 100000;
+
+boundaryField
+{
+    "(sides|frontAndBack)"
+    {
+        type            calculated;
+        value           $internalField;
+    }
+    region0_to_wallFilmRegion_wallFilmFaces
+    {
+        type            calculated;
+        value           $internalField;
+    }
+}
+
+// ************************************************************************* //

tutorials/lagrangian/reactingParcelFilmFoam/cylinder/0.org/p_rgh

@@ -11,7 +11,7 @@ FoamFile
     format      ascii;
     class       volScalarField;
     location    "0";
-    object      p;
+    object      p_rgh;
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -21,13 +21,13 @@ internalField   uniform 100000;
 
 boundaryField
 {
-    "(.*)"
+    "(sides|frontAndBack)"
     {
-        type            zeroGradient; // buoyantPressure;
+        type            buoyantPressure;
     }
     region0_to_wallFilmRegion_wallFilmFaces
     {
-        type            zeroGradient; // buoyantPressure;
+        type            buoyantPressure;
     }
 }
 

tutorials/lagrangian/reactingParcelFilmFoam/cylinder/0.org/wallFilmRegion/Tf

@@ -17,7 +17,7 @@ FoamFile
 
 dimensions      [0 0 0 1 0 0 0];
 
-internalField   uniform  288;
+internalField   uniform  300;
 
 boundaryField
 {
@@ -27,17 +27,11 @@ boundaryField
     }
 
     region0_to_wallFilmRegion_wallFilmFaces
-    {
-        type            fixedValue;
-        value           uniform 312.3;
-    }
-
-    left
     {
         type            zeroGradient;
     }
 
-    right
+    sides
     {
         type            zeroGradient;
     }
@@ -48,4 +42,5 @@ boundaryField
     }
 }
 
+
 // ************************************************************************* //

tutorials/lagrangian/reactingParcelFilmFoam/cylinder/0.org/wallFilmRegion/Uf

@@ -32,13 +32,7 @@ boundaryField
         value           uniform (0 0 0);
     }
 
-    left
+    sides
-    {
-        type            fixedValue;
-        value           uniform (0 0 0);
-    }
-
-    right
     {
         type            fixedValue;
         value           uniform (0 0 0);
@@ -46,9 +40,9 @@ boundaryField
 
     frontAndBack
     {
-        type            fixedValue;
+        type            slip;
-        value           uniform (0 0 0);
     }
 }
 
+
 // ************************************************************************* //

tutorials/lagrangian/reactingParcelFilmFoam/cylinder/0.org/wallFilmRegion/deltaf

@@ -0,0 +1,46 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    location    "0/wallFilmRegion";
+    object      deltaf;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 1 0 0 0 0 0];
+
+internalField   uniform 0;
+
+boundaryField
+{
+    wallFilmFaces_top
+    {
+        type            zeroGradient;
+    }
+
+    region0_to_wallFilmRegion_wallFilmFaces
+    {
+        type            zeroGradient;
+    }
+
+    sides
+    {
+        type            zeroGradient;
+    }
+
+    frontAndBack
+    {
+        type            zeroGradient;
+    }
+}
+
+
+// ************************************************************************* //

tutorials/lagrangian/reactingParcelFilmFoam/cylinder/0/0/H2O

@@ -0,0 +1,36 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    location    "0";
+    object      H2O;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+
+dimensions      [0 0 0 0 0 0 0];
+
+internalField   uniform 0.0;
+
+boundaryField
+{
+    "(sides|frontAndBack)"
+    {
+        type            zeroGradient;
+    }
+    region0_to_wallFilmRegion_wallFilmFaces
+    {
+        type            zeroGradient;
+    }
+}
+
+
+// ************************************************************************* //

tutorials/lagrangian/reactingParcelFilmFoam/cylinder/0/0/N2

@@ -0,0 +1,36 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    location    "0";
+    object      N2;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+
+dimensions      [0 0 0 0 0 0 0];
+
+internalField   uniform 0.79;
+
+boundaryField
+{
+    "(sides|frontAndBack)"
+    {
+        type            zeroGradient;
+    }
+    region0_to_wallFilmRegion_wallFilmFaces
+    {
+        type            zeroGradient;
+    }
+}
+
+
+// ************************************************************************* //

tutorials/lagrangian/reactingParcelFilmFoam/cylinder/0/0/O2

@@ -0,0 +1,36 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    location    "0";
+    object      O2;
+}
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 0 0 0 0 0 0];
+
+internalField   uniform 0.21;
+
+boundaryField
+{
+    "(sides|frontAndBack)"
+    {
+        type            zeroGradient;
+    }
+    region0_to_wallFilmRegion_wallFilmFaces
+    {
+        type            zeroGradient;
+    }
+}
+
+
+// ************************************************************************* //

tutorials/lagrangian/reactingParcelFilmFoam/cylinder/0/0/T

@@ -0,0 +1,36 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    location    "0";
+    class       volScalarField;
+    object      T;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 0 0 1 0 0 0];
+
+internalField   uniform 300;
+
+boundaryField
+{
+    "(sides|frontAndBack)"
+    {
+        type            fixedValue;
+        value           uniform 300;
+    }
+    region0_to_wallFilmRegion_wallFilmFaces
+    {
+        type            fixedValue;
+        value           uniform 300;
+    }
+}
+
+// ************************************************************************* //

tutorials/lagrangian/reactingParcelFilmFoam/cylinder/0/0/U

@@ -0,0 +1,37 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volVectorField;
+    location    "0";
+    object      U;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 1 -1 0 0 0 0];
+
+internalField   uniform (0 0 0); 
+
+boundaryField
+{
+    "(sides|frontAndBack)"
+    {
+        type            pressureInletOutletVelocity;
+        value           uniform (0 0 0);
+    }
+    region0_to_wallFilmRegion_wallFilmFaces
+    {
+        type            fixedValue;
+        value           uniform (0 0 0);
+    }
+}
+
+
+// ************************************************************************* //

tutorials/lagrangian/reactingParcelFilmFoam/cylinder/0/0/p

@@ -0,0 +1,36 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    location    "0";
+    object      p;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [1 -1 -2 0 0 0 0];
+
+internalField   uniform 100000;
+
+boundaryField
+{
+    "(sides|frontAndBack)"
+    {
+        type            calculated;
+        value           $internalField;
+    }
+    region0_to_wallFilmRegion_wallFilmFaces
+    {
+        type            calculated;
+        value           $internalField;
+    }
+}
+
+// ************************************************************************* //

tutorials/lagrangian/reactingParcelFilmFoam/cylinder/0/0/p_rgh

@@ -0,0 +1,34 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    location    "0";
+    object      p_rgh;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [1 -1 -2 0 0 0 0];
+
+internalField   uniform 100000;
+
+boundaryField
+{
+    "(sides|frontAndBack)"
+    {
+        type            buoyantPressure;
+    }
+    region0_to_wallFilmRegion_wallFilmFaces
+    {
+        type            buoyantPressure;
+    }
+}
+
+// ************************************************************************* //

tutorials/lagrangian/reactingParcelFilmFoam/cylinder/0/0/wallFilmRegion/Tf

@@ -0,0 +1,41 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    location    "0/wallFilmRegion";
+    object      Tf;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 0 0 1 0 0 0];
+
+internalField   uniform  300;
+
+boundaryField
+{
+    wallFilmFaces_top
+    {
+        type            zeroGradient;
+    }
+
+    region0_to_wallFilmRegion_wallFilmFaces
+    {
+        type            zeroGradient;
+    }
+
+    filmWalls
+    {
+        type            zeroGradient;
+    }
+}
+
+
+// ************************************************************************* //

tutorials/lagrangian/reactingParcelFilmFoam/cylinder/0/0/wallFilmRegion/Uf

@@ -10,8 +10,8 @@ FoamFile
     version     2.0;
     format      ascii;
     class       volVectorField;
-    location    "0";
+    location    "0/wallFilmRegion";
-    object      U;
+    object      Uf;
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -21,12 +21,18 @@ internalField   uniform (0 0 0);
 
 boundaryField
 {
-    walls
+    wallFilmFaces_top
+    {
+        type            slip;
+    }
+
+    region0_to_wallFilmRegion_wallFilmFaces
     {
         type            fixedValue;
         value           uniform (0 0 0);
     }
-    wallFilm
+
+    filmWalls
     {
         type            fixedValue;
         value           uniform (0 0 0);

tutorials/lagrangian/reactingParcelFilmFoam/cylinder/0/0/wallFilmRegion/deltaf

@@ -0,0 +1,41 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    location    "0/wallFilmRegion";
+    object      deltaf;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 1 0 0 0 0 0];
+
+internalField   uniform 0;
+
+boundaryField
+{
+    wallFilmFaces_top
+    {
+        type            zeroGradient;
+    }
+
+    region0_to_wallFilmRegion_wallFilmFaces
+    {
+        type            zeroGradient;
+    }
+
+    filmWalls
+    {
+        type            zeroGradient;
+    }
+}
+
+
+// ************************************************************************* //