surfaceFilmModels::kinematicSingleLayer: Resolved issue with hydrostatic pressure contribution at corners

src/regionModels/surfaceFilmModels/derivedFvPatchFields/inclinedFilmNusseltHeight/inclinedFilmNusseltHeightFvPatchScalarField.C

@@ -132,8 +132,14 @@ void Foam::inclinedFilmNusseltHeightFvPatchScalarField::updateCoeffs()
     // note: normal pointing into the domain
     const vectorField n(-patch().nf());
 
-    // TODO: currently re-evaluating the entire gTan field to return this patch
+    const scalarField gTan
-    const scalarField gTan(film.gTan()().boundaryField()[patchi] & n);
+    (
+        (
+            film.g().value()
+          - film.nHat().boundaryField()[patchi]
+           *(film.g().value() & film.nHat().boundaryField()[patchi])
+        ) & n
+    );
 
     if (patch().size() && (max(mag(gTan)) < small))
     {

src/regionModels/surfaceFilmModels/derivedFvPatchFields/inclinedFilmNusseltInletVelocity/inclinedFilmNusseltInletVelocityFvPatchVectorField.C

@@ -113,7 +113,7 @@ void Foam::inclinedFilmNusseltInletVelocityFvPatchVectorField::updateCoeffs()
 
     const label patchi = patch().index();
 
-    // retrieve the film region from the database
+    // Retrieve the film region from the database
 
     const regionModels::regionModel& region =
         db().time().lookupObject<regionModels::regionModel>
@@ -127,12 +127,18 @@ void Foam::inclinedFilmNusseltInletVelocityFvPatchVectorField::updateCoeffs()
             const regionModels::surfaceFilmModels::kinematicSingleLayer&
         >(region);
 
-    // calculate the vector tangential to the patch
+    // Calculate the vector tangential to the patch
     // note: normal pointing into the domain
     const vectorField n(-patch().nf());
 
-    // TODO: currently re-evaluating the entire gTan field to return this patch
+    const scalarField gTan
-    const scalarField gTan(film.gTan()().boundaryField()[patchi] & n);
+    (
+        (
+            film.g().value()
+          - film.nHat().boundaryField()[patchi]
+           *(film.g().value() & film.nHat().boundaryField()[patchi])
+        ) & n
+    );
 
     if (patch().size() && (max(mag(gTan)) < small))
     {

src/regionModels/surfaceFilmModels/kinematicSingleLayer/kinematicSingleLayer.C

@@ -145,7 +145,7 @@ void kinematicSingleLayer::transferPrimaryRegionSourceFields()
 }
 
 
-tmp<volScalarField> kinematicSingleLayer::pu()
+tmp<volScalarField> kinematicSingleLayer::Su()
 {
     tmp<volScalarField> tpSp
     (
@@ -163,7 +163,7 @@ tmp<volScalarField> kinematicSingleLayer::pu()
 
     return volScalarField::New
     (
-        IOobject::modelName("pu", typeName),
+        IOobject::modelName("Su", typeName),
         pPrimary_                      // Pressure (mapped from primary region)
       - tpSp                           // Accumulated particle impingement
       - fvc::laplacian(sigma_, delta_) // Surface tension
@@ -171,14 +171,9 @@ tmp<volScalarField> kinematicSingleLayer::pu()
 }
 
 
-tmp<volScalarField> kinematicSingleLayer::pp()
+tmp<volScalarField> kinematicSingleLayer::ph() const
 {
-     // Hydrostatic effect
+    return -rho_*min(nHat() & g_, dimensionedScalar(g_.dimensions(), 0))*VbyA();
-    return volScalarField::New
-    (
-        IOobject::modelName("pp", typeName),
-        -rho_*gNormClipped()*VbyA()
-    );
 }
 
 
@@ -209,12 +204,7 @@ void kinematicSingleLayer::predictDelta()
 {
     DebugInFunction << endl;
 
-    solve
+    solve(fvm::ddt(rho_, alpha_) + fvc::div(phi_) == -rhoSp_);
-    (
-        fvm::ddt(rho_, alpha_) + fvc::div(phi_)
-     ==
-       -rhoSp_
-    );
 
     // Bound film volume fraction
     alpha_.max(0);
@@ -284,8 +274,7 @@ void kinematicSingleLayer::updateSurfaceVelocities()
 
 tmp<Foam::fvVectorMatrix> kinematicSingleLayer::solveMomentum
 (
-    const volScalarField& pu,
+    const volScalarField& pu
-    const volScalarField& pp
 )
 {
     DebugInFunction << endl;
@@ -327,17 +316,11 @@ tmp<Foam::fvVectorMatrix> kinematicSingleLayer::solveMomentum
                 (
                     fvc::interpolate(alpha_)
                    *(
-                        regionMesh().magSf()
+                        (
-                       *(
                             fvc::snGrad(pu, "snGrad(p)")
-
+                          + fvc::interpolate(ph())*fvc::snGrad(alpha_)
-                          + fvc::interpolate(alpha_)
+                        )*regionMesh().magSf()
-                           *fvc::snGrad(pp, "snGrad(p)")
+                      - fvc::interpolate(rho_)*(g_ & regionMesh().Sf())
-
-                          + fvc::interpolate(pp)
-                           *fvc::snGrad(alpha_)
-                        )
-                      - fvc::flux(rho_*gTan())
                     ), 0
                 )
             )
@@ -355,9 +338,8 @@ tmp<Foam::fvVectorMatrix> kinematicSingleLayer::solveMomentum
 
 void kinematicSingleLayer::solveAlpha
 (
-    const volScalarField& pu,
+    const fvVectorMatrix& UEqn,
-    const volScalarField& pp,
+    const volScalarField& pu
-    const fvVectorMatrix& UEqn
 )
 {
     DebugInFunction << endl;
@@ -368,7 +350,7 @@ void kinematicSingleLayer::solveAlpha
     const surfaceScalarField alphaf(fvc::interpolate(alpha_));
     const surfaceScalarField rhof(fvc::interpolate(rho_));
     const surfaceScalarField alpharAUf(fvc::interpolate(alpha_*rAU));
-    const surfaceScalarField ppf(fvc::interpolate(pp));
+    const surfaceScalarField phf(fvc::interpolate(ph()));
 
     const surfaceScalarField phiu
     (
@@ -376,10 +358,7 @@ void kinematicSingleLayer::solveAlpha
         (
             constrainFilmField
             (
-                (
+                fvc::snGrad(pu, "snGrad(p)")*regionMesh().magSf()
-                    fvc::snGrad(pu, "snGrad(p)")
-                  + alphaf*fvc::snGrad(pp, "snGrad(p)")
-                )*regionMesh().magSf()
               - rhof*(g_ & regionMesh().Sf()),
                 0
             )
@@ -392,10 +371,10 @@ void kinematicSingleLayer::solveAlpha
         constrainFilmField(rhof*(fvc::flux(HbyA) - alpharAUf*phiu), 0)
     );
 
-    const surfaceScalarField ddrhorAUppf
+    const surfaceScalarField ddrhorAUphf
     (
         "alphaCoeff",
-        alphaf*rhof*alpharAUf*ppf
+        alphaf*rhof*alpharAUf*phf
     );
 
     regionMesh().setFluxRequired(alpha_.name());
@@ -407,7 +386,7 @@ void kinematicSingleLayer::solveAlpha
         (
             fvm::ddt(rho_, alpha_)
           + fvm::div(phid, alpha_)
-          - fvm::laplacian(ddrhorAUppf, alpha_)
+          - fvm::laplacian(ddrhorAUphf, alpha_)
          ==
            -rhoSp_
         );
@@ -422,7 +401,7 @@ void kinematicSingleLayer::solveAlpha
             (
                 constrainFilmField
                 (
-                    ppf*fvc::snGrad(alpha_)*regionMesh().magSf(),
+                    phf*fvc::snGrad(alpha_)*regionMesh().magSf(),
                     0
                 )
             );
@@ -940,21 +919,18 @@ void kinematicSingleLayer::evolveRegion()
     // Predict delta_ from continuity with updated source
     predictDelta();
 
-    // Implicit pressure source coefficient - constant
-    const volScalarField pp(this->pp());
-
     while (pimple_.loop())
     {
         // Explicit pressure source contribution - varies with delta
-        const volScalarField pu(this->pu());
+        const volScalarField Su(this->Su());
 
         // Solve for momentum
-        const fvVectorMatrix UEqn(solveMomentum(pu, pp));
+        const fvVectorMatrix UEqn(solveMomentum(Su));
 
         // Film thickness correction loop
         while (pimple_.correct())
         {
-            solveAlpha(pu, pp, UEqn);
+            solveAlpha(UEqn, Su);
         }
     }
 

src/regionModels/surfaceFilmModels/kinematicSingleLayer/kinematicSingleLayer.H

@@ -231,11 +231,11 @@ protected:
         //- Transfer source fields from the primary region to the film region
         virtual void transferPrimaryRegionSourceFields();
 
-        //- Explicit pressure source contribution
+        //- Hydrostatic pressure
-        virtual tmp<volScalarField> pu();
+        tmp<volScalarField> ph() const;
 
-        //- Implicit pressure source coefficient
+        //- Explicit momentum source
-        virtual tmp<volScalarField> pp();
+        virtual tmp<volScalarField> Su();
 
         //- Correct film coverage field
         virtual void correctCoverage();
@@ -270,16 +270,14 @@ protected:
             //- Solve for film velocity
             virtual tmp<fvVectorMatrix> solveMomentum
             (
-                const volScalarField& pu,
+                const volScalarField& pu
-                const volScalarField& pp
             );
 
             //- Solve for film volume fraction and thickness
             virtual void solveAlpha
             (
-                const volScalarField& pu,
+                const fvVectorMatrix& UEqn,
-                const volScalarField& pp,
+                const volScalarField& Su
-                const fvVectorMatrix& UEqn
             );
 
 
@@ -461,16 +459,6 @@ public:
             //- Return the change in film mass due to sources/sinks
             inline tmp<volScalarField::Internal> deltaMass() const;
 
-            //- Return the gravity normal-to-patch component contribution
-            inline tmp<volScalarField> gNorm() const;
-
-            //- Return the gravity normal-to-patch component contribution
-            //  Clipped so that only non-zero if g & nHat_ < 0
-            inline tmp<volScalarField> gNormClipped() const;
-
-            //- Return the gravity tangential component contributions
-            inline tmp<volVectorField> gTan() const;
-
 
         // Evolution
 

src/regionModels/surfaceFilmModels/kinematicSingleLayer/kinematicSingleLayerI.H

@@ -231,32 +231,6 @@ inline tmp<volScalarField::Internal> kinematicSingleLayer::deltaMass() const
 }
 
 
-inline tmp<volScalarField> kinematicSingleLayer::gNorm() const
-{
-    return volScalarField::New("gNorm", g_ & nHat());
-}
-
-
-inline tmp<volScalarField> kinematicSingleLayer::gNormClipped() const
-{
-    tmp<volScalarField> tgNormClipped
-    (
-        volScalarField::New("gNormClipped", g_ & nHat())
-    );
-
-    volScalarField& gNormClipped = tgNormClipped.ref();
-    gNormClipped.min(0.0);
-
-    return tgNormClipped;
-}
-
-
-inline tmp<volVectorField> kinematicSingleLayer::gTan() const
-{
-    return volVectorField::New("gTan", g_ - nHat()*gNorm());
-}
-
-
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 } // End namespace surfaceFilmModels

src/regionModels/surfaceFilmModels/submodels/kinematic/injectionModel/BrunDrippingInjection/BrunDrippingInjection.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2016-2018 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2016-2019 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -76,7 +76,7 @@ void BrunDrippingInjection::correct
         refCast<const kinematicSingleLayer>(this->film());
 
     // Calculate available dripping mass
-    tmp<volScalarField> tsinAlpha(film.gNorm()/mag(film.g()));
+    tmp<volScalarField> tsinAlpha((film.g()/mag(film.g())) & film.nHat());
     const scalarField& sinAlpha = tsinAlpha();
     const scalarField& magSf = film.magSf();
 

src/regionModels/surfaceFilmModels/submodels/kinematic/injectionModel/drippingInjection/drippingInjection.C

@@ -91,8 +91,7 @@ void drippingInjection::correct
     const scalar pi = constant::mathematical::pi;
 
     // calculate available dripping mass
-    tmp<volScalarField> tgNorm(film.gNorm());
+    const scalarField gNorm(film.g() & film.nHat()());
-    const scalarField& gNorm = tgNorm();
     const scalarField& magSf = film.magSf();
 
     const scalarField& delta = film.delta();

src/regionModels/surfaceFilmModels/thermoSingleLayer/thermoSingleLayer.C

@@ -631,16 +631,13 @@ void thermoSingleLayer::evolveRegion()
     // Predict delta_ from continuity with updated source
     predictDelta();
 
-    // Implicit pressure source coefficient - constant
-    const volScalarField pp(this->pp());
-
     while (pimple_.loop())
     {
-        // Explicit pressure source contribution - varies with delta
+        // Explicit momentum source
-        const volScalarField pu(this->pu());
+        const volScalarField Su(this->Su());
 
         // Solve for momentum for U_
-        const fvVectorMatrix UEqn(solveMomentum(pu, pp));
+        const fvVectorMatrix UEqn(solveMomentum(Su));
 
         // Solve energy for h_ - also updates thermo
         solveEnergy();
@@ -648,7 +645,7 @@ void thermoSingleLayer::evolveRegion()
         // Film thickness correction loop
         while (pimple_.correct())
         {
-            solveAlpha(pu, pp, UEqn);
+            solveAlpha(UEqn, Su);
         }
     }