fvModels::waveForcing: Added a continuity correction to momentum equation

to compensate for phase-change generated by the phase-fraction forcing.  This
removes spurious velocity generated around the wave.

src/waves/fvModels/forcing/forcing.C

@@ -116,57 +116,63 @@ void Foam::fv::forcing::readCoeffs()
 }
 
 
-Foam::tmp<Foam::volScalarField::Internal> Foam::fv::forcing::forceCoeff() const
+Foam::tmp<Foam::volScalarField::Internal> Foam::fv::forcing::scale() const
 {
-    tmp<volScalarField::Internal> tforceCoeff
+    tmp<volScalarField::Internal> tscale
     (
-        new volScalarField::Internal
+        volScalarField::Internal::New
         (
-            IOobject
-            (
-                typedName("forceCoeff"),
-                mesh().time().timeName(),
-                mesh()
-            ),
+            typedName("scale"),
             mesh(),
-            dimensionedScalar(lambda_.dimensions(), scale_.valid() ? 0 : 1)
+            dimensionedScalar(dimless, scale_.valid() ? 0 : 1)
         )
     );
-    scalarField& forceCoeff = tforceCoeff.ref();
+
+    scalarField& scale = tscale.ref();
 
     forAll(origins_, i)
     {
         const vectorField& c = mesh().cellCentres();
         const scalarField x((c - origins_[i]) & directions_[i]);
-        forceCoeff = max(forceCoeff, scale_->value(x));
+        scale = max(scale, scale_->value(x));
     }
 
-    forceCoeff *= lambda_.value();
-
     // Write out the force coefficient for debugging
     if (debug && mesh().time().writeTime())
     {
-        volScalarField vForceCoeff
-        (
-            IOobject
-            (
-                typedName("forceCoeff"),
-                mesh().time().timeName(),
-                mesh()
-            ),
-            mesh(),
-            lambda_.dimensions(),
-            zeroGradientFvPatchField<scalar>::typeName
-        );
-        vForceCoeff.primitiveFieldRef() = forceCoeff;
-        vForceCoeff.correctBoundaryConditions();
-        vForceCoeff.write();
+        tscale->write();
+    }
+
+    return tscale;
+}
+
+
+Foam::tmp<Foam::volScalarField::Internal> Foam::fv::forcing::forceCoeff
+(
+    const volScalarField::Internal& scale
+) const
+{
+    tmp<volScalarField::Internal> tforceCoeff
+    (
+        volScalarField::Internal::New(typedName("forceCoeff"), lambda_*scale)
+    );
+
+    // Write out the force coefficient for debugging
+    if (debug && mesh().time().writeTime())
+    {
+        tforceCoeff->write();
     }
 
     return tforceCoeff;
 }
 
 
+Foam::tmp<Foam::volScalarField::Internal> Foam::fv::forcing::forceCoeff() const
+{
+    return forceCoeff(scale());
+}
+
+
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 Foam::fv::forcing::forcing

src/waves/fvModels/forcing/forcing.H

@@ -80,6 +80,15 @@ protected:
         //- Non-virtual read
         void readCoeffs();
 
+        //- Return the scale distribution
+        tmp<volScalarField::Internal> scale() const;
+
+        //- Return the force coefficient given the scale distribution
+        tmp<volScalarField::Internal> forceCoeff
+        (
+            const volScalarField::Internal& scale
+        ) const;
+
         //- Return the force coefficient
         tmp<volScalarField::Internal> forceCoeff() const;
 

src/waves/fvModels/waveForcing/waveForcing.C

@@ -29,6 +29,9 @@ License
 #include "fvmSup.H"
 #include "addToRunTimeSelectionTable.H"
 
+#include "fvcDdt.H"
+#include "fvcDiv.H"
+
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
 namespace Foam
@@ -56,7 +59,7 @@ Foam::fv::waveForcing::waveForcing
     liquidPhaseName_(coeffs().lookup<word>("liquidPhase")),
     alphaName_(IOobject::groupName("alpha", liquidPhaseName_)),
     UName_(coeffs().lookupOrDefault<word>("U", "U")),
-    forceCoeff_(this->forceCoeff())
+    scale_(this->scale())
 {}
 
 
@@ -76,8 +79,10 @@ void Foam::fv::waveForcing::addSup
 {
     if (fieldName == alphaName_)
     {
-        eqn -= fvm::Sp(forceCoeff_(), eqn.psi());
-        eqn += forceCoeff_()*alphaWaves_();
+        const volScalarField::Internal forceCoeff(this->forceCoeff(scale_));
+
+        eqn -= fvm::Sp(forceCoeff, eqn.psi());
+        eqn += forceCoeff*alphaWaves_();
     }
 }
 
@@ -91,34 +96,45 @@ void Foam::fv::waveForcing::addSup
 {
     if (fieldName == UName_)
     {
-        eqn -= fvm::Sp(rho*forceCoeff_(), eqn.psi());
-        eqn += rho*forceCoeff_()*Uwaves_();
+        const volScalarField::Internal forceCoeff(rho*this->forceCoeff(scale_));
+
+        eqn -= fvm::Sp(forceCoeff, eqn.psi());
+        eqn += forceCoeff*Uwaves_();
+
+        const surfaceScalarField& rhoPhi =
+            mesh().lookupObject<surfaceScalarField>("rhoPhi");
+
+        eqn += fvm::Sp
+        (
+            scale()*(fvc::ddt(rho)()() + fvc::div(rhoPhi)()()),
+            eqn.psi()
+        );
     }
 }
 
 
 bool Foam::fv::waveForcing::movePoints()
 {
-    forceCoeff_ = this->forceCoeff();
+    scale_ = this->scale();
     return true;
 }
 
 
 void Foam::fv::waveForcing::topoChange(const polyTopoChangeMap&)
 {
-    forceCoeff_ = this->forceCoeff();
+    scale_ = this->scale();
 }
 
 
 void Foam::fv::waveForcing::mapMesh(const polyMeshMap& map)
 {
-    forceCoeff_ = this->forceCoeff();
+    scale_ = this->scale();
 }
 
 
 void Foam::fv::waveForcing::distribute(const polyDistributionMap&)
 {
-    forceCoeff_ = this->forceCoeff();
+    scale_ = this->scale();
 }
 
 

src/waves/fvModels/waveForcing/waveForcing.H

@@ -115,8 +115,8 @@ class waveForcing
         //- Velocity field calculated from the current wave form
         tmp<volVectorField::Internal> Uwaves_;
 
-        //- Forcing coefficient field
-        tmp<volScalarField::Internal> forceCoeff_;
+        //- Forcing coefficient scale field
+        tmp<volScalarField::Internal> scale_;
 
 
 public: