ENH: Adding shear stress to the film from the wall function

TUT: inclinedPlaneFilm/pitzDailyWithSprinklers: add shearStress model

src/regionFaModels/liquidFilm/liquidFilmBase.C

@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2020 OpenCFD Ltd.
+    Copyright (C) 2020-2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -551,6 +551,13 @@ scalar liquidFilmBase::pRef()
     return pRef_;
 }
 
+
+word liquidFilmBase::UName() const
+{
+    return UName_;
+}
+
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 } // End namespace areaSurfaceFilmModels

src/regionFaModels/liquidFilm/liquidFilmBase.H

@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2020 OpenCFD Ltd.
+    Copyright (C) 2020-2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -261,6 +261,9 @@ public:
             //- Access to pRef
             scalar pRef();
 
+            //- Name of the U field
+            word UName() const;
+
 
             // Transfer fields - to the primary region (lagragian injection)
 

src/regionFaModels/liquidFilm/subModels/kinematic/filmTurbulenceModel/filmTurbulenceModel/filmTurbulenceModel.C

@@ -27,6 +27,8 @@ License
 
 #include "filmTurbulenceModel.H"
 #include "gravityMeshObject.H"
+#include "turbulentTransportModel.H"
+#include "turbulentFluidThermoModel.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -55,6 +57,17 @@ filmTurbulenceModel::frictionMethodTypeNames_
 };
 
 
+const Enum
+<
+    filmTurbulenceModel::shearMethodType
+>
+filmTurbulenceModel::shearMethodTypeNames_
+{
+    { shearMethodType::msimple, "simple" },
+    { shearMethodType::mwallFunction, "wallFunction" }
+};
+
+
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 filmTurbulenceModel::filmTurbulenceModel
@@ -66,8 +79,16 @@ filmTurbulenceModel::filmTurbulenceModel
 :
     film_(film),
     dict_(dict.subDict(modelType + "Coeffs")),
-    method_(frictionMethodTypeNames_.get("friction", dict_))
+    method_(frictionMethodTypeNames_.get("friction", dict_)),
-{}
+    shearMethod_(shearMethodTypeNames_.get("shearStress", dict_)),
+    rhoName_(dict_.getOrDefault<word>("rho", "rho")),
+    rhoRef_(VGREAT)
+{
+    if (rhoName_ == "rhoInf")
+    {
+        rhoRef_ = dict_.get<scalar>("rhoInf");
+    }
+}
 
 
 // * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
@@ -168,6 +189,167 @@ tmp<areaScalarField> filmTurbulenceModel::Cw() const
 }
 
 
+tmp<faVectorMatrix> filmTurbulenceModel::primaryRegionFriction
+(
+    areaVectorField& U
+) const
+{
+    tmp<faVectorMatrix> tshearStress
+    (
+        new faVectorMatrix(U, sqr(U.dimensions())*sqr(dimLength))
+    );
+
+    switch (shearMethod_)
+    {
+        case msimple:
+        {
+            tmp<areaVectorField> Up = film_.Up();
+
+            const dimensionedScalar Cf
+            (
+                "Cf",
+                dimVelocity,
+                dict_.get<scalar>("Cf")
+            );
+
+            tshearStress.ref() += - fam::Sp(Cf, U) + Cf*Up();
+
+            break;
+        }
+        case mwallFunction:
+        {
+            tmp<volSymmTensorField> tdevRhoReff = devRhoReff();
+
+            const volSymmTensorField::Boundary& devRhoReffb
+                = tdevRhoReff().boundaryField();
+
+            const label patchi = film_.patchID();
+
+            const surfaceVectorField::Boundary& Sfb =
+                film_.primaryMesh().Sf().boundaryField();
+
+            vectorField fT(Sfb[patchi] & devRhoReffb[patchi]);
+
+            const vectorField& nHat =
+                film_.regionMesh().faceAreaNormals().internalField();
+
+            // Substract normal component
+            fT -= nHat*(fT & nHat);
+
+            auto taForce = tmp<areaVectorField>::New
+            (
+                IOobject
+                (
+                    "taForce",
+                    film_.primaryMesh().time().timeName(),
+                    film_.primaryMesh()
+                ),
+                film_.regionMesh(),
+                dimensionedVector(sqr(dimVelocity), Zero)
+            );
+            vectorField& aForce = taForce.ref().primitiveFieldRef();
+
+            // Map ft to surface
+            const vectorField afT(film_.vsm().mapToSurface(fT));
+
+            const DimensionedField<scalar, areaMesh>& magSf =
+                film_.regionMesh().S();
+
+            aForce = afT/(film_.rho().primitiveField()*magSf);
+
+            tshearStress.ref() += taForce();
+
+            if (film_.regionMesh().time().writeTime())
+            {
+                taForce().write();
+            }
+
+            break;
+        }
+    }
+
+    return tshearStress;
+}
+
+
+tmp<Foam::volSymmTensorField> filmTurbulenceModel::devRhoReff() const
+{
+    typedef compressible::turbulenceModel cmpTurbModel;
+    typedef incompressible::turbulenceModel icoTurbModel;
+
+    const fvMesh& m = film_.primaryMesh();
+
+    const auto& U = m.lookupObject<volVectorField>(film_.UName());
+
+    if (m.foundObject<cmpTurbModel>(cmpTurbModel::propertiesName))
+    {
+        const auto& turb =
+            m.lookupObject<cmpTurbModel>(cmpTurbModel::propertiesName);
+
+        return turb.devRhoReff();
+    }
+    else if (m.foundObject<icoTurbModel>(icoTurbModel::propertiesName))
+    {
+        const auto& turb =
+            m.lookupObject<icoTurbModel>(icoTurbModel::propertiesName);
+
+        return rho()*turb.devReff();
+    }
+    else if (m.foundObject<fluidThermo>(fluidThermo::dictName))
+    {
+        const auto& thermo =
+            m.lookupObject<fluidThermo>(fluidThermo::dictName);
+
+        return -thermo.mu()*dev(twoSymm(fvc::grad(U)));
+    }
+    else if (m.foundObject<transportModel>("transportProperties"))
+    {
+        const auto& laminarT =
+            m.lookupObject<transportModel>("transportProperties");
+
+        return -rho()*laminarT.nu()*dev(twoSymm(fvc::grad(U)));
+    }
+    else if (m.foundObject<dictionary>("transportProperties"))
+    {
+        const auto& transportProperties =
+            m.lookupObject<dictionary>("transportProperties");
+
+        const dimensionedScalar nu("nu", dimViscosity, transportProperties);
+
+        return -rho()*nu*dev(twoSymm(fvc::grad(U)));
+    }
+    else
+    {
+        FatalErrorInFunction
+            << "No valid model for viscous stress calculation"
+            << exit(FatalError);
+
+        return volSymmTensorField::null();
+    }
+}
+
+
+tmp<Foam::volScalarField> filmTurbulenceModel::rho() const
+{
+    const fvMesh& m = film_.primaryMesh();
+    if (rhoName_ == "rhoInf")
+    {
+        return tmp<volScalarField>::New
+        (
+            IOobject
+            (
+                "rho",
+                m.time().timeName(),
+                m
+            ),
+            m,
+            dimensionedScalar(dimDensity, rhoRef_)
+        );
+    }
+
+    return m.lookupObject<volScalarField>(rhoName_);
+}
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 } // End namespace areaSurfaceFilmModels

src/regionFaModels/liquidFilm/subModels/kinematic/filmTurbulenceModel/filmTurbulenceModel/filmTurbulenceModel.H

@@ -82,6 +82,13 @@ public:
             mManningStrickler
         };
 
+        //- Options for the shear stress models
+        enum shearMethodType
+        {
+            msimple,
+            mwallFunction
+        };
+
 
 protected:
 
@@ -93,12 +100,24 @@ protected:
         //- Names for friction models
         static const Enum<frictionMethodType> frictionMethodTypeNames_;
 
+        //- Names for shear stress models
+        static const Enum<shearMethodType> shearMethodTypeNames_;
+
         //- Model dictionary
         const dictionary dict_;
 
         //- Method used
         const frictionMethodType method_;
 
+        //- Shear method used
+        const shearMethodType shearMethod_;
+
+        //- Name of density field (optional)
+        word rhoName_;
+
+        //- Reference density needed for incompressible calculations
+        scalar rhoRef_;
+
 
 public:
 
@@ -154,7 +173,19 @@ public:
             const liquidFilmBase& film() const;
 
 
-        // Evaluation
+        // Turbulence
+
+            //- Return the effective viscous stress (laminar + turbulent)
+            tmp<volSymmTensorField> devRhoReff() const;
+
+            //- Return primary region friction
+            tmp<faVectorMatrix> primaryRegionFriction
+            (
+                areaVectorField& U
+            ) const;
+
+            //- Return rho if specified otherwise rhoRef
+            tmp<volScalarField> rho() const;
 
             //- Return the wall film surface friction
             virtual tmp<areaScalarField> Cw() const;
@@ -162,6 +193,9 @@ public:
             //- Return the film turbulence viscosity
             virtual tmp<areaScalarField> mut() const = 0;
 
+
+        // Evaluation
+
             //- Correct/update the model
             virtual void correct() = 0;
 

src/regionFaModels/liquidFilm/subModels/kinematic/filmTurbulenceModel/laminar/laminar.C

@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2020 OpenCFD Ltd.
+    Copyright (C) 2020-2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -51,14 +51,7 @@ laminar::laminar
     const dictionary& dict
 )
 :
-    filmTurbulenceModel(type(), film, dict),
+    filmTurbulenceModel(type(), film, dict)
-    Cf_(dict_.get<scalar>("Cf"))
-{}
-
-
-// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
-
-laminar::~laminar()
 {}
 
 
@@ -66,20 +59,19 @@ laminar::~laminar()
 
 tmp<areaScalarField> laminar::mut() const
 {
-    auto tmut =
+    auto tmut = tmp<areaScalarField>::New
-        tmp<areaScalarField>::New
+    (
+        IOobject
         (
-            IOobject
+            "mut",
-            (
+            film().primaryMesh().time().timeName(),
-                "mut",
+            film().primaryMesh(),
-                film().primaryMesh().time().timeName(),
+            IOobject::NO_READ,
-                film().primaryMesh(),
+            IOobject::NO_WRITE
-                IOobject::NO_READ,
+        ),
-                IOobject::NO_WRITE
+        film().regionMesh(),
-            ),
+        dimensionedScalar(dimMass/dimLength/dimTime)
-            film().regionMesh(),
+    );
-            dimensionedScalar(dimMass/dimLength/dimTime)
-        );
 
     return tmut;
 }
@@ -90,19 +82,19 @@ void laminar::correct()
 
 
 tmp<faVectorMatrix> laminar::Su(areaVectorField& U) const
+{
+    return primaryRegionFriction(U) + wallFriction(U);
+}
+
+
+tmp<faVectorMatrix> laminar::wallFriction(areaVectorField& U) const
 {
     // local references to film fields
     tmp<areaVectorField> Uw = film_.Uw();
-    tmp<areaVectorField> Up = film_.Up();
-
-    // employ simple coeff-based model
-    const dimensionedScalar Cf("Cf", dimVelocity, Cf_);
-
     tmp<areaScalarField> wf = Cw();
 
     return
     (
-       - fam::Sp(Cf, U) + Cf*Up()     // surface contribution
        - fam::Sp(wf(), U) + wf()*Uw() // wall contribution
     );
 }

src/regionFaModels/liquidFilm/subModels/kinematic/filmTurbulenceModel/laminar/laminar.H

@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2020 OpenCFD Ltd.
+    Copyright (C) 2020-2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -57,12 +57,6 @@ class laminar
 :
     public filmTurbulenceModel
 {
-    // Private Data
-
-        //- Surface roughness coefficient
-        scalar Cf_;
-
-
     // Private Member Functions
 
         //- No copy construct
@@ -85,21 +79,27 @@ public:
 
 
     //- Destructor
-    virtual ~laminar();
+    virtual ~laminar() = default;
 
 
     // Member Functions
 
-        // Evolution
+        // Turbulence
+
+            //- Wall friction
+            tmp<faVectorMatrix> wallFriction(areaVectorField& U) const;
 
             //- Return the film turbulence viscosity
             virtual tmp<areaScalarField> mut() const;
 
-            //- Correct/update the model
-            virtual void correct();
-
             //- Return the source for the film momentum equation
             virtual tmp<faVectorMatrix> Su(areaVectorField& U) const;
+
+
+        // Evaluation
+
+            //- Correct/update the model
+            virtual void correct();
 };
 
 

tutorials/incompressible/pimpleFoam/laminar/inclinedPlaneFilm/0/U

@@ -51,6 +51,7 @@ boundaryField
 
         laminarCoeffs
         {
+            shearStress     simple;
             friction        ManningStrickler;
             n               0.1; // Manning number
             Cf              0.9; // Gas friction

tutorials/lagrangian/kinematicParcelFoam/pitzDailyWithSprinklers/0/U

@@ -57,6 +57,7 @@ boundaryField
 
         laminarCoeffs
         {
+            shearStress     simple;
             friction        ManningStrickler; // Wall friction model
             n               0.005;  // Manning number
             Cf              0;      // Gas friction