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

src/finiteVolume/Make/files

@@ -204,8 +204,8 @@ $(interpolation)/interpolationCellPoint/makeInterpolationCellPoint.C
 $(interpolation)/interpolationCellPointFace/makeInterpolationCellPointFace.C
 $(interpolation)/interpolationCellPointWallModified/cellPointWeightWallModified/cellPointWeightWallModified.C
 $(interpolation)/interpolationCellPointWallModified/makeInterpolationCellPointWallModified.C
-$(interpolation)/interpolationPoint/pointMVCWeight.C
+$(interpolation)/interpolationPointMVC/pointMVCWeight.C
-$(interpolation)/interpolationPoint/makeInterpolationPoint.C
+$(interpolation)/interpolationPointMVC/makeInterpolationPointMVC.C
 
 volPointInterpolation = interpolation/volPointInterpolation
 $(volPointInterpolation)/volPointInterpolation.C

src/finiteVolume/interpolation/interpolation/interpolationPointMVC/interpolationPointMVC.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2004-2010 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2010-2011 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -23,13 +23,13 @@ License
 
 \*---------------------------------------------------------------------------*/
 
-#include "interpolationPoint.H"
+#include "interpolationPointMVC.H"
 #include "volPointInterpolation.H"
 
 // * * * * * * * * * * * * * * * * Constructor * * * * * * * * * * * * * * * //
 
 template<class Type>
-Foam::interpolationPoint<Type>::interpolationPoint
+Foam::interpolationPointMVC<Type>::interpolationPointMVC
 (
     const GeometricField<Type, fvPatchField, volMesh>& psi
 )

src/finiteVolume/interpolation/interpolation/interpolationPointMVC/interpolationPointMVC.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2004-2010 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2010-2011 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -22,7 +22,7 @@ License
     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
 
 Class
-    Foam::interpolationPoint
+    Foam::interpolationPointMVC
 
 Description
     Given cell centre values interpolates to vertices and uses these to
@@ -30,8 +30,8 @@ Description
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef interpolationPoint_H
+#ifndef interpolationPointMVC_H
-#define interpolationPoint_H
+#define interpolationPointMVC_H
 
 #include "interpolation.H"
 #include "pointMVCWeight.H"
@@ -42,11 +42,11 @@ namespace Foam
 {
 
 /*---------------------------------------------------------------------------*\
-                   Class interpolationPoint Declaration
+                   Class interpolationPointMVC Declaration
 \*---------------------------------------------------------------------------*/
 
 template<class Type>
-class interpolationPoint
+class interpolationPointMVC
 :
     public interpolation<Type>
 {
@@ -57,16 +57,17 @@ protected:
         //- Interpolated volfield
         const GeometricField<Type, pointPatchField, pointMesh> psip_;
 
+
 public:
 
     //- Runtime type information
-    TypeName("point");
+    TypeName("pointMVC");
 
 
     // Constructors
 
         //- Construct from components
-        interpolationPoint
+        interpolationPointMVC
         (
             const GeometricField<Type, fvPatchField, volMesh>& psi
         );
@@ -93,12 +94,12 @@ public:
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-#include "interpolationPointI.H"
+#include "interpolationPointMVCI.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 #ifdef NoRepository
-#   include "interpolationPoint.C"
+#   include "interpolationPointMVC.C"
 #endif
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

src/finiteVolume/interpolation/interpolation/interpolationPointMVC/interpolationPointMVCI.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2004-2010 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2010-2011 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -26,7 +26,7 @@ License
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 template<class Type>
-inline Type Foam::interpolationPoint<Type>::interpolate
+inline Type Foam::interpolationPointMVC<Type>::interpolate
 (
     const pointMVCWeight& cpw
 ) const
@@ -36,7 +36,7 @@ inline Type Foam::interpolationPoint<Type>::interpolate
 
 
 template<class Type>
-inline Type Foam::interpolationPoint<Type>::interpolate
+inline Type Foam::interpolationPointMVC<Type>::interpolate
 (
     const vector& position,
     const label cellI,

src/finiteVolume/interpolation/interpolation/interpolationPointMVC/makeInterpolationPointMVC.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2004-2010 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2010-2011 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -23,13 +23,13 @@ License
 
 \*---------------------------------------------------------------------------*/
 
-#include "interpolationPoint.H"
+#include "interpolationPointMVC.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 namespace Foam
 {
-    makeInterpolation(interpolationPoint);
+    makeInterpolation(interpolationPointMVC);
 }
 
 // ************************************************************************* //

src/finiteVolume/interpolation/interpolation/interpolationPointMVC/pointMVCWeight.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2004-2010 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2010-2011 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -56,7 +56,7 @@ void Foam::pointMVCWeight::calcWeights
     forAll(f, j)
     {
         label jPlus1 = f.fcIndex(j);
-        scalar l = mag(u[j]-u[jPlus1]);
+        scalar l = mag(u[j] - u[jPlus1]);
         theta[j] = 2.0*Foam::asin(l/2.0);
     }
 
@@ -68,7 +68,7 @@ void Foam::pointMVCWeight::calcWeights
         weights[pid] =
             1.0
           / dist[pid]
-          * (Foam::tan(theta[jMin1]/2.0)+Foam::tan(theta[j]/2.0));
+          * (Foam::tan(theta[jMin1]/2.0) + Foam::tan(theta[j]/2.0));
         sumWeight += weights[pid];
     }
 
@@ -132,11 +132,10 @@ void Foam::pointMVCWeight::calcWeights
             //Pout<< "    uj:" << u[j] << " ujPlus1:" << u[jPlus1]
             //    << " temp:" << temp << endl;
 
-            scalar l = mag(u[j]-u[jPlus1]);
+            scalar l = min(mag(u[j] - u[jPlus1]), 2.0);
             scalar angle = 2.0*Foam::asin(l/2.0);
 
-            //Pout<< "    j:" << j << " l:" << l
+            //Pout<< "    j:" << j << " l:" << l << " angle:" << angle << endl;
-            //    << " angle:" << angle << endl;
 
             v += 0.5*angle*temp;
         }
@@ -145,14 +144,14 @@ void Foam::pointMVCWeight::calcWeights
         v /= vNorm;
 
         // Make sure v points towards the polygon
-//if (((v&u[0]) < 0) != (mesh.faceOwner()[faceI] != cellIndex_))
+        //if (((v&u[0]) < 0) != (mesh.faceOwner()[faceI] != cellIndex_))
-//{
+        //{
-//    FatalErrorIn("pointMVCWeight::calcWeights(..)")
+        //    FatalErrorIn("pointMVCWeight::calcWeights(..)")
-//        << "v:" << v << " u[0]:" << u[0]
+        //        << "v:" << v << " u[0]:" << u[0]
-//        << exit(FatalError);
+        //        << exit(FatalError);
-//}
+        //}
 
-        if ((v&u[0]) < 0)
+        if ((v & u[0]) < 0)
         {
             v = -v;
         }
@@ -165,22 +164,22 @@ void Foam::pointMVCWeight::calcWeights
             label jPlus1 = f.fcIndex(j);
             //Pout<< "    uj:" << u[j] << " ujPlus1:" << u[jPlus1] << endl;
 
-            vector n0 = u[j] ^ v;
+            vector n0 = u[j]^v;
             n0 /= mag(n0);
-            vector n1 = u[jPlus1] ^ v;
+            vector n1 = u[jPlus1]^v;
             n1 /= mag(n1);
 
-            scalar l = mag(n0-n1);
+            scalar l = min(mag(n0 - n1), 2.0);
             //Pout<< "    l:" << l << endl;
             alpha(j) = 2.0*Foam::asin(l/2.0);
 
-            vector temp = n0 ^ n1;
+            vector temp = n0^n1;
             if ((temp&v) < 0.0)
             {
                 alpha[j] = -alpha[j];
             }
 
-            l = mag(u[j]-v);
+            l = min(mag(u[j] - v), 2.0);
             //Pout<< "    l:" << l << endl;
             theta(j) = 2.0*Foam::asin(l/2.0);
         }
@@ -211,7 +210,7 @@ void Foam::pointMVCWeight::calcWeights
             sum +=
                 1.0
               / Foam::tan(theta[j])
-              * (Foam::tan(alpha[j]/2.0)+Foam::tan(alpha[jMin1]/2.0));
+              * (Foam::tan(alpha[j]/2.0) + Foam::tan(alpha[jMin1]/2.0));
         }
 
         // The special case when x lies on the polygon, handle it using 2D mvc.
@@ -234,7 +233,7 @@ void Foam::pointMVCWeight::calcWeights
               / sum
               / dist[pid]
               / Foam::sin(theta[j])
-              * (Foam::tan(alpha[j]/2.0)+Foam::tan(alpha[jMin1]/2.0));
+              * (Foam::tan(alpha[j]/2.0) + Foam::tan(alpha[jMin1]/2.0));
         }
     }
 

src/finiteVolume/interpolation/interpolation/interpolationPointMVC/pointMVCWeight.H

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

src/finiteVolume/interpolation/interpolation/interpolationPointMVC/pointMVCWeightI.H

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

src/lagrangian/intermediate/submodels/CloudFunctionObjects/CloudFunctionObjectList/CloudFunctionObjectList.C

@@ -89,12 +89,12 @@ Foam::CloudFunctionObjectList<CloudType>::CloudFunctionObjectList
 template<class CloudType>
 Foam::CloudFunctionObjectList<CloudType>::CloudFunctionObjectList
 (
-    const CloudFunctionObjectList& ppm
+    const CloudFunctionObjectList& cfol
 )
 :
-    PtrList<CloudFunctionObject<CloudType> >(ppm),
+    PtrList<CloudFunctionObject<CloudType> >(cfol),
-    owner_(ppm.owner_),
+    owner_(cfol.owner_),
-    dict_(ppm.dict_)
+    dict_(cfol.dict_)
 {}
 
 

src/meshTools/indexedOctree/treeDataPrimitivePatch.C

@@ -307,7 +307,7 @@ getVolumeType
 
             forAll(eFaces, i)
             {
-                edgeNormal += patch_.faceNormal()[eFaces[i]];
+                edgeNormal += patch_.faceNormals()[eFaces[i]];
             }
 
             if (debug & 2)

src/regionModels/surfaceFilmModels/Make/files

@@ -10,6 +10,13 @@ thermoSingleLayer/thermoSingleLayer.C
 submodels/subModelBase.C
 
 KINEMATICMODELS=submodels/kinematic
+$(KINEMATICMODELS)/force/force/force.C
+$(KINEMATICMODELS)/force/force/forceNew.C
+$(KINEMATICMODELS)/force/forceList/forceList.C
+$(KINEMATICMODELS)/force/contactAngleForce/contactAngleForce.C
+$(KINEMATICMODELS)/force/surfaceShearForce/surfaceShearForce.C
+$(KINEMATICMODELS)/force/thermocapillaryForce/thermocapillaryForce.C
+
 $(KINEMATICMODELS)/injectionModel/injectionModel/injectionModel.C
 $(KINEMATICMODELS)/injectionModel/injectionModel/injectionModelNew.C
 $(KINEMATICMODELS)/injectionModel/injectionModelList/injectionModelList.C

src/regionModels/surfaceFilmModels/kinematicSingleLayer/kinematicSingleLayer.C

@@ -34,6 +34,10 @@ License
 #include "directMappedWallPolyPatch.H"
 #include "mapDistribute.H"
 
+#include "cachedRandom.H"
+#include "normal.H"
+#include "mathematicalConstants.H"
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 namespace Foam
@@ -61,8 +65,6 @@ bool kinematicSingleLayer::read()
         solution.lookup("nCorr") >> nCorr_;
         solution.lookup("nNonOrthCorr") >> nNonOrthCorr_;
 
-        coeffs_.lookup("Cf") >> Cf_;
-
         return true;
     }
     else
@@ -76,9 +78,10 @@ void kinematicSingleLayer::correctThermoFields()
 {
     if (thermoModel_ == tmConstant)
     {
-        rho_ == dimensionedScalar(coeffs_.lookup("rho0"));
+        const dictionary& constDict(coeffs_.subDict("constantThermoCoeffs"));
-        mu_ == dimensionedScalar(coeffs_.lookup("mu0"));
+        rho_ == dimensionedScalar(constDict.lookup("rho0"));
-        sigma_ == dimensionedScalar(coeffs_.lookup("sigma0"));
+        mu_ == dimensionedScalar(constDict.lookup("mu0"));
+        sigma_ == dimensionedScalar(constDict.lookup("sigma0"));
     }
     else
     {
@@ -273,25 +276,6 @@ void kinematicSingleLayer::updateSurfaceVelocities()
 }
 
 
-tmp<fvVectorMatrix> kinematicSingleLayer::tau(volVectorField& U) const
-{
-    // Calculate shear stress
-    volScalarField Cs("Cs", rho_*Cf_*mag(Us_ - U));
-    volScalarField Cw
-    (
-        "Cw",
-        mu_/(0.3333*(delta_ + dimensionedScalar("SMALL", dimLength, SMALL)))
-    );
-    Cw.min(1.0e+06);
-
-    return
-    (
-       - fvm::Sp(Cs, U) + Cs*Us_ // surface contribution
-       - fvm::Sp(Cw, U) + Cw*Uw_ // wall contribution
-    );
-}
-
-
 tmp<Foam::fvVectorMatrix> kinematicSingleLayer::solveMomentum
 (
     const volScalarField& pu,
@@ -312,9 +296,8 @@ tmp<Foam::fvVectorMatrix> kinematicSingleLayer::solveMomentum
       + fvm::div(phi_, U_)
      ==
       - USp_
-      + tau(U_)
-      + fvc::grad(sigma_)
       - fvm::SuSp(rhoSp_, U_)
+      + forces_.correct(U_)
     );
 
     fvVectorMatrix& UEqn = tUEqn();
@@ -459,8 +442,6 @@ kinematicSingleLayer::kinematicSingleLayer
 
     cumulativeContErr_(0.0),
 
-    Cf_(readScalar(coeffs().lookup("Cf"))),
-
     rho_
     (
         IOobject
@@ -773,6 +754,8 @@ kinematicSingleLayer::kinematicSingleLayer
 
     injection_(*this, coeffs_),
 
+    forces_(*this, coeffs_),
+
     addedMassTotal_(0.0)
 {
     if (readFields)

src/regionModels/surfaceFilmModels/kinematicSingleLayer/kinematicSingleLayer.H

@@ -42,6 +42,7 @@ SourceFiles
 #include "fvMatrices.H"
 
 #include "injectionModelList.H"
+#include "forceList.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -93,12 +94,6 @@ protected:
             scalar cumulativeContErr_;
 
 
-        // Model parameters
-
-            //- Skin frition coefficient for film/primary region interface
-            scalar Cf_;
-
-
         // Thermo properties
 
             // Fields
@@ -199,6 +194,9 @@ protected:
             //- Cloud injection
             injectionModelList injection_;
 
+            //- List of film forces
+            forceList forces_;
+
 
        // Checks
 
@@ -238,9 +236,6 @@ protected:
         //- Update film surface velocities
         virtual void updateSurfaceVelocities();
 
-        //- Return the stress term for the momentum equation
-        virtual tmp<fvVectorMatrix> tau(volVectorField& dU) const;
-
         //- Constrain a film region master/slave boundaries of a field to a
         //  given value
         template<class Type>
@@ -314,12 +309,6 @@ public:
             inline label nNonOrthCorr() const;
 
 
-        // Model parameters
-
-            //- Return the skin friction coefficient
-            inline scalar Cf() const;
-
-
         // Thermo properties
 
             //- Return const access to the dynamic viscosity / [Pa.s]

src/regionModels/surfaceFilmModels/kinematicSingleLayer/kinematicSingleLayerI.H

@@ -61,12 +61,6 @@ inline label kinematicSingleLayer::nNonOrthCorr() const
 }
 
 
-inline scalar kinematicSingleLayer::Cf() const
-{
-    return Cf_;
-}
-
-
 inline const volScalarField& kinematicSingleLayer::mu() const
 {
     return mu_;

src/regionModels/surfaceFilmModels/submodels/kinematic/force/contactAngleForce/contactAngleForce.C

@@ -0,0 +1,164 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2011 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 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 "contactAngleForce.H"
+#include "addToRunTimeSelectionTable.H"
+#include "fvcGrad.H"
+#include "unitConversion.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace regionModels
+{
+namespace surfaceFilmModels
+{
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+defineTypeNameAndDebug(contactAngleForce, 0);
+addToRunTimeSelectionTable(force, contactAngleForce, dictionary);
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+contactAngleForce::contactAngleForce
+(
+    const surfaceFilmModel& owner,
+    const dictionary& dict
+)
+:
+    force(typeName, owner, dict),
+    deltaWet_(readScalar(coeffs_.lookup("deltaWet"))),
+    Ccf_(readScalar(coeffs_.lookup("Ccf"))),
+    rndGen_(label(0), -1),
+    distribution_
+    (
+        distributionModels::distributionModel::New
+        (
+            coeffs_.subDict("contactAngleDistribution"),
+            rndGen_
+        )
+    )
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+contactAngleForce::~contactAngleForce()
+{}
+
+
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+
+tmp<fvVectorMatrix> contactAngleForce::correct(volVectorField& U)
+{
+    tmp<volVectorField> tForce
+    (
+        new volVectorField
+        (
+            IOobject
+            (
+                "contactForce",
+                owner_.time().timeName(),
+                owner_.regionMesh(),
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            owner_.regionMesh(),
+            dimensionedVector("zero", dimForce/dimArea, vector::zero)
+        )
+    );
+
+    vectorField& force = tForce().internalField();
+
+    const labelUList& own = owner_.regionMesh().owner();
+    const labelUList& nbr = owner_.regionMesh().neighbour();
+
+    const scalarField& magSf = owner_.magSf();
+
+    const volScalarField& delta = owner_.delta();
+    const volScalarField& sigma = owner_.sigma();
+
+    volScalarField alpha
+    (
+        "alpha",
+        pos(delta - dimensionedScalar("deltaWet", dimLength, deltaWet_))
+    );
+    volVectorField gradAlpha(fvc::grad(alpha));
+
+    scalarField nHits(force.size(), 0.0);
+
+    forAll(nbr, faceI)
+    {
+        const label cellO = own[faceI];
+        const label cellN = nbr[faceI];
+
+        label cellI = -1;
+        if ((delta[cellO] > deltaWet_) && (delta[cellN] < deltaWet_))
+        {
+            cellI = cellO;
+        }
+        else if ((delta[cellO] < deltaWet_) && (delta[cellN] > deltaWet_))
+        {
+            cellI = cellN;
+        }
+
+        if (cellI != -1)
+        {
+//            const scalar dx = Foam::sqrt(magSf[cellI]);
+            const scalar dx = owner_.regionMesh().deltaCoeffs()[faceI];
+            const vector n =
+                gradAlpha[cellI]/(mag(gradAlpha[cellI]) + ROOTVSMALL);
+            scalar theta = cos(degToRad(distribution_->sample()));
+            force[cellI] += Ccf_*n*sigma[cellI]*(1.0 - theta)/dx;
+            nHits[cellI]++;
+        }
+    }
+
+    nHits = max(nHits, 1.0);
+    force /= (nHits*magSf);
+
+    if (owner_.regionMesh().time().outputTime())
+    {
+        tForce().write();
+    }
+
+    tmp<fvVectorMatrix>
+        tfvm(new fvVectorMatrix(U, dimForce/dimArea*dimVolume));
+
+    tfvm() += tForce;
+
+    return tfvm;
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace surfaceFilmModels
+} // End namespace regionModels
+} // End namespace Foam
+
+// ************************************************************************* //

src/regionModels/surfaceFilmModels/submodels/kinematic/force/contactAngleForce/contactAngleForce.H

@@ -0,0 +1,125 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2011 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 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::contactAngleForce
+
+Description
+    Film contact angle force
+
+SourceFiles
+    contactAngleForce.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef contactAngleForce_H
+#define contactAngleForce_H
+
+#include "force.H"
+#include "distributionModel.H"
+#include "cachedRandom.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace regionModels
+{
+namespace surfaceFilmModels
+{
+
+/*---------------------------------------------------------------------------*\
+                      Class contactAngleForce Declaration
+\*---------------------------------------------------------------------------*/
+
+class contactAngleForce
+:
+    public force
+{
+private:
+
+    // Private Data
+
+        //- Threshold film thickness beyon which the film is 'wet'
+        scalar deltaWet_;
+
+        //- Coefficient applied to the contact angle force
+        scalar Ccf_;
+
+        //- Random number generator
+        cachedRandom rndGen_;
+
+        //- Parcel size PDF model
+        const autoPtr<distributionModels::distributionModel> distribution_;
+
+
+
+    // Private member functions
+
+        //- Disallow default bitwise copy construct
+        contactAngleForce(const contactAngleForce&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const contactAngleForce&);
+
+
+public:
+
+    //- Runtime type information
+    TypeName("contactAngle");
+
+
+    // Constructors
+
+        //- Construct from surface film model
+        contactAngleForce
+        (
+            const surfaceFilmModel& owner,
+            const dictionary& dict
+        );
+
+
+    //- Destructor
+    virtual ~contactAngleForce();
+
+
+    // Member Functions
+
+        // Evolution
+
+            //- Correct
+            virtual tmp<fvVectorMatrix> correct(volVectorField& U);
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace surfaceFilmModels
+} // End namespace regionModels
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

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

@@ -0,0 +1,73 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2011 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 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 "force.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace regionModels
+{
+namespace surfaceFilmModels
+{
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+defineTypeNameAndDebug(force, 0);
+defineRunTimeSelectionTable(force, dictionary);
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+force::force(const surfaceFilmModel& owner)
+:
+    subModelBase(owner)
+{}
+
+
+force::force
+(
+    const word& type,
+    const surfaceFilmModel& owner,
+    const dictionary& dict
+)
+:
+    subModelBase(type, owner, dict)
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+force::~force()
+{}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace surfaceFilmModels
+} // End namespace regionModels
+} // End namespace Foam
+
+// ************************************************************************* //

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

@@ -0,0 +1,138 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2011 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 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::force
+
+Description
+    Base class for film (stress-based) force models
+
+SourceFiles
+    force.C
+    forceNew.C
+\*---------------------------------------------------------------------------*/
+
+#ifndef force_H
+#define force_H
+
+#include "subModelBase.H"
+#include "runTimeSelectionTables.H"
+#include "fvMatrices.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace regionModels
+{
+namespace surfaceFilmModels
+{
+
+/*---------------------------------------------------------------------------*\
+                          Class force Declaration
+\*---------------------------------------------------------------------------*/
+
+class force
+:
+    public subModelBase
+{
+private:
+
+    // Private Member Functions
+
+        //- Disallow default bitwise copy construct
+        force(const force&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const force&);
+
+
+public:
+
+    //- Runtime type information
+    TypeName("force");
+
+
+    // Declare runtime constructor selection table
+
+         declareRunTimeSelectionTable
+         (
+             autoPtr,
+             force,
+             dictionary,
+             (
+                const surfaceFilmModel& owner,
+                const dictionary& dict
+             ),
+             (owner, dict)
+         );
+
+    // Constructors
+
+        //- Construct null
+        force(const surfaceFilmModel& owner);
+
+        //- Construct from type name, dictionary and surface film model
+        force
+        (
+            const word& type,
+            const surfaceFilmModel& owner,
+            const dictionary& dict
+        );
+
+
+    // Selectors
+
+        //- Return a reference to the selected force model
+        static autoPtr<force> New
+        (
+            const surfaceFilmModel& owner,
+            const dictionary& dict,
+            const word& mdoelType
+        );
+
+
+    //- Destructor
+    virtual ~force();
+
+
+    // Member Functions
+
+        // Evolution
+
+            //- Correct
+            virtual tmp<fvVectorMatrix> correct(volVectorField& U) = 0;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace surfaceFilmModels
+} // End namespace regionModels
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/regionModels/surfaceFilmModels/submodels/kinematic/force/force/forceNew.C

@@ -0,0 +1,77 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2011 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 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 "force.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace regionModels
+{
+namespace surfaceFilmModels
+{
+
+// * * * * * * * * * * * * * * * * Selectors * * * * * * * * * * * * * * * * //
+
+autoPtr<force> force::New
+(
+    const surfaceFilmModel& model,
+    const dictionary& dict,
+    const word& modelType
+)
+{
+    Info<< "        " << modelType << endl;
+
+    dictionaryConstructorTable::iterator cstrIter =
+        dictionaryConstructorTablePtr_->find(modelType);
+
+    if (cstrIter == dictionaryConstructorTablePtr_->end())
+    {
+        FatalErrorIn
+        (
+            "force::New"
+            "("
+                "const surfaceFilmModel&, "
+                "const dictionary&, "
+                "const word&"
+            ")"
+        )   << "Unknown force type " << modelType
+            << nl << nl << "Valid force types are:" << nl
+            << dictionaryConstructorTablePtr_->toc()
+            << exit(FatalError);
+    }
+
+    return autoPtr<force>(cstrIter()(model, dict));
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace surfaceFilmModels
+} // End namespace regionModels
+} // End namespace Foam
+
+// ************************************************************************* //

src/regionModels/surfaceFilmModels/submodels/kinematic/force/forceList/forceList.C

@@ -0,0 +1,101 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2011 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 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 "forceList.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace regionModels
+{
+namespace surfaceFilmModels
+{
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+forceList::forceList(const surfaceFilmModel& owner)
+:
+    PtrList<force>()
+{}
+
+
+forceList::forceList
+(
+    const surfaceFilmModel& owner,
+    const dictionary& dict
+)
+:
+    PtrList<force>()
+{
+    const wordList models(dict.lookup("forces"));
+
+    Info<< "    Selecting film force models" << endl;
+    if (models.size() > 0)
+    {
+        this->setSize(models.size());
+
+        forAll(models, i)
+        {
+            set(i, force::New(owner, dict, models[i]));
+        }
+    }
+    else
+    {
+        Info<< "        none" << endl;
+    }
+}
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+forceList::~forceList()
+{}
+
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+
+tmp<fvVectorMatrix> forceList::correct(volVectorField& U)
+{
+    tmp<fvVectorMatrix> tResult
+    (
+        new fvVectorMatrix(U, dimForce/dimArea*dimVolume)
+    );
+    fvVectorMatrix& result = tResult();
+
+    forAll(*this, i)
+    {
+        result += this->operator[](i).correct(U);
+    }
+
+    return tResult;
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace surfaceFilmModels
+} // End namespace regionModels
+} // End namespace Foam
+
+// ************************************************************************* //

src/regionModels/surfaceFilmModels/submodels/kinematic/force/forceList/forceList.H

@@ -0,0 +1,95 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2011 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 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::forceList
+
+Description
+    List container for film sources
+
+SourceFiles
+    forceList.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef forceList_H
+#define forceList_H
+
+#include "PtrList.H"
+#include "force.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace regionModels
+{
+namespace surfaceFilmModels
+{
+
+/*---------------------------------------------------------------------------*\
+                       Class forceList Declaration
+\*---------------------------------------------------------------------------*/
+
+class forceList
+:
+    public PtrList<force>
+{
+public:
+
+    // Constructors
+
+        //- Construct null
+        forceList(const surfaceFilmModel& owner);
+
+        //- Construct from type name, dictionary and surface film model
+        forceList
+        (
+            const surfaceFilmModel& owner,
+            const dictionary& dict
+        );
+
+
+    //- Destructor
+    virtual ~forceList();
+
+
+    // Member functions
+
+        //- Return (net) force system
+        tmp<fvVectorMatrix> correct(volVectorField& U);
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace surfaceFilmModels
+} // End namespace regionModels
+} // End namespace Foam
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/regionModels/surfaceFilmModels/submodels/kinematic/force/surfaceShearForce/surfaceShearForce.C

@@ -0,0 +1,100 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2011 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 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 "surfaceShearForce.H"
+#include "addToRunTimeSelectionTable.H"
+#include "fvmSup.H"
+#include "kinematicSingleLayer.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace regionModels
+{
+namespace surfaceFilmModels
+{
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+defineTypeNameAndDebug(surfaceShearForce, 0);
+addToRunTimeSelectionTable(force, surfaceShearForce, dictionary);
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+surfaceShearForce::surfaceShearForce
+(
+    const surfaceFilmModel& owner,
+    const dictionary& dict
+)
+:
+    force(typeName, owner, dict),
+    Cf_(readScalar(coeffs_.lookup("Cf")))
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+surfaceShearForce::~surfaceShearForce()
+{}
+
+
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+
+tmp<fvVectorMatrix> surfaceShearForce::correct(volVectorField& U)
+{
+    const kinematicSingleLayer& film =
+        static_cast<const kinematicSingleLayer&>(owner_);
+
+    const volScalarField& rho = film.rho();
+    const volScalarField& mu = film.mu();
+    const volVectorField& Us = film.Us();
+    const volVectorField& Uw = film.Uw();
+    const volScalarField& delta = film.delta();
+
+    // Calculate shear stress
+    volScalarField Cs("Cs", rho*Cf_*mag(Us - U));
+    volScalarField Cw
+    (
+        "Cw",
+        mu/(0.3333*(delta + dimensionedScalar("SMALL", dimLength, SMALL)))
+    );
+    Cw.min(1.0e+06);
+
+    return
+    (
+       - fvm::Sp(Cs, U) + Cs*Us // surface contribution
+       - fvm::Sp(Cw, U) + Cw*Uw // wall contribution
+    );
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace surfaceFilmModels
+} // End namespace regionModels
+} // End namespace Foam
+
+// ************************************************************************* //

src/regionModels/surfaceFilmModels/submodels/kinematic/force/surfaceShearForce/surfaceShearForce.H

@@ -0,0 +1,114 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2011 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 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::surfaceShearForce
+
+Description
+    Film surface shear force
+
+SourceFiles
+    surfaceShearForce.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef surfaceShearForce_H
+#define surfaceShearForce_H
+
+#include "force.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace regionModels
+{
+namespace surfaceFilmModels
+{
+
+/*---------------------------------------------------------------------------*\
+                      Class surfaceShearForce Declaration
+\*---------------------------------------------------------------------------*/
+
+class surfaceShearForce
+:
+    public force
+{
+private:
+
+    // Private Data
+
+        //- Surface roughness coefficient
+        scalar Cf_;
+
+
+
+    // Private member functions
+
+        //- Disallow default bitwise copy construct
+        surfaceShearForce(const surfaceShearForce&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const surfaceShearForce&);
+
+
+public:
+
+    //- Runtime type information
+    TypeName("surfaceShear");
+
+
+    // Constructors
+
+        //- Construct from surface film model
+        surfaceShearForce
+        (
+            const surfaceFilmModel& owner,
+            const dictionary& dict
+        );
+
+
+    //- Destructor
+    virtual ~surfaceShearForce();
+
+
+    // Member Functions
+
+        // Evolution
+
+            //- Correct
+            virtual tmp<fvVectorMatrix> correct(volVectorField& U);
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace surfaceFilmModels
+} // End namespace regionModels
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/regionModels/surfaceFilmModels/submodels/kinematic/force/thermocapillaryForce/thermocapillaryForce.C

@@ -0,0 +1,83 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2011 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 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 "thermocapillaryForce.H"
+#include "addToRunTimeSelectionTable.H"
+#include "fvcGrad.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace regionModels
+{
+namespace surfaceFilmModels
+{
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+defineTypeNameAndDebug(thermocapillaryForce, 0);
+addToRunTimeSelectionTable(force, thermocapillaryForce, dictionary);
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+thermocapillaryForce::thermocapillaryForce
+(
+    const surfaceFilmModel& owner,
+    const dictionary& dict
+)
+:
+    force(owner)
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+thermocapillaryForce::~thermocapillaryForce()
+{}
+
+
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+
+tmp<fvVectorMatrix> thermocapillaryForce::correct(volVectorField& U)
+{
+    const volScalarField& sigma = owner_.sigma();
+
+    tmp<fvVectorMatrix>
+        tfvm(new fvVectorMatrix(U, dimForce/dimArea*dimVolume));
+
+    tfvm() += fvc::grad(sigma);
+
+    return tfvm;
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace surfaceFilmModels
+} // End namespace regionModels
+} // End namespace Foam
+
+// ************************************************************************* //

src/regionModels/surfaceFilmModels/submodels/kinematic/force/thermocapillaryForce/thermocapillaryForce.H

@@ -0,0 +1,107 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2011 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 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::thermocapillaryForce
+
+Description
+    Thermocapillary force
+
+SourceFiles
+    thermocapillaryForce.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef thermocapillaryForce_H
+#define thermocapillaryForce_H
+
+#include "force.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace regionModels
+{
+namespace surfaceFilmModels
+{
+
+/*---------------------------------------------------------------------------*\
+                  Class thermocapillaryForce Declaration
+\*---------------------------------------------------------------------------*/
+
+class thermocapillaryForce
+:
+    public force
+{
+private:
+
+    // Private member functions
+
+        //- Disallow default bitwise copy construct
+        thermocapillaryForce(const thermocapillaryForce&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const thermocapillaryForce&);
+
+
+public:
+
+    //- Runtime type information
+    TypeName("thermocapillary");
+
+
+    // Constructors
+
+        //- Construct from surface film model
+        thermocapillaryForce
+        (
+            const surfaceFilmModel& owner,
+            const dictionary& dict
+        );
+
+
+    //- Destructor
+    virtual ~thermocapillaryForce();
+
+
+    // Member Functions
+
+        // Evolution
+
+            //- Correct
+            virtual tmp<fvVectorMatrix> correct(volVectorField& U);
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace surfaceFilmModels
+} // End namespace regionModels
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/regionModels/surfaceFilmModels/thermoSingleLayer/thermoSingleLayer.C

@@ -97,11 +97,13 @@ void thermoSingleLayer::correctThermoFields()
     {
         case tmConstant:
         {
-            rho_ == dimensionedScalar(coeffs_.lookup("rho0"));
+            const dictionary&
-            mu_ == dimensionedScalar(coeffs_.lookup("mu0"));
+                constDict(coeffs_.subDict("constantThermoCoeffs"));
-            sigma_ == dimensionedScalar(coeffs_.lookup("sigma0"));
+            rho_ == dimensionedScalar(constDict.lookup("rho0"));
-            Cp_ == dimensionedScalar(coeffs_.lookup("Cp0"));
+            mu_ == dimensionedScalar(constDict.lookup("mu0"));
-            kappa_ == dimensionedScalar(coeffs_.lookup("kappa0"));
+            sigma_ == dimensionedScalar(constDict.lookup("sigma0"));
+            Cp_ == dimensionedScalar(constDict.lookup("Cp0"));
+            kappa_ == dimensionedScalar(constDict.lookup("kappa0"));
 
             break;
         }

tutorials/lagrangian/reactingParcelFilmFoam/cylinder/constant/surfaceFilmProperties

@@ -15,34 +15,27 @@ FoamFile
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-//surfaceFilmModel kinematicSingleLayer;
 surfaceFilmModel thermoSingleLayer;
 
 regionName      wallFilmRegion;
 
 active          true;
 
-kinematicSingleLayerCoeffs
-{
-    thermoModel constant;
-    rho0        rho0 [1 -3 0 0 0] 1000;
-    mu0         mu0 [1 -1 -1 0 0] 1e-3;
-    sigma0      sigma0 [1 0 -2 0 0] 0.07;
-
-    deltaStable deltaStable [0 1 0 0 0] 0;
-    Cf          0.005;
-
-    injectionModels ();
-}
-
-
 thermoSingleLayerCoeffs
 {
     thermoModel singleComponent;
     liquid      H2O;
 
-    deltaStable deltaStable [0 1 0 0 0] 0;
+    forces
-    Cf          0.005;
+    (
+        surfaceShear
+        thermocapillary
+    );
+
+    surfaceShearCoeffs
+    {
+        Cf          0.005;
+    }
 
     injectionModels ();
 

tutorials/lagrangian/reactingParcelFilmFoam/hotBoxes/constant/surfaceFilmProperties

@@ -25,9 +25,8 @@ active          true;
 thermoSingleLayerCoeffs
 {
     thermoModel singleComponent; // constant
-    liquid      H2O;
 
-    Cf          0.005;
+    liquid      H2O;
 
     radiationModel  none;
 
@@ -52,6 +51,35 @@ thermoSingleLayerCoeffs
         }
     }
 
+    forces
+    (
+        surfaceShear
+        thermocapillary
+        contactAngle
+    );
+
+    surfaceShearCoeffs
+    {
+        Cf          0.005;
+    }
+
+    contactAngleCoeffs
+    {
+        deltaWet        1e-4;
+        Ccf             0.085;
+        contactAngleDistribution
+        {
+            type            normal;
+            normalDistribution
+            {
+                minValue        50;
+                maxValue        100;
+                expectation     75;
+                variance        100;
+            }
+        }
+    }
+
     injectionModels
     (
         curvatureSeparation

tutorials/lagrangian/reactingParcelFilmFoam/hotBoxes/system/wallFilmRegion.org/fvSchemes

@@ -42,6 +42,7 @@ gradSchemes
     snGradCorr(pp) Gauss linear;
     snGradCorr(pu) Gauss linear;
     grad(nHat)  Gauss linear;
+    grad(alpha) Gauss linear;
 }
 
 laplacianSchemes

tutorials/lagrangian/reactingParcelFilmFoam/splashPanel/constant/surfaceFilmProperties

@@ -15,32 +15,28 @@ FoamFile
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-//surfaceFilmModel kinematicSingleLayer;
 surfaceFilmModel thermoSingleLayer;
 
 regionName      wallFilmRegion;
 
 active          true;
 
-kinematicSingleLayerCoeffs
-{
-    thermoModel constant;
-    rho0        rho0 [1 -3 0 0 0] 1000;
-    mu0         mu0 [1 -1 -1 0 0] 1e-3;
-    sigma0      sigma0 [1 0 -2 0 0] 0.07;
-
-    Cf          0.005;
-
-    injectionModels ();
-}
-
-
 thermoSingleLayerCoeffs
 {
     thermoModel singleComponent;
+
     liquid      H2O;
 
-    Cf          0.005;
+    forces
+    (
+        surfaceShear
+        thermocapillary
+    );
+
+    surfaceShearCoeffs
+    {
+        Cf          0.005;
+    }
 
     injectionModels
     ();