ENH: Updated film surface shear force to include near-wall turbulence effects

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

@@ -27,6 +27,7 @@ License
 #include "addToRunTimeSelectionTable.H"
 #include "fvmSup.H"
 #include "kinematicSingleLayer.H"
+#include "turbulenceModel.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -65,27 +66,76 @@ surfaceShearForce::~surfaceShearForce()
 
 tmp<fvVectorMatrix> surfaceShearForce::correct(volVectorField& U)
 {
+    // local reference to film model
     const kinematicSingleLayer& film =
         static_cast<const kinematicSingleLayer&>(owner_);
 
-    const volScalarField& rho = film.rho();
+    // local references to film fields
     const volScalarField& mu = film.mu();
-    const volVectorField& Us = film.Us();
     const volVectorField& Uw = film.Uw();
     const volScalarField& delta = film.delta();
+    const volVectorField& Up = film.UPrimary();
 
-    // Calculate shear stress
+    // film surface linear coeff to apply to velocity
-    volScalarField Cs("Cs", rho*Cf_*mag(Us - U));
+    tmp<volScalarField> tCs;
-    volScalarField Cw
+
-    (
+    typedef compressible::turbulenceModel turbModel;
-        "Cw",
+    if (film.primaryMesh().foundObject<turbModel>("turbulenceProperties"))
-        mu/(0.3333*(delta + dimensionedScalar("SMALL", dimLength, SMALL)))
+    {
-    );
+        // local reference to turbulence model
+        const turbModel& turb =
+            film.primaryMesh().lookupObject<turbModel>("turbulenceProperties");
+
+        // calculate and store the stress on the primary region
+        const volSymmTensorField primaryReff(turb.devRhoReff());
+
+        // create stress field on film
+        // - note boundary condition types (mapped)
+        // - to map, the field name must be the same as the field on the
+        //   primary region
+        volSymmTensorField Reff
+        (
+            IOobject
+            (
+                primaryReff.name(),
+                film.regionMesh().time().timeName(),
+                film.regionMesh(),
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            film.regionMesh(),
+            dimensionedSymmTensor
+            (
+                "zero",
+                primaryReff.dimensions(),
+                symmTensor::zero
+            ),
+            film.mappedPushedFieldPatchTypes<symmTensor>()
+        );
+
+        // map stress from primary region to film region
+        Reff.correctBoundaryConditions();
+
+        dimensionedScalar U0("SMALL", U.dimensions(), SMALL);
+        tCs = Cf_*mag(-film.nHat() & Reff)/(mag(Up - U) + U0);
+    }
+    else
+    {
+        // laminar case - employ simple coeff-based model
+        const volScalarField& rho = film.rho();
+        tCs = Cf_*rho*mag(Up - U);
+    }
+    
+    dimensionedScalar d0("SMALL", delta.dimensions(), SMALL);
+
+    // linear coeffs to apply to velocity
+    const volScalarField& Cs = tCs();
+    volScalarField Cw("Cw", mu/(0.3333*(delta + d0)));
     Cw.min(1.0e+06);
 
     return
     (
-       - fvm::Sp(Cs, U) + Cs*Us // surface contribution
+       - fvm::Sp(Cs, U) + Cs*Up // surface contribution
        - fvm::Sp(Cw, U) + Cw*Uw // wall contribution
     );
 }