From e30cbfc6db320e2b34e4750545d1a3ed24627535 Mon Sep 17 00:00:00 2001
From: sergio <s.ferraris@opencfd.co.uk>
Date: Tue, 13 Jul 2021 13:12:38 +0100
Subject: [PATCH] ENH: liquidFilm subModels: add injection models

ENH: Several additions to thinFilm models

1) Changing the faSolution entry nonOrthoCorr to nFilmCorr entry.
   This entry specify the number of iterations on the film thickness
   Eq which does not correct for laplacian non-orthogonality.

2) Changing laminar turbulent model for the skin friction exerted by
   the flow to the film.

3) Adding optional relaxation factors for source terms coming from
   the added mass of the lagragian particles into the film such as
   rhoSp, pnSp and Usp. Similarly the film pressure (pf) can  be
   relaxed in the inner loops.

4) Adding optional entries to curvatureSepration injection model.
   minimum (fThreshold) force and minimum curvature (minInvR1) for
   separation were added to have more control on determining the
   points of film separation

ENH: ThermoSurfaceFilm in lagragian to search for finiteArea film models
---
 src/Allwmake                                  |   6 +-
 .../InjectionModel/InjectionModel.C           |   3 +-
 .../SurfaceFilmModel/SurfaceFilmModel.C       |  87 ++++-
 .../SurfaceFilmModel/SurfaceFilmModel.H       |  26 +-
 .../SurfaceFilmModel/SurfaceFilmModelI.H      |  13 +
 .../ThermoSurfaceFilm/ThermoSurfaceFilm.C     | 252 +++++++-------
 .../ThermoSurfaceFilm/ThermoSurfaceFilm.H     |  34 +-
 src/regionFaModels/Make/files                 |  16 +
 .../liquidFilm/liquidFilmBase.C               | 179 +++++-----
 .../liquidFilm/liquidFilmBase.H               |  70 +++-
 .../liquidFilmModel/liquidFilmModel.C         | 143 +++++++-
 .../liquidFilmModel/liquidFilmModel.H         |  55 ++-
 .../liquidFilm/subModels/filmSubModelBase.C   | 110 ++++++
 .../liquidFilm/subModels/filmSubModelBase.H   | 139 ++++++++
 .../liquidFilm/subModels/filmSubModelBaseI.H  |  57 +++
 .../subModels/filmSubModelBaseTemplates.C     |  58 ++++
 .../filmTurbulenceModel/filmTurbulenceModel.C |   4 +-
 .../filmTurbulenceModel/filmTurbulenceModel.H |   4 +-
 .../filmTurbulenceModel/laminar/laminar.C     |   4 +-
 .../filmTurbulenceModel/laminar/laminar.H     |   6 +-
 .../curvatureSeparation/curvatureSeparation.C | 324 ++++++++++++++++++
 .../curvatureSeparation/curvatureSeparation.H | 167 +++++++++
 .../injectionModel/injectionModel.C           | 106 ++++++
 .../injectionModel/injectionModel.H           | 169 +++++++++
 .../injectionModel/injectionModelNew.C        |  73 ++++
 .../injectionModelList/injectionModelList.C   | 169 +++++++++
 .../injectionModelList/injectionModelList.H   | 127 +++++++
 27 files changed, 2117 insertions(+), 284 deletions(-)
 create mode 100644 src/regionFaModels/liquidFilm/subModels/filmSubModelBase.C
 create mode 100644 src/regionFaModels/liquidFilm/subModels/filmSubModelBase.H
 create mode 100644 src/regionFaModels/liquidFilm/subModels/filmSubModelBaseI.H
 create mode 100644 src/regionFaModels/liquidFilm/subModels/filmSubModelBaseTemplates.C
 create mode 100644 src/regionFaModels/liquidFilm/subModels/kinematic/injectionModel/curvatureSeparation/curvatureSeparation.C
 create mode 100644 src/regionFaModels/liquidFilm/subModels/kinematic/injectionModel/curvatureSeparation/curvatureSeparation.H
 create mode 100644 src/regionFaModels/liquidFilm/subModels/kinematic/injectionModel/injectionModel/injectionModel.C
 create mode 100644 src/regionFaModels/liquidFilm/subModels/kinematic/injectionModel/injectionModel/injectionModel.H
 create mode 100644 src/regionFaModels/liquidFilm/subModels/kinematic/injectionModel/injectionModel/injectionModelNew.C
 create mode 100644 src/regionFaModels/liquidFilm/subModels/kinematic/injectionModel/injectionModelList/injectionModelList.C
 create mode 100644 src/regionFaModels/liquidFilm/subModels/kinematic/injectionModel/injectionModelList/injectionModelList.H

diff --git a/src/Allwmake b/src/Allwmake
index e488f93978..0af702324e 100755
--- a/src/Allwmake
+++ b/src/Allwmake
@@ -83,14 +83,18 @@ wmake $targetType ODE
 thermophysicalModels/Allwmake $targetType $*
 TurbulenceModels/Allwmake $targetType $*
 wmake $targetType combustionModels
+
+wmakeLnInclude -u regionFaModels
+wmakeLnInclude -u faOptions
 regionModels/Allwmake $targetType $*
-lagrangian/Allwmake $targetType $*
 
 wmake $targetType fvOptions
 wmake $targetType faOptions
 wmake $targetType fvMotionSolver
 wmake $targetType regionFaModels
 
+lagrangian/Allwmake $targetType $*
+
 wmake $targetType overset
 
 # snappyHexMesh uses overset voxelMesh
diff --git a/src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/InjectionModel/InjectionModel.C b/src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/InjectionModel/InjectionModel.C
index 923af2c002..94e4961ca8 100644
--- a/src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/InjectionModel/InjectionModel.C
+++ b/src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/InjectionModel/InjectionModel.C
@@ -5,8 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2011-2017 OpenFOAM Foundation
-    Copyright (C) 2016-2020 OpenCFD Ltd.
+    Copyright (C) 2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
diff --git a/src/lagrangian/intermediate/submodels/Kinematic/SurfaceFilmModel/SurfaceFilmModel/SurfaceFilmModel.C b/src/lagrangian/intermediate/submodels/Kinematic/SurfaceFilmModel/SurfaceFilmModel/SurfaceFilmModel.C
index 4696ac7021..d9e474202b 100644
--- a/src/lagrangian/intermediate/submodels/Kinematic/SurfaceFilmModel/SurfaceFilmModel/SurfaceFilmModel.C
+++ b/src/lagrangian/intermediate/submodels/Kinematic/SurfaceFilmModel/SurfaceFilmModel/SurfaceFilmModel.C
@@ -41,13 +41,16 @@ Foam::SurfaceFilmModel<CloudType>::SurfaceFilmModel(CloudType& owner)
     CloudSubModelBase<CloudType>(owner),
     g_(owner.g()),
     ejectedParcelType_(0),
+    injectionOffset_(1.1),
+    minDiameter_(0),
     massParcelPatch_(0),
     diameterParcelPatch_(0),
     UFilmPatch_(0),
     rhoFilmPatch_(0),
     deltaFilmPatch_(0),
     nParcelsTransferred_(0),
-    nParcelsInjected_(0)
+    nParcelsInjected_(0),
+    totalMassTransferred_(0)
 {}
 
 
@@ -63,7 +66,15 @@ Foam::SurfaceFilmModel<CloudType>::SurfaceFilmModel
     g_(owner.g()),
     ejectedParcelType_
     (
-        this->coeffDict().getOrDefault("ejectedParcelType", -1)
+        this->coeffDict().template getOrDefault<label>("ejectedParcelType", -1)
+    ),
+    injectionOffset_
+    (
+        this->coeffDict().template getOrDefault<scalar>("injectionOffset", 1.1)
+    ),
+    minDiameter_
+    (
+        this->coeffDict().template getOrDefault<scalar>("minDiameter", -1)
     ),
     massParcelPatch_(0),
     diameterParcelPatch_(0),
@@ -71,7 +82,8 @@ Foam::SurfaceFilmModel<CloudType>::SurfaceFilmModel
     rhoFilmPatch_(0),
     deltaFilmPatch_(owner.mesh().boundary().size()),
     nParcelsTransferred_(0),
-    nParcelsInjected_(0)
+    nParcelsInjected_(0),
+    totalMassTransferred_()
 {}
 
 
@@ -84,13 +96,16 @@ Foam::SurfaceFilmModel<CloudType>::SurfaceFilmModel
     CloudSubModelBase<CloudType>(sfm),
     g_(sfm.g_),
     ejectedParcelType_(sfm.ejectedParcelType_),
+    injectionOffset_(sfm.injectionOffset_),
+    minDiameter_(sfm.minDiameter_),
     massParcelPatch_(sfm.massParcelPatch_),
     diameterParcelPatch_(sfm.diameterParcelPatch_),
     UFilmPatch_(sfm.UFilmPatch_),
     rhoFilmPatch_(sfm.rhoFilmPatch_),
     deltaFilmPatch_(sfm.deltaFilmPatch_),
     nParcelsTransferred_(sfm.nParcelsTransferred_),
-    nParcelsInjected_(sfm.nParcelsInjected_)
+    nParcelsInjected_(sfm.nParcelsInjected_),
+    totalMassTransferred_(sfm.totalMassTransferred_)
 {}
 
 
@@ -130,7 +145,7 @@ void Foam::SurfaceFilmModel<CloudType>::injectParticles
                     diameterParcelPatch_[j],
                     deltaFilmPatch_[primaryPatchi][j]
                 );
-            const point pos = Cf[j] - 1.1*offset*Sf[j]/magSf[j];
+            const point pos = Cf[j] - injectionOffset_*offset*Sf[j]/magSf[j];
 
             // Create a new parcel
             parcelType* pPtr =
@@ -144,7 +159,6 @@ void Foam::SurfaceFilmModel<CloudType>::injectParticles
             if (pPtr->nParticle() > 0.001)
             {
                 // Check new parcel properties
-//                cloud.checkParcelProperties(*pPtr, 0.0, true);
                 cloud.checkParcelProperties(*pPtr, 0.0, false);
 
                 // Add the new parcel to the cloud
@@ -215,7 +229,7 @@ void Foam::SurfaceFilmModel<CloudType>::inject(TrackCloudType& cloud)
             >(names[i]);
 
         // Check that it is a type areaFilm
-        if (regionFa && isA<areaFilm>(*regionFa))
+        if (regionFa && isA<areaFilm>(*regionFa) && regionFa->active())
         {
             areaFilm& film =
                 const_cast<areaFilm&>(refCast<const areaFilm>(*regionFa));
@@ -227,6 +241,22 @@ void Foam::SurfaceFilmModel<CloudType>::inject(TrackCloudType& cloud)
             cacheFilmFields(patchId, film);
 
             injectParticles(patchId, injectorCellsPatch, cloud);
+
+            forAll(injectorCellsPatch, facei)
+            {
+                if (diameterParcelPatch_[facei] > 0)
+                {
+                    film.addSources
+                    (
+                        patchId,
+                        facei,
+                       -massParcelPatch_[facei],    // mass
+                        Zero,                       // tangential momentum
+                        Zero,                       // impingement
+                        Zero                        // energy
+                    );
+                }
+            }
         }
     }
 }
@@ -266,22 +296,26 @@ void Foam::SurfaceFilmModel<CloudType>::cacheFilmFields
     const regionModels::areaSurfaceFilmModels::liquidFilmBase& filmModel
 )
 {
-    // WIP: the finite area film does not support splashing so far
-    //massParcelPatch_ = filmModel.cloudMassTrans().boundaryField()[filmPatchi];
-    //filmModel.toPrimary(filmPatchi, massParcelPatch_);
-    massParcelPatch_.setSize(filmModel.Uf().size(), Zero);
-    diameterParcelPatch_.setSize(filmModel.Uf().size(), Zero);
-    //    filmModel.cloudDiameterTrans().boundaryField()[filmPatchi];
-    //filmModel.toPrimary(filmPatchi, diameterParcelPatch_, maxEqOp<scalar>());
     const volSurfaceMapping& map = filmModel.region().vsm();
 
-    UFilmPatch_.setSize(filmModel.Uf().size(), Zero);
+    massParcelPatch_.setSize(filmModel.Uf().size(), Zero);
+
+    const scalarField& massParcelPatch =
+        filmModel.cloudMassTrans().boundaryField()[filmPatchi];
+
+    map.mapToField(massParcelPatch, massParcelPatch_);
+
+    const scalarField&  diameterParcelPatch =
+        filmModel.cloudDiameterTrans().boundaryField()[filmPatchi];
+
+    diameterParcelPatch_.setSize(filmModel.Uf().size(), Zero);
+    map.mapToField(diameterParcelPatch, diameterParcelPatch_);
+
+    UFilmPatch_.setSize(filmModel.Uf().size(), vector::zero);
     map.mapToField(filmModel.Uf(), UFilmPatch_);
-    //filmModel.toPrimary(filmPatchi, UFilmPatch_);
 
     rhoFilmPatch_.setSize(UFilmPatch_.size(), Zero);
     map.mapToField(filmModel.rho(), rhoFilmPatch_);
-    //filmModel.toPrimary(filmPatchi, rhoFilmPatch_);
 
     deltaFilmPatch_[filmPatchi].setSize(UFilmPatch_.size(), Zero);
     map.mapToField(filmModel.h(), deltaFilmPatch_[filmPatchi]);
@@ -303,6 +337,14 @@ void Foam::SurfaceFilmModel<CloudType>::setParcelProperties
 
     p.nParticle() = massParcelPatch_[filmFacei]/p.rho()/vol;
 
+    if (minDiameter_ != -1)
+    {
+        if (p.d() < minDiameter_)
+        {
+            p.nParticle() = 0;
+        }
+    }
+
     if (ejectedParcelType_ >= 0)
     {
         p.typeId() = ejectedParcelType_;
@@ -319,22 +361,33 @@ void Foam::SurfaceFilmModel<CloudType>::info(Ostream& os)
     label nInject0 =
         this->template getModelProperty<label>("nParcelsInjected");
 
+    scalar massTransferred0 =
+        this->template getModelProperty<scalar>("massTransferred");
+
     label nTransTotal =
         nTrans0 + returnReduce(nParcelsTransferred_, sumOp<label>());
 
     label nInjectTotal =
         nInject0 + returnReduce(nParcelsInjected_, sumOp<label>());
 
+    scalar massTransferredTotal =
+        massTransferred0 + returnReduce(totalMassTransferred_, sumOp<scalar>());
+
+
     os  << "    Surface film:" << nl
         << "      - parcels absorbed            = " << nTransTotal << nl
+        << "      - mass absorbed               = " << massTransferredTotal << nl
         << "      - parcels ejected             = " << nInjectTotal << endl;
 
     if (this->writeTime())
     {
         this->setModelProperty("nParcelsTransferred", nTransTotal);
         this->setModelProperty("nParcelsInjected", nInjectTotal);
+        this->setModelProperty("massTransferred", massTransferredTotal);
+
         nParcelsTransferred_ = 0;
         nParcelsInjected_ = 0;
+        totalMassTransferred_ = 0;
     }
 }
 
diff --git a/src/lagrangian/intermediate/submodels/Kinematic/SurfaceFilmModel/SurfaceFilmModel/SurfaceFilmModel.H b/src/lagrangian/intermediate/submodels/Kinematic/SurfaceFilmModel/SurfaceFilmModel/SurfaceFilmModel.H
index 0e5a7a4cef..410dd63e23 100644
--- a/src/lagrangian/intermediate/submodels/Kinematic/SurfaceFilmModel/SurfaceFilmModel/SurfaceFilmModel.H
+++ b/src/lagrangian/intermediate/submodels/Kinematic/SurfaceFilmModel/SurfaceFilmModel/SurfaceFilmModel.H
@@ -6,6 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2017 OpenFOAM Foundation
+    Copyright (C) 2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -46,6 +47,7 @@ SourceFiles
 #include "runTimeSelectionTables.H"
 #include "CloudSubModelBase.H"
 
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 namespace Foam
@@ -96,14 +98,20 @@ protected:
         //  type
         label ejectedParcelType_;
 
+        //- Injection offset position
+        scalar injectionOffset_;
+
+        //- Minimum diameter particle injection
+        scalar minDiameter_;
+
 
         // Cached injector fields per film patch
 
             //- Parcel mass / patch face
-            scalarList massParcelPatch_;
+            scalarField massParcelPatch_;
 
             //- Parcel diameter / patch face
-            scalarList diameterParcelPatch_;
+            scalarField diameterParcelPatch_;
 
             //- Film velocity / patch face
             Field<vector> UFilmPatch_;
@@ -124,6 +132,12 @@ protected:
             label nParcelsInjected_;
 
 
+        // Total mass info
+
+            //- Total mass transferred to the film
+            scalar totalMassTransferred_;
+
+
     // Protected functions
 
         //- Cache the film fields in preparation for injection
@@ -138,7 +152,7 @@ protected:
         virtual void cacheFilmFields
         (
             const label primaryPatchi,
-            const regionModels::areaSurfaceFilmModels::liquidFilmBase&
+            const areaFilm&
         );
 
         //- Inject particles in cloud
@@ -232,6 +246,12 @@ public:
             //  the film model
             inline label& nParcelsInjected();
 
+            //- Return non-const total mass transferred
+            inline scalar& totalMassTransferred();
+
+            //- Return consr access to mass transferred
+            inline scalar totalMassTransferred() const;
+
 
         // Member Functions
 
diff --git a/src/lagrangian/intermediate/submodels/Kinematic/SurfaceFilmModel/SurfaceFilmModel/SurfaceFilmModelI.H b/src/lagrangian/intermediate/submodels/Kinematic/SurfaceFilmModel/SurfaceFilmModel/SurfaceFilmModelI.H
index 67420350c6..4240d00ea3 100644
--- a/src/lagrangian/intermediate/submodels/Kinematic/SurfaceFilmModel/SurfaceFilmModel/SurfaceFilmModelI.H
+++ b/src/lagrangian/intermediate/submodels/Kinematic/SurfaceFilmModel/SurfaceFilmModel/SurfaceFilmModelI.H
@@ -63,5 +63,18 @@ Foam::label Foam::SurfaceFilmModel<CloudType>::nParcelsInjected() const
     return nParcelsInjected_;
 }
 
+template<class CloudType>
+Foam::scalar& Foam::SurfaceFilmModel<CloudType>::totalMassTransferred()
+{
+    return totalMassTransferred_;
+}
+
+
+template<class CloudType>
+Foam::scalar Foam::SurfaceFilmModel<CloudType>::totalMassTransferred() const
+{
+    return totalMassTransferred_;
+}
+
 
 // ************************************************************************* //
diff --git a/src/lagrangian/intermediate/submodels/Thermodynamic/SurfaceFilmModel/ThermoSurfaceFilm/ThermoSurfaceFilm.C b/src/lagrangian/intermediate/submodels/Thermodynamic/SurfaceFilmModel/ThermoSurfaceFilm/ThermoSurfaceFilm.C
index 70594de693..bfe1f3fa9b 100644
--- a/src/lagrangian/intermediate/submodels/Thermodynamic/SurfaceFilmModel/ThermoSurfaceFilm/ThermoSurfaceFilm.C
+++ b/src/lagrangian/intermediate/submodels/Thermodynamic/SurfaceFilmModel/ThermoSurfaceFilm/ThermoSurfaceFilm.C
@@ -131,57 +131,51 @@ Foam::vector Foam::ThermoSurfaceFilm<CloudType>::splashDirection
 
 
 template<class CloudType>
-void Foam::ThermoSurfaceFilm<CloudType>::absorbInteraction
-(
-    const parcelType& p,
-    const polyPatch& pp,
-    const label facei,
-    const scalar mass,
-    bool& keepParticle
-)
+void Foam::ThermoSurfaceFilm<CloudType>::init()
 {
-    if (debug)
-    {
-        Info<< "Parcel " << p.origId() << " absorbInteraction" << endl;
-    }
+    const fvMesh& mesh = this->owner().mesh();
 
-    // Patch face normal
-    const vector& nf = pp.faceNormals()[facei];
-
-    // Patch velocity
-    const vector& Up = this->owner().U().boundaryField()[pp.index()][facei];
-
-    // Relative parcel velocity
-    const vector Urel = p.U() - Up;
-
-    // Parcel normal velocity
-    const vector Un = nf*(Urel & nf);
-
-    // Parcel tangential velocity
-    const vector Ut = Urel - Un;
-
-    if (filmModel_)
-    {
-        filmModel_->addSources
+    // set up filmModel pointer
+    filmModel_ =
+        const_cast<regionFilm*>
         (
-            pp.index(),
-            facei,
-            mass,                           // mass
-            mass*Ut,                        // tangential momentum
-            mass*mag(Un),                   // impingement pressure
-            mass*p.hs()                     // energy
+            mesh.time().objectRegistry::template findObject
+            <
+                regionFilm
+            >
+            (
+                "surfaceFilmProperties"
+            )
         );
 
-        this->nParcelsTransferred()++;
+    // set up areaFilms
+    const wordList names =
+        mesh.time().objectRegistry::template
+            sortedNames<regionModels::regionFaModel>();
 
-        keepParticle = false;
+    forAll(names, i)
+    {
+        const regionModels::regionFaModel* regionFa =
+            mesh.time().objectRegistry::template findObject
+            <
+                regionModels::regionFaModel
+            >(names[i]);
+
+        if (regionFa && isA<areaFilm>(*regionFa))
+        {
+            areaFilm& film =
+                const_cast<areaFilm&>(refCast<const areaFilm>(*regionFa));
+            areaFilms_.append(&film);
+        }
     }
 }
 
+
 template<class CloudType>
+template<class filmType>
 void Foam::ThermoSurfaceFilm<CloudType>::absorbInteraction
 (
-    areaFilm& filmModel,
+    filmType& film,
     const parcelType& p,
     const polyPatch& pp,
     const label facei,
@@ -189,7 +183,6 @@ void Foam::ThermoSurfaceFilm<CloudType>::absorbInteraction
     bool& keepParticle
 )
 {
-
     if (debug)
     {
         Info<< "Parcel " << p.origId() << " absorbInteraction" << endl;
@@ -209,21 +202,22 @@ void Foam::ThermoSurfaceFilm<CloudType>::absorbInteraction
 
     // Parcel tangential velocity
     const vector Ut = Urel - Un;
-DebugVar(mass)
-    filmModel.addSources
+
+    film.addSources
     (
         pp.index(),
         facei,
-        mass,                           // mass
-        mass*Ut,                        // tangential momentum
-        mass*mag(Un),                   // impingement pressure
-        mass*p.hs()                     // energy
+        mass,                             // mass
+        mass*Ut,                          // tangential momentum
+        mass*mag(Un),                     // impingement pressure
+        mass*p.hs()                       // energy
     );
 
     this->nParcelsTransferred()++;
 
-    keepParticle = false;
+    this->totalMassTransferred() += mass;
 
+    keepParticle = false;
 }
 
 
@@ -258,9 +252,10 @@ void Foam::ThermoSurfaceFilm<CloudType>::bounceInteraction
 
 
 template<class CloudType>
+template<class filmType>
 void Foam::ThermoSurfaceFilm<CloudType>::drySplashInteraction
 (
-    //regionModels::surfaceFilmModels::surfaceFilmRegionModel& filmModel,
+    filmType& filmModel,
     const parcelType& p,
     const polyPatch& pp,
     const label facei,
@@ -301,22 +296,24 @@ void Foam::ThermoSurfaceFilm<CloudType>::drySplashInteraction
 
     if (We < Wec) // Adhesion - assume absorb
     {
-        absorbInteraction(p, pp, facei, m, keepParticle);
+        absorbInteraction<filmType>
+            (filmModel, p, pp, facei, m, keepParticle);
     }
     else // Splash
     {
         // Ratio of incident mass to splashing mass
         const scalar mRatio = 0.2 + 0.6*rndGen_.sample01<scalar>();
-        splashInteraction
-            (p, pp, facei, mRatio, We, Wec, sigma, keepParticle);
+        splashInteraction<filmType>
+            (filmModel, p, pp, facei, mRatio, We, Wec, sigma, keepParticle);
     }
 }
 
 
 template<class CloudType>
+template<class filmType>
 void Foam::ThermoSurfaceFilm<CloudType>::wetSplashInteraction
 (
-    //regionModels::surfaceFilmModels::surfaceFilmRegionModel& filmModel,
+    filmType& filmModel,
     parcelType& p,
     const polyPatch& pp,
     const label facei,
@@ -359,7 +356,8 @@ void Foam::ThermoSurfaceFilm<CloudType>::wetSplashInteraction
 
     if (We < 2) // Adhesion - assume absorb
     {
-        absorbInteraction(p, pp, facei, m, keepParticle);
+        absorbInteraction<filmType>
+            (filmModel, p, pp, facei, m, keepParticle);
     }
     else if ((We >= 2) && (We < 20)) // Bounce
     {
@@ -377,23 +375,25 @@ void Foam::ThermoSurfaceFilm<CloudType>::wetSplashInteraction
     }
     else if ((We >= 20) && (We < Wec)) // Spread - assume absorb
     {
-        absorbInteraction(p, pp, facei, m, keepParticle);
+        absorbInteraction<filmType>
+            (filmModel, p, pp, facei, m, keepParticle);
     }
     else    // Splash
     {
         // Ratio of incident mass to splashing mass
         // splash mass can be > incident mass due to film entrainment
         const scalar mRatio = 0.2 + 0.9*rndGen_.sample01<scalar>();
-        splashInteraction
-            (p, pp, facei, mRatio, We, Wec, sigma, keepParticle);
+        splashInteraction<filmType>
+            (filmModel, p, pp, facei, mRatio, We, Wec, sigma, keepParticle);
     }
 }
 
 
 template<class CloudType>
+template<class filmType>
 void Foam::ThermoSurfaceFilm<CloudType>::splashInteraction
 (
-    //regionModels::surfaceFilmModels::surfaceFilmRegionModel& filmModel,
+    filmType& filmModel,
     const parcelType& p,
     const polyPatch& pp,
     const label facei,
@@ -466,7 +466,8 @@ void Foam::ThermoSurfaceFilm<CloudType>::splashInteraction
     // Switch to absorb if insufficient energy for splash
     if (EKs <= 0)
     {
-        absorbInteraction(p, pp, facei, m, keepParticle);
+        absorbInteraction<filmType>
+            (filmModel, p, pp, facei, m, keepParticle);
         return;
     }
 
@@ -521,7 +522,8 @@ void Foam::ThermoSurfaceFilm<CloudType>::splashInteraction
     // Transfer remaining part of parcel to film 0 - splashMass can be -ve
     // if entraining from the film
     const scalar mDash = m - mSplash;
-    absorbInteraction(p, pp, facei, mDash, keepParticle);
+    absorbInteraction<filmType>
+        (filmModel, p, pp, facei, mDash, keepParticle);
 }
 
 
@@ -541,6 +543,7 @@ Foam::ThermoSurfaceFilm<CloudType>::ThermoSurfaceFilm
         owner.db().objectRegistry::template lookupObject<SLGThermo>("SLGThermo")
     ),
     filmModel_(nullptr),
+    areaFilms_(0),
     TFilmPatch_(0),
     CpFilmPatch_(0),
     interactionType_
@@ -570,17 +573,7 @@ Foam::ThermoSurfaceFilm<CloudType>::ThermoSurfaceFilm
         this->coeffDict().readEntry("Cf", Cf_);
     }
 
-    filmModel_ =
-        const_cast<regionModels::surfaceFilmModels::surfaceFilmRegionModel*>
-        (
-            this->owner().mesh().time().objectRegistry::template findObject
-            <
-                regionModels::surfaceFilmModels::surfaceFilmRegionModel
-            >
-            (
-                "surfaceFilmProperties"
-            )
-        );
+    init();
 }
 
 
@@ -594,6 +587,7 @@ Foam::ThermoSurfaceFilm<CloudType>::ThermoSurfaceFilm
     rndGen_(sfm.rndGen_),
     thermo_(sfm.thermo_),
     filmModel_(nullptr),
+    areaFilms_(0),
     TFilmPatch_(sfm.TFilmPatch_),
     CpFilmPatch_(sfm.CpFilmPatch_),
     interactionType_(sfm.interactionType_),
@@ -605,17 +599,7 @@ Foam::ThermoSurfaceFilm<CloudType>::ThermoSurfaceFilm
     Cf_(sfm.Cf_),
     nParcelsSplashed_(sfm.nParcelsSplashed_)
 {
-     filmModel_ =
-        const_cast<regionModels::surfaceFilmModels::surfaceFilmRegionModel*>
-        (
-            this->owner().mesh().time().objectRegistry::template findObject
-            <
-                regionModels::surfaceFilmModels::surfaceFilmRegionModel
-            >
-            (
-                "surfaceFilmProperties"
-            )
-        );
+    init();
 }
 
 
@@ -636,10 +620,10 @@ bool Foam::ThermoSurfaceFilm<CloudType>::transferParcel
     bool& keepParticle
 )
 {
-    const fvMesh& mesh = this->owner().mesh();
-
     const label patchi = pp.index();
 
+    bool bInteraction(false);
+
     // Check the singleLayer film models
     if (filmModel_ && filmModel_->isRegionPatch(patchi))
     {
@@ -656,7 +640,9 @@ bool Foam::ThermoSurfaceFilm<CloudType>::transferParcel
             case itAbsorb:
             {
                 const scalar m = p.nParticle()*p.mass();
-                absorbInteraction(p, pp, facei, m, keepParticle);
+
+                absorbInteraction<regionFilm>
+                    (*filmModel_, p, pp, facei, m, keepParticle);
 
                 break;
             }
@@ -666,11 +652,13 @@ bool Foam::ThermoSurfaceFilm<CloudType>::transferParcel
 
                 if (dry)
                 {
-                    drySplashInteraction(p, pp, facei, keepParticle);
+                    drySplashInteraction<regionFilm>
+                        (*filmModel_, p, pp, facei, keepParticle);
                 }
                 else
                 {
-                    wetSplashInteraction( p, pp, facei, keepParticle);
+                    wetSplashInteraction<regionFilm>
+                        (*filmModel_, p, pp, facei, keepParticle);
                 }
 
                 break;
@@ -684,59 +672,66 @@ bool Foam::ThermoSurfaceFilm<CloudType>::transferParcel
         }
 
         // Transfer parcel/parcel interactions complete
-        return true;
+        bInteraction = true;
     }
-    else
+
+
+    for (areaFilm& film : areaFilms_)
     {
-
-        // Check the finite area films
-        wordList names =
-            mesh.time().objectRegistry::template
-                sortedNames<regionModels::regionFaModel>();
-
-        forAll(names, i)
+        if (patchi == film.patchID())
         {
-            const regionModels::regionFaModel* regionFa =
-                mesh.time().objectRegistry::template findObject
-                <
-                    regionModels::regionFaModel
-                >(names[i]);
+            const label facei = pp.whichFace(p.face());
 
-            // Check that regionFa is a areaFilm
-            if (regionFa && isA<areaFilm>(*regionFa))
+            switch (interactionType_)
             {
-                areaFilm& film =
-                    const_cast<areaFilm&>(refCast<const areaFilm>(*regionFa));
-
-                const label facei = pp.whichFace(p.face());
-
-                switch (interactionType_)
+                // It only supports absorp model
+                case itAbsorb:
                 {
-                    // It only supports absorp model
-                    case itAbsorb:
-                    {
-                        const scalar m = p.nParticle()*p.mass();
-                        absorbInteraction
-                        (
-                            film, p, pp, facei, m, keepParticle
-                        );
-                        break;
-                    }
-                    default:
-                    {
-                        FatalErrorInFunction
-                            << "Unknown interaction type enumeration"
-                            << abort(FatalError);
-                    }
+                    const scalar m = p.nParticle()*p.mass();
+
+                    absorbInteraction<areaFilm>
+                    (
+                        film, p, pp, facei, m, keepParticle
+                    );
+                    break;
+                }
+                case itBounce:
+                {
+                    bounceInteraction(p, pp, facei, keepParticle);
+
+                    break;
+                }
+                case itSplashBai:
+                {
+                    bool dry = film.h()[facei] < deltaWet_;
+
+                    if (dry)
+                    {
+                        drySplashInteraction<areaFilm>
+                            (film, p, pp, facei, keepParticle);
+                    }
+                    else
+                    {
+                        wetSplashInteraction<areaFilm>
+                            (film, p, pp, facei, keepParticle);
+                    }
+
+                    break;
+                }
+                default:
+                {
+                    FatalErrorInFunction
+                        << "Unknown interaction type enumeration"
+                        << abort(FatalError);
                 }
-                // Transfer parcel/parcel interactions complete
-                return true;
             }
+            // Transfer parcel/parcel interactions complete
+            bInteraction = true;
         }
     }
 
     // Parcel not interacting with film
-    return false;
+    return bInteraction;
 }
 
 
@@ -779,7 +774,8 @@ void Foam::ThermoSurfaceFilm<CloudType>::cacheFilmFields
 
     TFilmPatch_.setSize(filmModel.Tf().size(), Zero);
     map.mapToField(filmModel.Tf(), TFilmPatch_);
-    //map.mapToField(filmModel.Cp(), CpFilmPatch_);
+    CpFilmPatch_.setSize(filmModel.Tf().size(), Zero);
+    map.mapToField(filmModel.Cp(), CpFilmPatch_);
 }
 
 
diff --git a/src/lagrangian/intermediate/submodels/Thermodynamic/SurfaceFilmModel/ThermoSurfaceFilm/ThermoSurfaceFilm.H b/src/lagrangian/intermediate/submodels/Thermodynamic/SurfaceFilmModel/ThermoSurfaceFilm/ThermoSurfaceFilm.H
index b1b821a294..5e5221062b 100644
--- a/src/lagrangian/intermediate/submodels/Thermodynamic/SurfaceFilmModel/ThermoSurfaceFilm/ThermoSurfaceFilm.H
+++ b/src/lagrangian/intermediate/submodels/Thermodynamic/SurfaceFilmModel/ThermoSurfaceFilm/ThermoSurfaceFilm.H
@@ -55,6 +55,7 @@ SourceFiles
 #define ThermoSurfaceFilm_H
 
 #include "SurfaceFilmModel.H"
+#include "UPtrList.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -121,6 +122,9 @@ protected:
         typedef typename
             regionModels::areaSurfaceFilmModels::liquidFilmBase areaFilm;
 
+        typedef typename
+            regionModels::surfaceFilmModels::surfaceFilmRegionModel regionFilm;
+
         //- Reference to the cloud random number generator
         Random& rndGen_;
 
@@ -128,8 +132,10 @@ protected:
         const SLGThermo& thermo_;
 
         //- Pointer to filmModel
-        regionModels::surfaceFilmModels::surfaceFilmRegionModel* filmModel_;
+        regionFilm* filmModel_;
 
+        //- UPointers to area films
+        UPtrList<areaFilm> areaFilms_;
 
         // Cached injector fields per film patch
 
@@ -189,24 +195,17 @@ protected:
             const vector& nf
         ) const;
 
+        //- Init films pointers
+        void init();
+
 
         // Interaction models
 
             //- Absorb parcel into film
+            template<class filmType>
             void absorbInteraction
             (
-                //regionModels::surfaceFilmModels::surfaceFilmRegionModel&,
-                const parcelType& p,
-                const polyPatch& pp,
-                const label facei,
-                const scalar mass,
-                bool& keepParticle
-            );
-
-            //- Absorb parcel into area film
-            void absorbInteraction
-            (
-                areaFilm&,
+                filmType&,
                 const parcelType& p,
                 const polyPatch& pp,
                 const label facei,
@@ -224,9 +223,10 @@ protected:
             ) const;
 
             //- Parcel interaction with dry surface
+            template<class filmType>
             void drySplashInteraction
             (
-                //regionModels::surfaceFilmModels::surfaceFilmRegionModel&,
+                filmType&,
                 const parcelType& p,
                 const polyPatch& pp,
                 const label facei,
@@ -234,9 +234,10 @@ protected:
             );
 
             //- Parcel interaction with wetted surface
+            template<class filmType>
             void wetSplashInteraction
             (
-                //regionModels::surfaceFilmModels::surfaceFilmRegionModel&,
+                filmType&,
                 parcelType& p,
                 const polyPatch& pp,
                 const label facei,
@@ -244,9 +245,10 @@ protected:
             );
 
             //- Bai parcel splash interaction model
+            template<class filmType>
             void splashInteraction
             (
-                //regionModels::surfaceFilmModels::surfaceFilmRegionModel&,
+                filmType&,
                 const parcelType& p,
                 const polyPatch& pp,
                 const label facei,
diff --git a/src/regionFaModels/Make/files b/src/regionFaModels/Make/files
index 70521f391a..85634c54ee 100644
--- a/src/regionFaModels/Make/files
+++ b/src/regionFaModels/Make/files
@@ -15,9 +15,25 @@ derivedFvPatchFields/vibrationShell/vibrationShellFvPatchScalarField.C
 
 liquidFilm/liquidFilmBase.C
 liquidFilm/liquidFilmBaseNew.C
+
+/* Sub-Model */
+
 liquidFilm/subModels/kinematic/filmTurbulenceModel/filmTurbulenceModel/filmTurbulenceModel.C
 liquidFilm/subModels/kinematic/filmTurbulenceModel/filmTurbulenceModel/filmTurbulenceModelNew.C
 liquidFilm/subModels/kinematic/filmTurbulenceModel/laminar/laminar.C
+
+liquidFilm/subModels/kinematic/injectionModel/injectionModelList/injectionModelList.C
+
+liquidFilm/subModels/kinematic/injectionModel/injectionModel/injectionModel.C
+liquidFilm/subModels/kinematic/injectionModel/injectionModel/injectionModelNew.C
+
+liquidFilm/subModels/kinematic/injectionModel/curvatureSeparation/curvatureSeparation.C
+
+liquidFilm/subModels/filmSubModelBase.C
+
+liquidFilm/liquidFilmBase.C
+liquidFilm/liquidFilmBaseNew.C
+
 liquidFilm/liquidFilmModel/liquidFilmModel.C
 liquidFilm/kinematicThinFilm/kinematicThinFilm.C
 derivedFvPatchFields/filmShell/velocityFilmShellFvPatchVectorField.C
diff --git a/src/regionFaModels/liquidFilm/liquidFilmBase.C b/src/regionFaModels/liquidFilm/liquidFilmBase.C
index 540bbe9da9..a1d36f1f55 100644
--- a/src/regionFaModels/liquidFilm/liquidFilmBase.C
+++ b/src/regionFaModels/liquidFilm/liquidFilmBase.C
@@ -28,7 +28,6 @@ License
 #include "liquidFilmBase.H"
 #include "faMesh.H"
 #include "faCFD.H"
-#include "uniformDimensionedFields.H"
 #include "gravityMeshObject.H"
 #include "movingWallVelocityFvPatchVectorField.H"
 #include "turbulentFluidThermoModel.H"
@@ -60,11 +59,11 @@ bool liquidFilmBase::read(const dictionary& dict)
     regionFaModel::read(dict);
     if (active_)
     {
-        const dictionary& solution = this->solution().subDict("PISO");
+        const dictionary& solution = this->solution().subDict("PIMPLE");
         solution.readEntry("momentumPredictor", momentumPredictor_);
-        solution.readIfPresent("nOuterCorr", nOuterCorr_);
+        solution.readEntry("nOuterCorr", nOuterCorr_);
         solution.readEntry("nCorr", nCorr_);
-        solution.readEntry("nNonOrthCorr", nNonOrthCorr_);
+        solution.readEntry("nFilmCorr", nFilmCorr_);
     }
     return true;
 }
@@ -73,28 +72,24 @@ bool liquidFilmBase::read(const dictionary& dict)
 scalar liquidFilmBase::CourantNumber() const
 {
     scalar CoNum = 0.0;
-    scalar meanCoNum = 0.0;
     scalar velMag = 0.0;
 
-    if (regionMesh().nInternalEdges())
-    {
-        edgeScalarField SfUfbyDelta
-        (
-            regionMesh().edgeInterpolation::deltaCoeffs()*mag(phif_)
-        );
+    edgeScalarField SfUfbyDelta
+    (
+        regionMesh().edgeInterpolation::deltaCoeffs()*mag(phif_)
+    );
 
-        CoNum = max(SfUfbyDelta/regionMesh().magLe())
-            .value()*time().deltaT().value();
+    CoNum = max(SfUfbyDelta/regionMesh().magLe())
+        .value()*time().deltaT().value();
 
-        meanCoNum = (sum(SfUfbyDelta)/sum(regionMesh().magLe()))
-            .value()*time().deltaT().value();
+    velMag = max(mag(phif_)/regionMesh().magLe()).value();
 
-        velMag = max(mag(phif_)/regionMesh().magLe()).value();
-    }
+    reduce(CoNum, maxOp<scalar>());
+    reduce(velMag, maxOp<scalar>());
 
-    Info<< "Film Courant Number mean: " << meanCoNum
+    Info<< "Film Courant Number: "
         << " max: " << CoNum
-        << " Film velocity magnitude: " << velMag << endl;
+        << " Film velocity magnitude: (h)" << velMag << endl;
 
     return CoNum;
 }
@@ -114,19 +109,21 @@ liquidFilmBase::liquidFilmBase
 
     momentumPredictor_
     (
-        this->solution().subDict("PISO").lookup("momentumPredictor")
+        this->solution().subDict("PIMPLE").get<bool>("momentumPredictor")
     ),
     nOuterCorr_
     (
-        this->solution().subDict("PISO").lookupOrDefault("nOuterCorr", 1)
+        this->solution().subDict("PIMPLE").get<label>("nOuterCorr")
     ),
-    nCorr_(this->solution().subDict("PISO").get<label>("nCorr")),
-    nNonOrthCorr_
+    nCorr_(this->solution().subDict("PIMPLE").get<label>("nCorr")),
+    nFilmCorr_
     (
-        this->solution().subDict("PISO").get<label>("nNonOrthCorr")
+        this->solution().subDict("PIMPLE").get<label>("nFilmCorr")
     ),
 
-    h0_("h0", dimLength, SMALL),
+    h0_("h0", dimLength, 1e-7, dict),
+
+    deltaWet_("deltaWet", dimLength, 1e-4, dict),
 
     UName_(dict.get<word>("U")),
 
@@ -211,18 +208,6 @@ liquidFilmBase::liquidFilmBase
         fac::interpolate(h_*Uf_) & regionMesh().Le()
     ),
 
-    Tf_
-    (
-        IOobject
-        (
-            "Tf_" + regionName_,
-            primaryMesh().time().timeName(),
-            primaryMesh(),
-            IOobject::MUST_READ,
-            IOobject::AUTO_WRITE
-        ),
-        regionMesh()
-    ),
 
     gn_
     (
@@ -238,6 +223,8 @@ liquidFilmBase::liquidFilmBase
         dimensionedScalar(dimAcceleration, Zero)
     ),
 
+    g_(meshObjects::gravity::New(primaryMesh().time())),
+
     massSource_
     (
         IOobject
@@ -294,10 +281,9 @@ liquidFilmBase::liquidFilmBase
 
     faOptions_(Foam::fa::options::New(p))
 {
-    const uniformDimensionedVectorField& g =meshObjects::gravity::New(time());
-
-    gn_ = (g & regionMesh().faceAreaNormals());
+    const areaVectorField& ns = regionMesh().faceAreaNormals();
 
+    gn_ = g_ & ns;
 
     if (!faOptions_.optionList::size())
     {
@@ -407,37 +393,7 @@ Foam::tmp<Foam::areaVectorField> liquidFilmBase::Up() const
     );
 
     areaVectorField& Up = tUp.ref();
-/*
-    typedef compressible::turbulenceModel cmpTurbModelType;
-    typedef incompressible::turbulenceModel incmpTurbModelType;
 
-    word turbName
-    (
-        IOobject::groupName
-        (
-            turbulenceModel::propertiesName,
-            Up_.internalField().group()
-        )
-    );
-
-    scalarField nu(patch_.size(), Zero);
-    if (primaryMesh().foundObject<cmpTurbModelType>(turbName))
-    {
-        const cmpTurbModelType& turbModel =
-            primaryMesh().lookupObject<cmpTurbModelType>(turbName);
-
-        //const basicThermo& thermo = turbModel.transport();
-
-        nu = turbModel.nu(patchi);
-    }
-    if (primaryMesh().foundObject<incmpTurbModelType>(turbName))
-    {
-        const incmpTurbModelType& turbModel =
-            primaryMesh().lookupObject<incmpTurbModelType>(turbName);
-
-        nu = turbModel.nu(patchi);
-    }
-*/
     scalarField hp(patch_.size(), Zero);
 
     // map areas h to hp on primary
@@ -481,7 +437,10 @@ tmp<areaScalarField> liquidFilmBase::pg() const
         const volScalarField& pp =
             primaryMesh().lookupObject<volScalarField>(pName_);
 
-        const volScalarField::Boundary& pw = pp.boundaryField();
+        volScalarField::Boundary& pw =
+            const_cast<volScalarField::Boundary&>(pp.boundaryField());
+
+        //pw -= pRef_;
 
         pfg.primitiveFieldRef() = vsmPtr_->mapInternalToSurface<scalar>(pw)();
     }
@@ -490,6 +449,30 @@ tmp<areaScalarField> liquidFilmBase::pg() const
 }
 
 
+tmp<areaScalarField> liquidFilmBase::alpha() const
+{
+    tmp<areaScalarField> talpha
+    (
+        new areaScalarField
+        (
+            IOobject
+            (
+                "talpha",
+                primaryMesh().time().timeName(),
+                primaryMesh()
+            ),
+            regionMesh(),
+            dimensionedScalar(dimless, Zero)
+        )
+    );
+    areaScalarField& alpha = talpha.ref();
+
+    alpha = pos0(h_ - deltaWet_);
+
+    return talpha;
+}
+
+
 void liquidFilmBase::addSources
 (
     const label patchi,
@@ -507,46 +490,65 @@ void liquidFilmBase::addSources
 
     momentumSource_.boundaryFieldRef()[patchi][facei] += momentumSource;
 
-    addedMassTotal_ += massSource;
-
-    if (debug)
-    {
-        InfoInFunction
-            << "\nSurface film: " << type() << ": adding to film source:" << nl
-            << "    mass     = " << massSource << nl
-            << "    momentum = " << momentumSource << nl
-            << "    pressure = " << pressureSource << endl;
-    }
 }
 
 
+void liquidFilmBase::preEvolveRegion()
+{
+    regionFaModel::preEvolveRegion();
+}
+
+void liquidFilmBase::postEvolveRegion()
+{
+    if (debug && primaryMesh().time().writeTime())
+    {
+        massSource_.write();
+        pnSource_.write();
+        momentumSource_.write();
+    }
+
+    massSource_.boundaryFieldRef() = Zero;
+    pnSource_.boundaryFieldRef() = Zero;
+    momentumSource_.boundaryFieldRef() = Zero;
+
+    regionFaModel::postEvolveRegion();
+}
+
 Foam::fa::options& liquidFilmBase::faOptions()
 {
      return faOptions_;
 }
 
-
 const areaVectorField& liquidFilmBase::Uf() const
 {
      return Uf_;
 }
 
-
-const areaScalarField& liquidFilmBase::Tf() const
-{
-     return Tf_;
-}
-
 const areaScalarField& liquidFilmBase::gn() const
 {
      return gn_;
 }
 
+const uniformDimensionedVectorField& liquidFilmBase::g() const
+{
+    return g_;
+}
+
 const areaScalarField& liquidFilmBase::h() const
 {
      return h_;
 }
 
+const edgeScalarField& liquidFilmBase::phif() const
+{
+    return phif_;
+}
+
+const edgeScalarField& liquidFilmBase::phi2s() const
+{
+    return phi2s_;
+}
+
 const dimensionedScalar& liquidFilmBase::h0() const
 {
      return h0_;
@@ -557,6 +559,7 @@ const regionFaModel& liquidFilmBase::region() const
     return *this;
 }
 
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 } // End namespace areaSurfaceFilmModels
diff --git a/src/regionFaModels/liquidFilm/liquidFilmBase.H b/src/regionFaModels/liquidFilm/liquidFilmBase.H
index e0f9522fee..172f15dbfc 100644
--- a/src/regionFaModels/liquidFilm/liquidFilmBase.H
+++ b/src/regionFaModels/liquidFilm/liquidFilmBase.H
@@ -41,10 +41,10 @@ SourceFiles
 #include "autoPtr.H"
 #include "areaFieldsFwd.H"
 #include "volFieldsFwd.H"
+#include "uniformDimensionedFields.H"
 #include "regionFaModel.H"
 #include "faOptions.H"
 
-
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 namespace Foam
@@ -54,7 +54,6 @@ namespace regionModels
 namespace areaSurfaceFilmModels
 {
 
-
 /*---------------------------------------------------------------------------*\
                          Class liquidFilmBase Declaration
 \*---------------------------------------------------------------------------*/
@@ -92,15 +91,19 @@ protected:
             //- Number of PISO-like correctors
             label nCorr_;
 
-            //- Number of non-orthogonal correctors
-            label nNonOrthCorr_;
+            //- Number of film thickness correctors
+            label nFilmCorr_;
 
             //- Cumulative continuity error
             scalar cumulativeContErr_;
 
-            //- Small delta
+            //- Smallest numerical thickness
             dimensionedScalar h0_;
 
+            //- Delta wet for sub-models
+            dimensionedScalar deltaWet_;
+
+
         //- Name of the velocity field
         word UName_;
 
@@ -131,14 +134,14 @@ protected:
             //- Film height flux
             edgeScalarField phi2s_;
 
-            //- Film temperature
-            areaScalarField Tf_;
-
             //- Normal gravity field
             areaScalarField gn_;
 
+            //- Gravity
+            uniformDimensionedVectorField g_;
 
-            // Mass exchanage fields
+
+            // Mass exchanage fields from the primary region (lagragian)
 
                 //- Mass
                 volScalarField massSource_;
@@ -156,10 +159,12 @@ protected:
                 scalar addedMassTotal_;
 
 
+
         //- Pointer to faOptions
         Foam::fa::options& faOptions_;
 
 
+
     // Protected Member Functions
 
         //- Read control parameters from dictionary
@@ -235,6 +240,10 @@ public:
             //- Map primary static pressure
             tmp<areaScalarField> pg() const;
 
+            //- Wet indicator using h0
+            tmp<areaScalarField> alpha() const;
+
+
         // Access functions
 
             //- Return faOptions
@@ -243,15 +252,22 @@ public:
             //- Access const reference Uf
             const areaVectorField& Uf() const;
 
-            //- Access const reference Tf
-            const areaScalarField& Tf() const;
-
             //- Access const reference gn
             const areaScalarField& gn() const;
 
+            //- Gravity
+            const uniformDimensionedVectorField& g() const;
+
             //- Access const reference h
             const areaScalarField& h() const;
 
+            //- Access to momentum flux
+            const edgeScalarField& phif() const;
+
+            //- Access continuity flux
+            const edgeScalarField& phi2s() const;
+
+
             //- Return h0
             const dimensionedScalar& h0() const;
 
@@ -259,6 +275,28 @@ public:
             const regionFaModel& region() const;
 
 
+            // Transfer fields - to the primary region (lagragian injection)
+
+                //- Return mass transfer source - Eulerian phase only
+                //virtual tmp<volScalarField> primaryMassTrans() const = 0;
+
+                //- Return the film mass available for transfer to cloud
+                virtual const volScalarField& cloudMassTrans() const = 0;
+
+                //- Return the parcel diameters originating from film to cloud
+                virtual const volScalarField& cloudDiameterTrans() const = 0;
+
+
+
+        // Evolution
+
+            //- Pre-evolve film
+            virtual void preEvolveRegion();
+
+            //- Post-evolve film
+            virtual void postEvolveRegion();
+
+
         // Thermo variables
 
              //- Access const reference mu
@@ -270,11 +308,16 @@ public:
             //- Access const reference sigma
             virtual const areaScalarField& sigma() const = 0;
 
+            //- Access const reference Tf
+            virtual const areaScalarField& Tf() const = 0;
+
+            //- Access const reference Cp
+            virtual const areaScalarField& Cp() const = 0;
+
 
         // External hook to add sources and mass exchange variables
 
 
-
             //- Add sources
             virtual void addSources
             (
@@ -285,7 +328,6 @@ public:
                 const scalar pressureSource,   // [kg.m/s] (normal momentum)
                 const scalar energySource = 0  // [J]
             );
-
 };
 
 
diff --git a/src/regionFaModels/liquidFilm/liquidFilmModel/liquidFilmModel.C b/src/regionFaModels/liquidFilm/liquidFilmModel/liquidFilmModel.C
index b4a07ededc..306a09b895 100644
--- a/src/regionFaModels/liquidFilm/liquidFilmModel/liquidFilmModel.C
+++ b/src/regionFaModels/liquidFilm/liquidFilmModel/liquidFilmModel.C
@@ -31,6 +31,7 @@ License
 #include "gravityMeshObject.H"
 #include "volFields.H"
 
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 namespace Foam
@@ -61,6 +62,7 @@ void liquidFilmModel::correctThermoFields()
         rho_[faceI] = thermo_.rho(pRef_, Tf_[faceI], X);
         mu_[faceI] = thermo_.mu(pRef_, Tf_[faceI], X);
         sigma_[faceI] = thermo_.sigma(pRef_, Tf_[faceI], X);
+        Cp_[faceI] = thermo_.Cp(pRef_, Tf_[faceI], X);
     }
 
     forAll (regionMesh().boundary(), patchI)
@@ -70,14 +72,19 @@ void liquidFilmModel::correctThermoFields()
         scalarField& patchRho = rho_.boundaryFieldRef()[patchI];
         scalarField& patchmu = mu_.boundaryFieldRef()[patchI];
         scalarField& patchsigma = sigma_.boundaryFieldRef()[patchI];
+        scalarField& patchCp = Cp_.boundaryFieldRef()[patchI];
 
         forAll(patchRho, edgeI)
         {
             patchRho[edgeI] = thermo_.rho(pRef_, patchTf[edgeI], X);
             patchmu[edgeI] = thermo_.mu(pRef_, patchTf[edgeI], X);
             patchsigma[edgeI] = thermo_.sigma(pRef_, patchTf[edgeI], X);
+            patchCp[edgeI] = thermo_.Cp(pRef_, patchTf[edgeI], X);
         }
     }
+
+    //Initialize pf_
+    pf_ = rho_*gn_*h_ - sigma_*fac::laplacian(h_);
 }
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
@@ -119,6 +126,32 @@ liquidFilmModel::liquidFilmModel
         regionMesh(),
         dimensionedScalar(dimViscosity, Zero)
     ),
+    Tf_
+    (
+        IOobject
+        (
+            "Tf_" + regionName_,
+            primaryMesh().time().timeName(),
+            primaryMesh(),
+            IOobject::READ_IF_PRESENT,
+            IOobject::AUTO_WRITE
+        ),
+        regionMesh(),
+        dimensionedScalar(dimTemperature, Zero)
+    ),
+    Cp_
+    (
+        IOobject
+        (
+            "Cp_" + regionName_,
+            primaryMesh().time().timeName(),
+            primaryMesh(),
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        regionMesh(),
+        dimensionedScalar(dimEnergy/dimTemperature, Zero)
+    ),
     sigma_
     (
         IOobject
@@ -151,7 +184,9 @@ liquidFilmModel::liquidFilmModel
         (
             "rhoSp",
             primaryMesh().time().timeName(),
-            primaryMesh()
+            primaryMesh(),
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
         ),
         regionMesh(),
         dimensionedScalar(dimVelocity, Zero)
@@ -178,8 +213,49 @@ liquidFilmModel::liquidFilmModel
         regionMesh(),
         dimensionedScalar(dimPressure, Zero)
     ),
-    turbulence_(filmTurbulenceModel::New(*this, dict))
+    cloudMassTrans_
+    (
+        IOobject
+        (
+            "cloudMassTrans",
+            primaryMesh().time().timeName(),
+            primaryMesh(),
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        primaryMesh(),
+        dimensionedScalar(dimMass, Zero),
+        calculatedFvPatchField<scalar>::typeName
+    ),
+
+    cloudDiameterTrans_
+    (
+        IOobject
+        (
+            "cloudDiameterTrans",
+            primaryMesh().time().timeName(),
+            primaryMesh(),
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        primaryMesh(),
+        dimensionedScalar(dimLength, Zero),
+        calculatedFvPatchField<scalar>::typeName
+    ),
+
+    turbulence_(filmTurbulenceModel::New(*this, dict)),
+
+    availableMass_(regionMesh().faces().size(), Zero),
+
+    injection_(*this, dict),
+
+    forces_(*this, dict)
 {
+
+    if (dict.found("T0"))
+    {
+        Tf_ = dimensionedScalar("T0", dimTemperature, dict);
+    }
     correctThermoFields();
 }
 
@@ -208,8 +284,39 @@ const areaScalarField& liquidFilmModel::sigma() const
      return sigma_;
 }
 
+
+const areaScalarField& liquidFilmModel::Tf() const
+{
+     return Tf_;
+}
+
+
+const areaScalarField& liquidFilmModel::Cp() const
+{
+     return Cp_;
+}
+
+const volScalarField& liquidFilmModel::cloudMassTrans() const
+{
+    return cloudMassTrans_;
+}
+
+const volScalarField& liquidFilmModel::cloudDiameterTrans() const
+{
+    return cloudDiameterTrans_;
+}
+
+
+
 void liquidFilmModel::preEvolveRegion()
 {
+    liquidFilmBase::preEvolveRegion();
+
+
+
+    cloudMassTrans_ == dimensionedScalar(dimMass, Zero);
+    cloudDiameterTrans_ == dimensionedScalar(dimLength, Zero);
+
     const scalar deltaT = primaryMesh().time().deltaTValue();
     const scalarField rAreaDeltaT = 1/deltaT/regionMesh().S().field();
 
@@ -219,6 +326,7 @@ void liquidFilmModel::preEvolveRegion()
     // [kg.m/s]
     USp_.primitiveFieldRef() =
         vsm().mapToSurface(momentumSource_.boundaryField()[patchID()]);
+
     pnSp_.primitiveFieldRef() =
         vsm().mapToSurface(pnSource_.boundaryField()[patchID()]);
 
@@ -227,13 +335,38 @@ void liquidFilmModel::preEvolveRegion()
     rhoSp_.primitiveFieldRef() *= rAreaDeltaT/rho_;
     USp_.primitiveFieldRef() *= rAreaDeltaT/rho_;
     pnSp_.primitiveFieldRef() *= rAreaDeltaT/rho_;
+
+    rhoSp_.relax();
+    pnSp_.relax();
+    USp_.relax();
 }
 
 void liquidFilmModel::postEvolveRegion()
 {
-    massSource_ == dimensionedScalar(massSource_.dimensions(), Zero);
-    momentumSource_ == dimensionedVector(momentumSource_.dimensions(), Zero);
-    pnSource_ == dimensionedScalar(pnSource_.dimensions(), Zero);
+    availableMass_ = (h() - h0_)*rho()*regionMesh().S();
+    injection_.correct(availableMass_, cloudMassTrans_, cloudDiameterTrans_);
+    liquidFilmBase::postEvolveRegion();
+}
+
+
+void liquidFilmModel::info()
+{
+    Info<< "\nSurface film: " << type() << " on patch: " << patchID() << endl;
+
+    const DimensionedField<scalar, areaMesh>& sf = regionMesh().S();
+
+    Info<< indent << "min/max(mag(Uf))    = " << gMin(mag(Uf_.field())) << ", "
+        << gMax(mag(Uf_.field())) << nl
+        << indent << "min/max(delta)     = " << gMin(h_.field()) << ", " << gMax(h_.field()) << nl
+        << indent << "coverage           = "
+        << gSum(alpha()().field()*mag(sf.field()))/gSum(mag(sf.field())) <<  nl
+        << indent << "total mass         = "
+        << gSum(availableMass_) << nl;
+
+
+    Info<< indent << CourantNumber() << endl;
+
+    injection_.info(Info);
 }
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
diff --git a/src/regionFaModels/liquidFilm/liquidFilmModel/liquidFilmModel.H b/src/regionFaModels/liquidFilm/liquidFilmModel/liquidFilmModel.H
index 0aeebf6118..89225b70ca 100644
--- a/src/regionFaModels/liquidFilm/liquidFilmModel/liquidFilmModel.H
+++ b/src/regionFaModels/liquidFilm/liquidFilmModel/liquidFilmModel.H
@@ -42,8 +42,16 @@ SourceFiles
 #include "volFieldsFwd.H"
 #include "liquidFilmBase.H"
 #include "faMesh.H"
+
 #include "filmTurbulenceModel.H"
+
 #include "liquidMixtureProperties.H"
+#include "injectionModelList.H"
+
+#include "faCFD.H"
+
+//#include "transferModelList.H"
+#include "forceList.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -54,6 +62,7 @@ namespace regionModels
 namespace areaSurfaceFilmModels
 {
 
+
 /*---------------------------------------------------------------------------*\
                         Class liquidFilmModel Declaration
 \*---------------------------------------------------------------------------*/
@@ -94,6 +103,12 @@ protected:
                 //- Dynamic viscosity [Pa.s]
                 areaScalarField mu_;
 
+                //- Film temperature
+                areaScalarField Tf_;
+
+                 //- Film Heat capacity
+                areaScalarField Cp_;
+
                 //- Surface tension [m/s2]
                 areaScalarField sigma_;
 
@@ -112,6 +127,15 @@ protected:
                     areaScalarField pnSp_;
 
 
+                 // Transfer fields
+
+                    //- Film mass  for transfer to cloud
+                    volScalarField cloudMassTrans_;
+
+                    //- Parcel diameters originating from film to cloud
+                    volScalarField cloudDiameterTrans_;
+
+
     // General properties
 
 
@@ -120,6 +144,21 @@ protected:
             autoPtr<filmTurbulenceModel> turbulence_;
 
 
+            // Sub-models
+
+                //- Available mass for transfer via sub-models
+                scalarField availableMass_;
+
+                //- Cloud injection
+                injectionModelList injection_;
+
+                //- Transfer with the continuous phase
+                //transferModelList transfer_;
+
+                //- List of film forces
+                forceList forces_;
+
+
     // Protected member functions
 
 
@@ -169,6 +208,20 @@ public:
             //- Access const reference sigma
             const areaScalarField& sigma() const;
 
+            //- Access const reference Tf
+            const areaScalarField& Tf() const;
+
+            //- Access const reference Cp
+            const areaScalarField& Cp() const;
+
+
+            // Transfer fields - to the primary region (lagragian injection)
+
+                //- Return the film mass available for transfer to cloud
+                virtual const volScalarField& cloudMassTrans() const;
+
+                //- Return the parcel diameters originating from film to cloud
+                virtual const volScalarField& cloudDiameterTrans() const;
 
         // Evolution
 
@@ -185,7 +238,7 @@ public:
        // I-O
 
             //- Provide some feedback
-            //virtual void info();
+            virtual void info();
 };
 
 
diff --git a/src/regionFaModels/liquidFilm/subModels/filmSubModelBase.C b/src/regionFaModels/liquidFilm/subModels/filmSubModelBase.C
new file mode 100644
index 0000000000..9b072df6db
--- /dev/null
+++ b/src/regionFaModels/liquidFilm/subModels/filmSubModelBase.C
@@ -0,0 +1,110 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2021 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "filmSubModelBase.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace regionModels
+{
+namespace areaSurfaceFilmModels
+{
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+filmSubModelBase::filmSubModelBase(liquidFilmBase& film)
+:
+    subModelBase(film.outputProperties()),
+    filmModel_(film)
+{}
+
+
+filmSubModelBase::filmSubModelBase
+(
+    liquidFilmBase& film,
+    const dictionary& dict,
+    const word& baseName,
+    const word& modelType,
+    const word& dictExt
+)
+:
+    subModelBase
+    (
+        film.outputProperties(),
+        dict,
+        baseName,
+        modelType,
+        dictExt
+    ),
+    filmModel_(film)
+{}
+
+
+filmSubModelBase::filmSubModelBase
+(
+    const word& modelName,
+    liquidFilmBase& film,
+    const dictionary& dict,
+    const word& baseName,
+    const word& modelType
+)
+:
+    subModelBase
+    (
+        modelName,
+        film.outputProperties(),
+        dict,
+        baseName,
+        modelType
+    ),
+    filmModel_(film)
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+filmSubModelBase::~filmSubModelBase()
+{}
+
+
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+
+bool filmSubModelBase::writeTime() const
+{
+    return active() && filmModel_.time().writeTime();
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace areaSurfaceFilmModels
+} // End namespace regionModels
+} // End namespace Foam
+
+// ************************************************************************* //
diff --git a/src/regionFaModels/liquidFilm/subModels/filmSubModelBase.H b/src/regionFaModels/liquidFilm/subModels/filmSubModelBase.H
new file mode 100644
index 0000000000..12afa909d5
--- /dev/null
+++ b/src/regionFaModels/liquidFilm/subModels/filmSubModelBase.H
@@ -0,0 +1,139 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2021 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::regionModels::surfaceFilmModels::filmSubModelBase
+
+Description
+    Base class for surface film sub-models
+
+SourceFiles
+    filmSubModelBaseI.H
+    filmSubModelBase.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef filmSubModelBase_H
+#define filmSubModelBase_H
+
+#include "liquidFilmBase.H"
+#include "subModelBase.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace regionModels
+{
+namespace areaSurfaceFilmModels
+{
+
+/*---------------------------------------------------------------------------*\
+                      Class filmSubModelBase Declaration
+\*---------------------------------------------------------------------------*/
+
+class filmSubModelBase
+:
+    public subModelBase
+{
+protected:
+
+    // Protected Data
+
+        //- Reference to the film surface film model
+        liquidFilmBase& filmModel_;
+
+
+public:
+
+    // Constructors
+
+        //- Construct null
+        filmSubModelBase(liquidFilmBase& film);
+
+        //- Construct from film film without name
+        filmSubModelBase
+        (
+            liquidFilmBase& film,
+            const dictionary& dict,
+            const word& baseName,
+            const word& modelType,
+            const word& dictExt = "Coeffs"
+        );
+
+        //- Construct from film film with name
+        filmSubModelBase
+        (
+            const word& modelName,
+            liquidFilmBase& film,
+            const dictionary& dict,
+            const word& baseName,
+            const word& modelType
+        );
+
+
+    //- Destructor
+    virtual ~filmSubModelBase();
+
+
+    // Member Functions
+
+        // Access
+
+            //- Flag to indicate when to write a property
+            virtual bool writeTime() const;
+
+            //- Return const access to the film surface film model
+            inline const liquidFilmBase& film() const;
+
+            //- Return the reference to the film surface film model
+            inline liquidFilmBase& film();
+
+            template<class FilmType>
+            inline const FilmType& filmType() const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace areaSurfaceFilmModels
+} // End namespace regionModels
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#include "filmSubModelBaseI.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+    #include "filmSubModelBaseTemplates.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/regionFaModels/liquidFilm/subModels/filmSubModelBaseI.H b/src/regionFaModels/liquidFilm/subModels/filmSubModelBaseI.H
new file mode 100644
index 0000000000..70f07eb03f
--- /dev/null
+++ b/src/regionFaModels/liquidFilm/subModels/filmSubModelBaseI.H
@@ -0,0 +1,57 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2021 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace regionModels
+{
+namespace areaSurfaceFilmModels
+{
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+inline const liquidFilmBase& filmSubModelBase::film() const
+{
+    return filmModel_;
+}
+
+
+inline liquidFilmBase& filmSubModelBase::film()
+{
+    return filmModel_;
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace areaSurfaceFilmModels
+} // End namespace regionModels
+} // End namespace Foam
+
+// ************************************************************************* //
diff --git a/src/regionFaModels/liquidFilm/subModels/filmSubModelBaseTemplates.C b/src/regionFaModels/liquidFilm/subModels/filmSubModelBaseTemplates.C
new file mode 100644
index 0000000000..355f1a8851
--- /dev/null
+++ b/src/regionFaModels/liquidFilm/subModels/filmSubModelBaseTemplates.C
@@ -0,0 +1,58 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2021 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+namespace Foam
+{
+namespace regionModels
+{
+namespace areaSurfaceFilmModels
+{
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+template<class FilmType>
+const FilmType& filmSubModelBase::filmType() const
+{
+    if (!isA<FilmType>(filmModel_))
+    {
+        FatalErrorInFunction
+            << "Model " << this->modelType() << " requested film type "
+            << FilmType::typeName << " but film is type " << filmModel_.type()
+            << abort(FatalError);
+    }
+
+    return refCast<const FilmType>(filmModel_);
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace areaSurfaceFilmModels
+} // End namespace regionModels
+} // End namespace Foam
+
+// ************************************************************************* //
diff --git a/src/regionFaModels/liquidFilm/subModels/kinematic/filmTurbulenceModel/filmTurbulenceModel/filmTurbulenceModel.C b/src/regionFaModels/liquidFilm/subModels/kinematic/filmTurbulenceModel/filmTurbulenceModel/filmTurbulenceModel.C
index 6c0d5a2731..cec1b7acb9 100644
--- a/src/regionFaModels/liquidFilm/subModels/kinematic/filmTurbulenceModel/filmTurbulenceModel/filmTurbulenceModel.C
+++ b/src/regionFaModels/liquidFilm/subModels/kinematic/filmTurbulenceModel/filmTurbulenceModel/filmTurbulenceModel.C
@@ -106,7 +106,7 @@ tmp<areaScalarField> filmTurbulenceModel::Cw() const
             const scalarField& mu = film_.mu().primitiveField();
             const scalarField& h = film_.h().primitiveField();
             const scalar h0 = film_.h0().value();
-            Cw.primitiveFieldRef() = mu/((1.0/3.0)*(h + h0));
+            Cw.primitiveFieldRef() = 3*mu/(h + h0);
             Cw.min(5000.0);
             break;
         }
@@ -141,7 +141,7 @@ tmp<areaScalarField> filmTurbulenceModel::Cw() const
             const scalar h0 = film_.h0().value();
 
             Cw.primitiveFieldRef() =
-                sqr(n)*mag(g.value())*mag(Uf)/cbrt(h+ h0);
+                sqr(n)*mag(g.value())*mag(Uf)/cbrt(h+h0);
             break;
         }
         default:
diff --git a/src/regionFaModels/liquidFilm/subModels/kinematic/filmTurbulenceModel/filmTurbulenceModel/filmTurbulenceModel.H b/src/regionFaModels/liquidFilm/subModels/kinematic/filmTurbulenceModel/filmTurbulenceModel/filmTurbulenceModel.H
index ce2b620af9..5277983eec 100644
--- a/src/regionFaModels/liquidFilm/subModels/kinematic/filmTurbulenceModel/filmTurbulenceModel/filmTurbulenceModel.H
+++ b/src/regionFaModels/liquidFilm/subModels/kinematic/filmTurbulenceModel/filmTurbulenceModel/filmTurbulenceModel.H
@@ -35,8 +35,8 @@ SourceFiles
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef filmTurbulenceModel_H
-#define filmTurbulenceModel_H
+#ifndef regionFaModels_filmTurbulenceModel_H
+#define regionFaModels_filmTurbulenceModel_H
 
 #include "areaFieldsFwd.H"
 #include "runTimeSelectionTables.H"
diff --git a/src/regionFaModels/liquidFilm/subModels/kinematic/filmTurbulenceModel/laminar/laminar.C b/src/regionFaModels/liquidFilm/subModels/kinematic/filmTurbulenceModel/laminar/laminar.C
index 472f80dee2..7630e25466 100644
--- a/src/regionFaModels/liquidFilm/subModels/kinematic/filmTurbulenceModel/laminar/laminar.C
+++ b/src/regionFaModels/liquidFilm/subModels/kinematic/filmTurbulenceModel/laminar/laminar.C
@@ -96,13 +96,13 @@ tmp<faVectorMatrix> laminar::Su(areaVectorField& U) const
     tmp<areaVectorField> Up = film_.Up();
 
     // employ simple coeff-based model
-    areaScalarField Cs("Cs", Cf_*mag(Up() - U));
+    const dimensionedScalar Cf("Cf", dimVelocity, Cf_);
 
     tmp<areaScalarField> wf = Cw();
 
     return
     (
-       - fam::Sp(Cs, U) + Cs*Up()     // surface contribution
+       - fam::Sp(Cf, U) + Cf*Up()     // surface contribution
        - fam::Sp(wf(), U) + wf()*Uw() // wall contribution
     );
 }
diff --git a/src/regionFaModels/liquidFilm/subModels/kinematic/filmTurbulenceModel/laminar/laminar.H b/src/regionFaModels/liquidFilm/subModels/kinematic/filmTurbulenceModel/laminar/laminar.H
index a405929e23..b0e6b70541 100644
--- a/src/regionFaModels/liquidFilm/subModels/kinematic/filmTurbulenceModel/laminar/laminar.H
+++ b/src/regionFaModels/liquidFilm/subModels/kinematic/filmTurbulenceModel/laminar/laminar.H
@@ -34,10 +34,10 @@ SourceFiles
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef laminar_H
-#define laminar_H
+#ifndef regionFaModels_surfaceFilmModels_laminar_H
+#define regionFaModels_surfaceFilmModels_laminar_H
 
-#include "filmTurbulenceModel.H"
+#include "liquidFilm/subModels/kinematic/filmTurbulenceModel/filmTurbulenceModel/filmTurbulenceModel.H"
 #include "faCFD.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
diff --git a/src/regionFaModels/liquidFilm/subModels/kinematic/injectionModel/curvatureSeparation/curvatureSeparation.C b/src/regionFaModels/liquidFilm/subModels/kinematic/injectionModel/curvatureSeparation/curvatureSeparation.C
new file mode 100644
index 0000000000..ee40e48424
--- /dev/null
+++ b/src/regionFaModels/liquidFilm/subModels/kinematic/injectionModel/curvatureSeparation/curvatureSeparation.C
@@ -0,0 +1,324 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2011-2017 OpenFOAM Foundation
+    Copyright (C) 2020 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "curvatureSeparation.H"
+#include "addToRunTimeSelectionTable.H"
+#include "Time.H"
+#include "stringListOps.H"
+#include "cyclicPolyPatch.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace regionModels
+{
+namespace areaSurfaceFilmModels
+{
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+defineTypeNameAndDebug(curvatureSeparation, 0);
+addToRunTimeSelectionTable
+(
+    injectionModel,
+    curvatureSeparation,
+    dictionary
+);
+
+// * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
+
+tmp<areaScalarField> curvatureSeparation::calcInvR1
+(
+    const areaVectorField& U,
+    const scalarField& calcCosAngle
+) const
+{
+    const dimensionedScalar smallU(dimVelocity, ROOTVSMALL);
+    const areaVectorField UHat(U/(mag(U) + smallU));
+    tmp<areaScalarField> tinvR1
+    (
+        new areaScalarField("invR1", UHat & (UHat & -gradNHat_))
+    );
+
+    scalarField& invR1 = tinvR1.ref().primitiveFieldRef();
+
+    // apply defined patch radii
+    const scalar rMin = 1e-6;
+    const scalar definedInvR1 = 1.0/max(rMin, definedPatchRadii_);
+
+    if (definedPatchRadii_ > 0)
+    {
+        invR1 = definedInvR1;
+    }
+
+    // filter out large radii
+    const scalar rMax = 1e6;
+    forAll(invR1, i)
+    {
+        if ((mag(invR1[i]) < 1/rMax))
+        {
+            invR1[i] = -1.0;
+        }
+    }
+
+    return tinvR1;
+}
+
+
+tmp<scalarField> curvatureSeparation::calcCosAngle
+(
+    const edgeScalarField& phi
+) const
+{
+    const areaVectorField& U = film().Uf();
+    const dimensionedScalar smallU(dimVelocity, ROOTVSMALL);
+    const areaVectorField UHat(U/(mag(U) + smallU));
+
+    const faMesh& mesh = film().regionMesh();
+    const labelUList& own = mesh.edgeOwner();
+    const labelUList& nbr = mesh.edgeNeighbour();
+
+    scalarField phiMax(mesh.nFaces(), -GREAT);
+    scalarField cosAngle(UHat.size(), Zero);
+
+    const scalarField invR1(calcInvR1(U, cosAngle));
+
+    forAll(nbr, edgei)
+    {
+        const label cellO = own[edgei];
+        const label cellN = nbr[edgei];
+
+        if (phi[edgei] > phiMax[cellO])
+        {
+            phiMax[cellO] = phi[edgei];
+            cosAngle[cellO] = -gHat_ & UHat[cellN];
+        }
+        if (-phi[edgei] > phiMax[cellN])
+        {
+            phiMax[cellN] = -phi[edgei];
+            cosAngle[cellN] = -gHat_ & UHat[cellO];
+        }
+    }
+
+    cosAngle *= pos(invR1);
+
+    // checks
+    if (debug && mesh.time().writeTime())
+    {
+        areaScalarField volCosAngle
+        (
+            IOobject
+            (
+                "cosAngle",
+                film().primaryMesh().time().timeName(),
+                film().primaryMesh(),
+                IOobject::NO_READ
+            ),
+            film().regionMesh(),
+            dimensionedScalar(dimless, Zero)
+        );
+        volCosAngle.primitiveFieldRef() = cosAngle;
+        volCosAngle.correctBoundaryConditions();
+        volCosAngle.write();
+    }
+
+    return max(min(cosAngle, scalar(1)), scalar(-1));
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+curvatureSeparation::curvatureSeparation
+(
+    liquidFilmBase& film,
+    const dictionary& dict
+)
+:
+    injectionModel(type(), film, dict),
+    gradNHat_(fac::grad(film.regionMesh().faceAreaNormals())),
+    deltaByR1Min_(coeffDict_.getOrDefault<scalar>("deltaByR1Min", 0)),
+    definedPatchRadii_
+    (
+        coeffDict_.getOrDefault<scalar>("definedPatchRadii", 0)
+    ),
+    magG_(mag(film.g().value())),
+    gHat_(Zero),
+    fThreshold_
+    (
+        coeffDict_.getOrDefault<scalar>("fThreshold", 1e-8)
+    ),
+    minInvR1_
+    (
+        coeffDict_.getOrDefault<scalar>("minInvR1", 5)
+    )
+{
+    if (magG_ < ROOTVSMALL)
+    {
+        FatalErrorInFunction
+            << "Acceleration due to gravity must be non-zero"
+            << exit(FatalError);
+    }
+
+    gHat_ = film.g().value()/magG_;
+}
+
+
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+
+void curvatureSeparation::correct
+(
+    scalarField& availableMass,
+    scalarField& massToInject,
+    scalarField& diameterToInject
+)
+{
+    const faMesh& mesh = film().regionMesh();
+
+    const areaScalarField& delta = film().h();
+    const areaVectorField& U = film().Uf();
+    const edgeScalarField& phi = film().phi2s();
+    const areaScalarField& rho = film().rho();
+    const scalarField magSqrU(magSqr(film().Uf()));
+    const areaScalarField& sigma = film().sigma();
+
+    const scalarField cosAngle(calcCosAngle(phi));
+    const scalarField invR1(calcInvR1(U, cosAngle));
+
+
+    // calculate force balance
+    scalarField Fnet(mesh.nFaces(), Zero);
+    scalarField separated(mesh.nFaces(), Zero);
+
+    forAll(invR1, i)
+    {
+        if ((invR1[i] > minInvR1_) && (delta[i]*invR1[i] > deltaByR1Min_))
+        {
+            const scalar R1 = 1.0/(invR1[i] + ROOTVSMALL);
+            const scalar R2 = R1 + delta[i];
+
+            // inertial force
+            const scalar Fi = -delta[i]*rho[i]*magSqrU[i]*72.0/60.0*invR1[i];
+
+            // body force
+            const scalar Fb =
+               - 0.5*rho[i]*magG_*invR1[i]*(sqr(R1) - sqr(R2))*cosAngle[i];
+
+            // surface force
+            const 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;
+
+    addToInjectedMass(sum(massToInject));
+
+    if (debug && mesh.time().writeTime())
+    {
+        areaScalarField volFnet
+        (
+            IOobject
+            (
+                "Fnet",
+                film().primaryMesh().time().timeName(),
+                film().primaryMesh(),
+                IOobject::NO_READ
+            ),
+            mesh,
+            dimensionedScalar(dimForce, Zero)
+        );
+        volFnet.primitiveFieldRef() = Fnet;
+        volFnet.write();
+
+        areaScalarField areaSeparated
+        (
+            IOobject
+            (
+                "separated",
+                film().primaryMesh().time().timeName(),
+                film().primaryMesh(),
+                IOobject::NO_READ
+            ),
+            mesh,
+            dimensionedScalar(dimMass, Zero)
+        );
+        areaSeparated.primitiveFieldRef() = separated;
+        areaSeparated.write();
+
+        areaScalarField areaMassToInject
+        (
+            IOobject
+            (
+                "massToInject",
+                film().primaryMesh().time().timeName(),
+                film().primaryMesh(),
+                IOobject::NO_READ
+            ),
+            mesh,
+            dimensionedScalar(dimMass, Zero)
+        );
+        areaMassToInject.primitiveFieldRef() = massToInject;
+        areaMassToInject.write();
+
+        areaScalarField areaInvR1
+        (
+            IOobject
+            (
+                "InvR1",
+                film().primaryMesh().time().timeName(),
+                film().primaryMesh(),
+                IOobject::NO_READ
+            ),
+            mesh,
+            dimensionedScalar(inv(dimLength), Zero)
+        );
+        areaInvR1.primitiveFieldRef() = invR1;
+        areaInvR1.write();
+    }
+
+    injectionModel::correct();
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace areaSurfaceFilmModels
+} // End namespace regionModels
+} // End namespace Foam
+
+// ************************************************************************* //
diff --git a/src/regionFaModels/liquidFilm/subModels/kinematic/injectionModel/curvatureSeparation/curvatureSeparation.H b/src/regionFaModels/liquidFilm/subModels/kinematic/injectionModel/curvatureSeparation/curvatureSeparation.H
new file mode 100644
index 0000000000..0ec037e432
--- /dev/null
+++ b/src/regionFaModels/liquidFilm/subModels/kinematic/injectionModel/curvatureSeparation/curvatureSeparation.H
@@ -0,0 +1,167 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2011-2017 OpenFOAM Foundation
+    Copyright (C) 2021 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::regionModels::areaSurfaceFilmModels::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 "faCFD.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace regionModels
+{
+namespace areaSurfaceFilmModels
+{
+
+/*---------------------------------------------------------------------------*\
+                   Class curvatureSeparation Declaration
+\*---------------------------------------------------------------------------*/
+
+class curvatureSeparation
+:
+    public injectionModel
+{
+    // Private Member Functions
+
+        //- No copy construct
+        curvatureSeparation(const curvatureSeparation&) = delete;
+
+        //- No copy assignment
+        void operator=(const curvatureSeparation&) = delete;
+
+
+protected:
+
+    // Protected Data
+
+        //- Gradient of surface normals
+        areaTensorField 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
+        scalar definedPatchRadii_;
+
+        //- Magnitude of gravity vector
+        scalar magG_;
+
+        //- Direction of gravity vector
+        vector gHat_;
+
+        //- Threshold force for separation
+        scalar fThreshold_;
+
+        //- Minimum inv R1 for separation
+        scalar minInvR1_;
+
+
+    // Protected Member Functions
+
+        //- Calculate local (inverse) radius of curvature
+        tmp<areaScalarField> calcInvR1
+        (
+            const areaVectorField& U,
+            const scalarField& calcCosAngle
+        ) const;
+
+        //- Calculate the cosine of the angle between gravity vector and
+        //- cell out flow direction
+        tmp<scalarField> calcCosAngle(const edgeScalarField& phi) const;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("curvatureSeparation");
+
+
+    // Constructors
+
+        //- Construct from surface film model
+        curvatureSeparation
+        (
+            liquidFilmBase& film,
+            const dictionary& dict
+        );
+
+
+    //- Destructor
+    virtual ~curvatureSeparation() = default;
+
+
+    // Member Functions
+
+        // Evolution
+
+            //- Correct
+            virtual void correct
+            (
+                scalarField& availableMass,
+                scalarField& massToInject,
+                scalarField& diameterToInject
+            );
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace areaSurfaceFilmModels
+} // End namespace regionModels
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/regionFaModels/liquidFilm/subModels/kinematic/injectionModel/injectionModel/injectionModel.C b/src/regionFaModels/liquidFilm/subModels/kinematic/injectionModel/injectionModel/injectionModel.C
new file mode 100644
index 0000000000..b1e5475a85
--- /dev/null
+++ b/src/regionFaModels/liquidFilm/subModels/kinematic/injectionModel/injectionModel/injectionModel.C
@@ -0,0 +1,106 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2021 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "injectionModel.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace regionModels
+{
+namespace areaSurfaceFilmModels
+{
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+defineTypeNameAndDebug(injectionModel, 0);
+defineRunTimeSelectionTable(injectionModel, dictionary);
+
+// * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
+
+void injectionModel::addToInjectedMass(const scalar dMass)
+{
+    injectedMass_ += dMass;
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+injectionModel::injectionModel(liquidFilmBase& film)
+:
+    filmSubModelBase(film),
+    injectedMass_(0.0)
+{}
+
+
+injectionModel::injectionModel
+(
+    const word& modelType,
+    liquidFilmBase& film,
+    const dictionary& dict
+)
+:
+    filmSubModelBase(film, dict, typeName, modelType),
+    injectedMass_(0.0)
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+injectionModel::~injectionModel()
+{}
+
+
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+
+void injectionModel::correct()
+{
+    if (writeTime())
+    {
+        scalar injectedMass0 = getModelProperty<scalar>("injectedMass");
+        injectedMass0 += returnReduce(injectedMass_, sumOp<scalar>());
+        setModelProperty<scalar>("injectedMass", injectedMass0);
+        injectedMass_ = 0.0;
+    }
+}
+
+
+scalar injectionModel::injectedMassTotal() const
+{
+    const scalar injectedMass0 = getModelProperty<scalar>("injectedMass");
+    return injectedMass0 + returnReduce(injectedMass_, sumOp<scalar>());
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace areaSurfaceFilmModels
+} // End namespace regionModels
+} // End namespace Foam
+
+// ************************************************************************* //
diff --git a/src/regionFaModels/liquidFilm/subModels/kinematic/injectionModel/injectionModel/injectionModel.H b/src/regionFaModels/liquidFilm/subModels/kinematic/injectionModel/injectionModel/injectionModel.H
new file mode 100644
index 0000000000..5bf117f7a9
--- /dev/null
+++ b/src/regionFaModels/liquidFilm/subModels/kinematic/injectionModel/injectionModel/injectionModel.H
@@ -0,0 +1,169 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2021 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::regionModels::areaSurfaceFilmModels::injectionModel
+
+Description
+    Base class for film injection models, handling mass transfer from the
+    film.
+
+SourceFiles
+    injectionModel.C
+    injectionModelNew.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef injectionModel_H
+#define injectionModel_H
+
+#include "filmSubModelBase.H"
+#include "runTimeSelectionTables.H"
+#include "scalarField.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace regionModels
+{
+namespace areaSurfaceFilmModels
+{
+
+/*---------------------------------------------------------------------------*\
+                      Class injectionModel Declaration
+\*---------------------------------------------------------------------------*/
+
+class injectionModel
+:
+    public filmSubModelBase
+{
+    // Private Data
+
+        //- Injected mass
+        scalar injectedMass_;
+
+
+    // Private Member Functions
+
+        //- No copy construct
+        injectionModel(const injectionModel&) = delete;
+
+        //- No copy assignment
+        void operator=(const injectionModel&) = delete;
+
+
+protected:
+
+    // Protected Member Functions
+
+        //- Add to injected mass
+        void addToInjectedMass(const scalar dMass);
+
+        //- Correct
+        void correct();
+
+
+public:
+
+    //- Runtime type information
+    TypeName("injectionModel");
+
+
+    // Declare runtime constructor selection table
+
+         declareRunTimeSelectionTable
+         (
+             autoPtr,
+             injectionModel,
+             dictionary,
+             (
+                liquidFilmBase& film,
+                const dictionary& dict
+             ),
+             (film, dict)
+         );
+
+
+    // Constructors
+
+        //- Construct null
+        injectionModel(liquidFilmBase& film);
+
+        //- Construct from type name, dictionary and surface film model
+        injectionModel
+        (
+            const word& modelType,
+            liquidFilmBase& film,
+            const dictionary& dict
+        );
+
+
+    // Selectors
+
+        //- Return a reference to the selected injection model
+        static autoPtr<injectionModel> New
+        (
+            liquidFilmBase& film,
+            const dictionary& dict,
+            const word& mdoelType
+        );
+
+
+    //- Destructor
+    virtual ~injectionModel();
+
+
+    // Member Functions
+
+        //- Correct
+        virtual void correct
+        (
+            scalarField& availableMass,
+            scalarField& massToInject,
+            scalarField& diameterToInject
+        ) = 0;
+
+        //- Return the total mass injected
+        virtual scalar injectedMassTotal() const;
+
+        //- Accumulate the total mass injected for the patches into the
+        //- scalarField provided
+        virtual void patchInjectedMassTotals(scalar& patchMasses) const
+        {}
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace areaSurfaceFilmModels
+} // End namespace regionModels
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/regionFaModels/liquidFilm/subModels/kinematic/injectionModel/injectionModel/injectionModelNew.C b/src/regionFaModels/liquidFilm/subModels/kinematic/injectionModel/injectionModel/injectionModelNew.C
new file mode 100644
index 0000000000..e449744693
--- /dev/null
+++ b/src/regionFaModels/liquidFilm/subModels/kinematic/injectionModel/injectionModel/injectionModelNew.C
@@ -0,0 +1,73 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2021 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "injectionModel.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace regionModels
+{
+namespace areaSurfaceFilmModels
+{
+
+// * * * * * * * * * * * * * * * * Selectors * * * * * * * * * * * * * * * * //
+
+autoPtr<injectionModel> injectionModel::New
+(
+    liquidFilmBase& model,
+    const dictionary& dict,
+    const word& modelType
+)
+{
+    Info<< "        " << modelType << endl;
+
+    auto cstrIter = dictionaryConstructorTablePtr_->cfind(modelType);
+
+    if (!cstrIter.found())
+    {
+        FatalIOErrorInLookup
+        (
+            dict,
+            "injectionModel",
+            modelType,
+            *dictionaryConstructorTablePtr_
+        ) << exit(FatalIOError);
+    }
+
+    return autoPtr<injectionModel>(cstrIter()(model, dict));
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace areaSurfaceFilmModels
+} // End namespace regionModels
+} // End namespace Foam
+
+// ************************************************************************* //
diff --git a/src/regionFaModels/liquidFilm/subModels/kinematic/injectionModel/injectionModelList/injectionModelList.C b/src/regionFaModels/liquidFilm/subModels/kinematic/injectionModel/injectionModelList/injectionModelList.C
new file mode 100644
index 0000000000..01900da387
--- /dev/null
+++ b/src/regionFaModels/liquidFilm/subModels/kinematic/injectionModel/injectionModelList/injectionModelList.C
@@ -0,0 +1,169 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2020 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "injectionModelList.H"
+#include "volFields.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace regionModels
+{
+namespace areaSurfaceFilmModels
+{
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+injectionModelList::injectionModelList(liquidFilmBase& film)
+:
+    PtrList<injectionModel>(),
+    filmSubModelBase(film)
+{}
+
+
+injectionModelList::injectionModelList
+(
+    liquidFilmBase& film,
+    const dictionary& dict
+)
+:
+    PtrList<injectionModel>(),
+    filmSubModelBase
+    (
+        "injectionModelList",
+        film,
+        dict,
+        "injectionModelList",
+        "injectionModelList"
+    ),
+    massInjected_(Zero)
+{
+    const wordList activeModels(dict.lookup("injectionModels"));
+
+    wordHashSet models(activeModels);
+
+    Info<< "    Selecting film injection models" << endl;
+    if (models.size())
+    {
+        this->setSize(models.size());
+
+        label i = 0;
+        for (const word& model : models)
+        {
+            set(i, injectionModel::New(film, dict, model));
+            i++;
+        }
+    }
+    else
+    {
+        Info<< "        none" << endl;
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+injectionModelList::~injectionModelList()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void injectionModelList::correct
+(
+    scalarField& availableMass,
+    volScalarField& massToInject,
+    volScalarField& diameterToInject
+)
+{
+    const label patchi = film().patchID();
+
+    // Correct models that accumulate mass and diameter transfers
+    forAll(*this, i)
+    {
+        injectionModel& im = operator[](i);
+        im.correct
+        (
+            availableMass,
+            massToInject.boundaryFieldRef()[patchi],
+            diameterToInject.boundaryFieldRef()[patchi]
+        );
+    }
+
+    massInjected_ += gSum(massToInject.boundaryField()[patchi]);
+}
+
+
+void injectionModelList::info(Ostream& os)
+{
+    const polyBoundaryMesh& pbm = film().primaryMesh().boundaryMesh();
+
+    scalar injectedMass = 0;
+    scalar patchInjectedMasses = 0;
+
+    forAll(*this, i)
+    {
+        const injectionModel& im = operator[](i);
+        injectedMass += im.injectedMassTotal();
+        im.patchInjectedMassTotals(patchInjectedMasses);
+    }
+
+    os  << indent << "injected mass      = " << injectedMass << nl;
+
+    const label patchi = film().patchID();
+
+    if (mag(patchInjectedMasses) > VSMALL)
+    {
+        os  << indent << indent << "from patch " << pbm[patchi].name()
+            << " = " << patchInjectedMasses << nl;
+    }
+
+    scalar mass0(Zero);
+    this->getBaseProperty("massInjected", mass0);
+
+    scalar mass(massInjected_);
+    mass += mass0;
+
+    Info<< indent << "  - patch: " << pbm[patchi].name() << "  "
+        << mass << endl;
+
+    if (film().primaryMesh().time().writeTime())
+    {
+        setBaseProperty("massInjected", mass);
+        massInjected_ = 0.0;
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace areaSurfaceFilmModels
+} // End namespace regionModels
+} // End namespace Foam
+
+// ************************************************************************* //
diff --git a/src/regionFaModels/liquidFilm/subModels/kinematic/injectionModel/injectionModelList/injectionModelList.H b/src/regionFaModels/liquidFilm/subModels/kinematic/injectionModel/injectionModelList/injectionModelList.H
new file mode 100644
index 0000000000..b159f9486a
--- /dev/null
+++ b/src/regionFaModels/liquidFilm/subModels/kinematic/injectionModel/injectionModelList/injectionModelList.H
@@ -0,0 +1,127 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2020 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::regionModels::surfaceFilmModels::injectionModelList
+
+Description
+    List container for film injection models
+
+SourceFiles
+    injectionModelList.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef injectionModelList_H
+#define injectionModelList_H
+
+#include "PtrList.H"
+#include "injectionModel.H"
+#include "filmSubModelBase.H"
+#include "scalarField.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace regionModels
+{
+namespace areaSurfaceFilmModels
+{
+
+/*---------------------------------------------------------------------------*\
+                    Class injectionModelList Declaration
+\*---------------------------------------------------------------------------*/
+
+class injectionModelList
+:
+    public PtrList<injectionModel>,
+    public filmSubModelBase
+{
+    // Private Data
+
+        //- List of mass injected per patch
+        scalar massInjected_;
+
+
+    // Private Member Functions
+
+        //- No copy construct
+        injectionModelList(const injectionModelList&) = delete;
+
+        //- No copy assignment
+        void operator=(const injectionModelList&) = delete;
+
+
+public:
+
+    // Constructors
+
+        //- Construct null
+        injectionModelList(liquidFilmBase& film);
+
+        //- Construct from type name, dictionary and surface film model
+        injectionModelList
+        (
+            liquidFilmBase& film,
+            const dictionary& dict
+        );
+
+
+    //- Destructor
+    virtual ~injectionModelList();
+
+
+    // Member Functions
+
+        // Evolution
+
+            //- Correct
+            virtual void correct
+            (
+                scalarField& availableMass,
+                volScalarField& massToInject,
+                volScalarField& diameterToInject
+            );
+
+
+        // I-O
+
+            //- Provide some info
+            virtual void info(Ostream& os);
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace areaSurfaceFilmModels
+} // End namespace regionModels
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //