fv::massSource: Support mass sinks

If a negative mass flow rate is specified, the mass source fvModel will
now remove mass by adding implicit sources to the transport equations.
Properties are thereby removed at their current value. This is stable,
and is analogous to a zero-gradient outlet boundary condition.

src/fvModels/derived/massSource/massSource.C

@@ -98,16 +98,28 @@ void Foam::fv::massSource::addGeneralSupType
     const word& fieldName
 ) const
 {
-    const scalar massFlowRate = this->massFlowRate();
-    const Type value =
-        fieldValues_[fieldName]->value<Type>(mesh().time().userTimeValue());
-
     const labelUList cells = set_.cells();
 
-    forAll(cells, i)
+    const scalar massFlowRate = this->massFlowRate();
+
+    if (massFlowRate > 0)
     {
-        eqn.source()[cells[i]] -=
-            mesh().V()[cells[i]]/set_.V()*massFlowRate*value;
+        const Type value =
+            fieldValues_[fieldName]->value<Type>(mesh().time().userTimeValue());
+
+        forAll(cells, i)
+        {
+            eqn.source()[cells[i]] -=
+                mesh().V()[cells[i]]/set_.V()*massFlowRate*value;
+        }
+    }
+    else
+    {
+        forAll(cells, i)
+        {
+            eqn.diag()[cells[i]] +=
+                mesh().V()[cells[i]]/set_.V()*massFlowRate;
+        }
     }
 }
 
@@ -143,31 +155,52 @@ void Foam::fv::massSource::addSupType
     }
     else if (fieldName == heName_ && fieldValues_.found(TName_))
     {
-        if (fieldValues_.found(heName_))
-        {
-            WarningInFunction
-                << "Source " << name() << " defined for both field " << heName_
-                << " and " << TName_ << ". Only one of these should be present."
-                << endl;
-        }
-
         const scalar massFlowRate = this->massFlowRate();
-        const scalar T =
-            fieldValues_[TName_]->value<scalar>(mesh().time().userTimeValue());
-        const basicThermo& thermo =
-            mesh().lookupObject<basicThermo>
-            (
-                IOobject::groupName(physicalProperties::typeName, phaseName_)
-            );
-        const scalarField hs
-        (
-            thermo.hs(scalarField(cells.size(), T), cells)
-        );
 
-        forAll(cells, i)
+        if (massFlowRate > 0)
         {
-            eqn.source()[cells[i]] -=
-                mesh().V()[cells[i]]/set_.V()*massFlowRate*hs[i];
+            if (fieldValues_.found(heName_))
+            {
+                WarningInFunction
+                    << "Source " << name() << " defined for both field "
+                    << heName_ << " and " << TName_
+                    << ". Only one of these should be present." << endl;
+            }
+
+            const basicThermo& thermo =
+                mesh().lookupObject<basicThermo>
+                (
+                    IOobject::groupName
+                    (
+                        physicalProperties::typeName,
+                        phaseName_
+                    )
+                );
+
+            const scalar T =
+                fieldValues_[TName_]->value<scalar>
+                (
+                    mesh().time().userTimeValue()
+                );
+
+            const scalarField hs
+            (
+                thermo.hs(scalarField(cells.size(), T), cells)
+            );
+
+            forAll(cells, i)
+            {
+                eqn.source()[cells[i]] -=
+                    mesh().V()[cells[i]]/set_.V()*massFlowRate*hs[i];
+            }
+        }
+        else
+        {
+            forAll(cells, i)
+            {
+                eqn.diag()[cells[i]] +=
+                    mesh().V()[cells[i]]/set_.V()*massFlowRate;
+            }
         }
     }
     else
@@ -235,6 +268,7 @@ bool Foam::fv::massSource::addsSupToField(const word& fieldName) const
     if
     (
         isThisPhase
+     && massFlowRate() > 0
      && !(fieldName == rhoName_)
      && !(fieldName == heName_ && fieldValues_.found(TName_))
      && !fieldValues_.found(fieldName)

src/fvModels/derived/massSource/massSource.H

@@ -28,6 +28,10 @@ Description
     This fvModel applies a mass source to the continuity equation and to all
     field equations.
 
+    If the mass flow rate is positive then user-supplied fixed property values
+    are introduced to the field equations. If the mass flow rate is negative
+    then properties are removed at their current value.
+
 Usage
     Example usage:
     \verbatim
@@ -50,10 +54,10 @@ Usage
     }
     \endverbatim
 
-    Values should be provided for all solved for fields. Warnings will be
-    issued if values are not provided for fields for which transport equations
-    are solved. Warnings will also be issued if values are provided for fields
-    which are not solved for.
+    If the mass flow rate is positive then values should be provided for all
+    solved for fields. Warnings will be issued if values are not provided for
+    fields for which transport equations are solved. Warnings will also be
+    issued if values are provided for fields which are not solved for.
 
 SourceFiles
     massSource.C

test/multiphase/multiphaseEuler/populationBalance/drift/constant/fvModels

@@ -26,9 +26,7 @@ waterSink
     rho             rho.water;
 
     fieldValues
-    {
-        U.water     (0 0 0);
-    }
+    {}
 }
 
 // ************************************************************************* //