isothermalFilm: Corrected impingement pressure transfer

The previous implementation was dimensionally inconsistent and was
missing a factor of the VbyA field. This change will, in most cases,
reduce the total impingement pressure contribution.

applications/modules/isothermalFilm/fvModels/filmCloudTransfer/cloudFilmTransfer/CloudFilmTransfer/CloudFilmTransfer.C

@@ -83,21 +83,12 @@ void Foam::CloudFilmTransfer<CloudType>::absorbInteraction
     const parcelThermo& thermo =
         static_cast<const ThermoCloud<CloudType>&>(this->owner()).thermo();
 
-    // 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;
-
     const liquidProperties& liq = thermo.liquids().properties()[0];
 
     // Local pressure
@@ -107,8 +98,7 @@ void Foam::CloudFilmTransfer<CloudType>::absorbInteraction
     (
         facei,
         mass,                           // mass
-        mass*Ut,                        // tangential momentum
+        mass*Urel,                      // momentum
-        mass*mag(Un),                   // impingement pressure
         mass*liq.Hs(pc, p.T())          // energy
     );
 

applications/modules/isothermalFilm/fvModels/filmCloudTransfer/filmCloudTransfer.C

@@ -75,7 +75,6 @@ Foam::wordList Foam::fv::filmCloudTransfer::addSupFields() const
 {
     return wordList
     {
-        "pi",
         film_.alpha.name(),
         film_.thermo.he().name(),
         film_.U.name()
@@ -136,36 +135,6 @@ inline Foam::fv::filmCloudTransfer::CloudToFilmTransferRate
 }
 
 
-void Foam::fv::filmCloudTransfer::addSup
-(
-    fvMatrix<scalar>& eqn,
-    const word& fieldName
-) const
-{
-    if (debug)
-    {
-        Info<< type() << ": applying source to " << eqn.psi().name() << endl;
-    }
-
-    // Droplet impingement pressure
-    if (fieldName == "pi")
-    {
-        eqn +=
-            CloudToFilmTransferRate<scalar>
-            (
-                pressureFromCloud_,
-                dimPressure*dimVolume
-            );
-    }
-    else
-    {
-        FatalErrorInFunction
-            << "Support for field " << fieldName << " is not implemented"
-            << exit(FatalError);
-    }
-}
-
-
 void Foam::fv::filmCloudTransfer::addSup
 (
     const volScalarField& rho,
@@ -240,11 +209,20 @@ void Foam::fv::filmCloudTransfer::addSup
         Info<< type() << ": applying source to " << eqn.psi().name() << endl;
     }
 
-    eqn += CloudToFilmTransferRate<vector>(momentumFromCloud_, dimMomentum);
+    if (fieldName == film_.U.name())
-
-    if (ejection_.valid())
     {
-        eqn -= fvm::Sp(alpha()*rho()*ejection_->rate(), eqn.psi());
+        eqn += CloudToFilmTransferRate<vector>(momentumFromCloud_, dimMomentum);
+
+        if (ejection_.valid())
+        {
+            eqn -= fvm::Sp(alpha()*rho()*ejection_->rate(), eqn.psi());
+        }
+    }
+    else
+    {
+        FatalErrorInFunction
+            << "Support for field " << fieldName << " is not implemented"
+            << exit(FatalError);
     }
 }
 
@@ -260,13 +238,11 @@ void Foam::fv::filmCloudTransfer::resetFromCloudFields()
     {
         massFromCloud_.setSize(nCloudPatchFaces);
         momentumFromCloud_.setSize(nCloudPatchFaces);
-        pressureFromCloud_.setSize(nCloudPatchFaces);
         energyFromCloud_.setSize(nCloudPatchFaces);
     }
 
     massFromCloud_ = 0;
     momentumFromCloud_ = Zero;
-    pressureFromCloud_ = 0;
     energyFromCloud_ = 0;
 
     cloudFieldsTransferred_ = true;
@@ -281,13 +257,11 @@ void Foam::fv::filmCloudTransfer::parcelFromCloud
     const label facei,
     const scalar mass,
     const vector& momentum,
-    const scalar pressure,
     const scalar energy
 )
 {
     massFromCloud_[facei] += mass;
     momentumFromCloud_[facei] += momentum;
-    pressureFromCloud_[facei] += pressure;
     energyFromCloud_[facei] += energy;
 }
 

applications/modules/isothermalFilm/fvModels/filmCloudTransfer/filmCloudTransfer.H

@@ -75,7 +75,6 @@ class filmCloudTransfer
 
             scalarField massFromCloud_;
             vectorField momentumFromCloud_;
-            scalarField pressureFromCloud_;
             scalarField energyFromCloud_;
 
         //- Switch to ensure the ejection rate is not updated until
@@ -147,13 +146,6 @@ public:
 
         // Add explicit and implicit contributions to compressible equation
 
-            //- Add explicit droplet impingement contribution to pressure field
-            virtual void addSup
-            (
-                fvMatrix<scalar>& eqn,
-                const word& fieldName
-            ) const;
-
             //- Add explicit contribution to phase continuity
             virtual void addSup
             (
@@ -193,7 +185,6 @@ public:
                 const label facei,
                 const scalar mass,
                 const vector& momentum,
-                const scalar pressure,
                 const scalar energy
             );
 

applications/modules/isothermalFilm/momentumPredictor.C

@@ -76,8 +76,9 @@ Foam::solvers::isothermalFilm::pe() const
     // Update the pressure, mapping from the fluid region as required
     p.correctBoundaryConditions();
 
-    // Add the droplet impingement pressure
+    // Add the pressure caused normal momentum sources (e.g., parcels impinging
-    p.ref() += mesh.time().deltaT()*fvModels().source(p, "pi")().Su();
+    // with a normal velocity)
+    p.ref() += VbyA*(nHat & (fvModels().source(alpha, rho, U) & U));
 
     return p;
 }
@@ -90,6 +91,10 @@ void Foam::solvers::isothermalFilm::momentumPredictor()
     // Calculate the surface tension coefficient
     const volScalarField sigma(this->sigma());
 
+    // Get the momentum source and remove any normal components
+    fvVectorMatrix alphaRhoUsource(fvModels().source(alpha, rho, U));
+    alphaRhoUsource.source() -= nHat*(nHat & alphaRhoUsource.source());
+
     tUEqn =
     (
         fvm::ddt(alpha, rho, U) + fvm::div(alphaRhoPhi, U)
@@ -97,7 +102,7 @@ void Foam::solvers::isothermalFilm::momentumPredictor()
       + momentumTransport->divDevTau(U)
      ==
         contactForce(sigma)
-      + fvModels().source(alpha, rho, U)
+      + alphaRhoUsource
     );
     fvVectorMatrix& UEqn = tUEqn.ref();