ENH: Updated surface film functionality

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();
@@ -475,15 +453,11 @@ kinematicSingleLayer::kinematicSingleLayer
     momentumPredictor_(solution().subDict("PISO").lookup("momentumPredictor")),
     nOuterCorr_(readLabel(solution().subDict("PISO").lookup("nOuterCorr"))),
     nCorr_(readLabel(solution().subDict("PISO").lookup("nCorr"))),
-    nNonOrthCorr_
+    nNonOrthCorr_(readLabel(solution().subDict("PISO").lookup("nNonOrthCorr"))),
-    (
-        readLabel(solution().subDict("PISO").lookup("nNonOrthCorr"))
-    ),
 
     cumulativeContErr_(0.0),
 
     Cf_(readScalar(coeffs().lookup("Cf"))),
-    deltaStable_(coeffs().lookup("deltaStable")),
 
     rho_
     (
@@ -607,11 +581,11 @@ kinematicSingleLayer::kinematicSingleLayer
         dimLength*dimMass/dimTime
     ),
 
-    massForPrimary_
+    primaryMassTrans_
     (
         IOobject
         (
-            "massForPrimary",
+            "primaryMassTrans",
             time().timeName(),
             regionMesh(),
             IOobject::NO_READ,
@@ -621,11 +595,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 +623,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 +767,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 +811,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 +821,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 +932,12 @@ const volVectorField& kinematicSingleLayer::Uw() const
 }
 
 
+const surfaceScalarField& kinematicSingleLayer::phi() const
+{
+    return phi_;
+}
+
+
 const volScalarField& kinematicSingleLayer::rho() const
 {
     return rho_;
@@ -1002,21 +999,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 +1041,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 +1069,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 +1083,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,352 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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,58 @@ 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.C

@@ -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) 2009-2010 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License

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

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2010-2011 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2010-2010 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -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

@@ -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) 2009-2010 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -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
                 (