ENH: fvOptions::PhaseLimitStabilization: New fvOption to stabilize phase transport equations

in the limit of the phase fraction -> 0

Can be used with the Maxwell non-Newtonian laminar stress model when used in
multiphase solvers.

src/TurbulenceModels/turbulenceModels/laminar/Maxwell/Maxwell.C

@@ -25,6 +25,7 @@ License
 
 #include "Maxwell.H"
 #include "fvOptions.H"
+#include "uniformDimensionedFields.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -205,7 +206,17 @@ void Maxwell<BasicTurbulenceModel>::correct()
 
     tmp<volTensorField> tgradU(fvc::grad(U));
     const volTensorField& gradU = tgradU();
-    dimensionedScalar rLambda = 1.0/(lambda_);
+
+    uniformDimensionedScalarField rLambda
+    (
+        IOobject
+        (
+            IOobject::groupName("rLambda", this->alphaRhoPhi_.group()),
+            this->runTime_.constant(),
+            this->mesh_
+        ),
+        1.0/(lambda_)
+    );
 
     // Note sigma is positive on lhs of momentum eqn
     volSymmTensorField P
@@ -220,7 +231,7 @@ void Maxwell<BasicTurbulenceModel>::correct()
         fvm::ddt(alpha, rho, sigma)
       + fvm::div(alphaRhoPhi, sigma)
       + fvm::Sp(alpha*rho*rLambda, sigma)
-      ==
+     ==
         alpha*rho*P
       + fvOptions(alpha, rho, sigma)
     );

src/fvOptions/Make/files

@@ -27,6 +27,7 @@ $(derivedSources)/jouleHeatingSource/jouleHeatingSourceIO.C
 $(derivedSources)/meanVelocityForce/meanVelocityForce.C
 $(derivedSources)/meanVelocityForce/meanVelocityForceIO.C
 $(derivedSources)/meanVelocityForce/patchMeanVelocityForce/patchMeanVelocityForce.C
+$(derivedSources)/phaseLimitStabilization/phaseLimitStabilization.C
 $(derivedSources)/radialActuationDiskSource/radialActuationDiskSource.C
 $(derivedSources)/rotorDiskSource/rotorDiskSource.C
 $(derivedSources)/rotorDiskSource/bladeModel/bladeModel.C

src/fvOptions/sources/derived/phaseLimitStabilization/PhaseLimitStabilization.C

@@ -0,0 +1,88 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2017 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "PhaseLimitStabilization.H"
+#include "fvMatrices.H"
+#include "fvmSup.H"
+#include "uniformDimensionedFields.H"
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+template<class Type>
+Foam::fv::PhaseLimitStabilization<Type>::PhaseLimitStabilization
+(
+    const word& name,
+    const word& modelType,
+    const dictionary& dict,
+    const fvMesh& mesh
+)
+:
+    option(name, modelType, dict, mesh),
+    fieldName_(coeffs_.lookup("field")),
+    rateName_(coeffs_.lookup("rate")),
+    residualAlpha_(readScalar(coeffs_.lookup("residualAlpha")))
+{
+    fieldNames_.setSize(1, fieldName_);
+    applied_.setSize(1, false);
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class Type>
+void Foam::fv::PhaseLimitStabilization<Type>::addSup
+(
+    const volScalarField& alpha,
+    const volScalarField& rho,
+    fvMatrix<Type>& eqn,
+    const label fieldi
+)
+{
+    const GeometricField<Type, fvPatchField, volMesh>& psi = eqn.psi();
+
+    uniformDimensionedScalarField& rate =
+        mesh_.lookupObjectRef<uniformDimensionedScalarField>(rateName_);
+
+    eqn -= fvm::Sp(max(residualAlpha_ - alpha, scalar(0))*rho*rate, psi);
+}
+
+
+template<class Type>
+bool Foam::fv::PhaseLimitStabilization<Type>::read(const dictionary& dict)
+{
+    if (option::read(dict))
+    {
+        coeffs_.lookup("residualAlpha") >> residualAlpha_;
+
+        return true;
+    }
+    else
+    {
+        return false;
+    }
+}
+
+
+// ************************************************************************* //

src/fvOptions/sources/derived/phaseLimitStabilization/PhaseLimitStabilization.H

@@ -0,0 +1,153 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2017 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::fv::PhaseLimitStabilization
+
+Description
+    Stabilization source for phase transport equations
+
+    Applies an implicit source to the phase transport equation for the specified
+    \c field when the phase volume fraction is below \c residualAlpha.  The
+    stabilization rate is provided by the registered
+    uniformDimensionedScalarField \c rate, which could be extended to also
+    support volScalarField and volScalarField::Internal field types.  The \c
+    field is currently stabilized towards zero in the limit of the phase volume
+    fraction approaching zero but this could be extended to support a
+    specified value or a value or field looked-up from the database.
+
+Usage
+    Example usage:
+    \verbatim
+    stabilization
+    {
+        type            symmTensorPhaseLimitStabilization;
+
+        field           sigma.liquid;
+        rate            rLambda.liquid;
+        residualAlpha   1e-3;
+    }
+    \endverbatim
+
+SourceFiles
+    PhaseLimitStabilization.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef PhaseLimitStabilization_H
+#define PhaseLimitStabilization_H
+
+#include "fvOption.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace fv
+{
+
+/*---------------------------------------------------------------------------*\
+                Class PhaseLimitStabilization Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class Type>
+class PhaseLimitStabilization
+:
+    public option
+{
+    // Private data
+
+        //- Field name
+        word fieldName_;
+
+        //- Rate field name
+        word rateName_;
+
+        //- Residual alpha value below which stabilization is applied
+        scalar residualAlpha_;
+
+
+    // Private Member Functions
+
+        //- Disallow default bitwise copy construct
+        PhaseLimitStabilization(const PhaseLimitStabilization&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const PhaseLimitStabilization&);
+
+
+public:
+
+    //- Runtime type information
+    TypeName("PhaseLimitStabilization");
+
+
+    // Constructors
+
+        //- Construct from components
+        PhaseLimitStabilization
+        (
+            const word& name,
+            const word& modelType,
+            const dictionary& dict,
+            const fvMesh& mesh
+        );
+
+
+    //- Destructor
+    virtual ~PhaseLimitStabilization()
+    {}
+
+
+    // Member Functions
+
+        //- Source term to compressible phase equation
+        virtual void addSup
+        (
+            const volScalarField& alpha,
+            const volScalarField& rho,
+            fvMatrix<Type>& eqn,
+            const label fieldi
+        );
+
+        //- Read dictionary
+        virtual bool read(const dictionary& dict);
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace fv
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+    #include "PhaseLimitStabilization.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/fvOptions/sources/derived/phaseLimitStabilization/phaseLimitStabilization.C

@@ -0,0 +1,38 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2017 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "makeFvOption.H"
+#include "PhaseLimitStabilization.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+makeFvOption(PhaseLimitStabilization, scalar);
+makeFvOption(PhaseLimitStabilization, vector);
+makeFvOption(PhaseLimitStabilization, sphericalTensor);
+makeFvOption(PhaseLimitStabilization, symmTensor);
+makeFvOption(PhaseLimitStabilization, tensor);
+
+
+// ************************************************************************* //