porosityModels::solidification: Added optional phase-fraction for VoF solvers etc.

Description
    Simple solidification porosity model

    This is a simple approximation to solidification where the solid phase
    is represented as a porous blockage with the drag-coefficient evaluated from

        \f[
            S = - \alpha \rho D(T) U
        \f]

    where
    \vartable
        \alpha  | Optional phase-fraction of solidifying phase
        D(T)    | User-defined drag-coefficient as function of temperature
    \endvartable

    Note that the latent heat of solidification is not included and the
    temperature is unchanged by the modelled change of phase.

    Example of the solidification model specification:
    \verbatim
        type            solidification;

        solidificationCoeffs
        {
            // Solidify between 330K and 330.5K
            D table
            (
                (330.0     10000) // Solid below 330K
                (330.5     0)     // Liquid above 330.5K
            );

            // Optional phase-fraction of solidifying phase
            alpha alpha.liquid;

            // Solidification porosity is isotropic
            // use the global coordinate system
            coordinateSystem
            {
                type    cartesian;
                origin  (0 0 0);
                coordinateRotation
                {
                    type    axesRotation;
                    e1      (1 0 0);
                    e2      (0 1 0);
                }
            }
        }
    \endverbatim

src/finiteVolume/cfdTools/general/porosityModel/solidification/solidification.C

@@ -53,6 +53,7 @@ Foam::porosityModels::solidification::solidification
 :
     porosityModel(name, modelType, mesh, dict, cellZoneName),
     TName_(coeffs_.lookupOrDefault<word>("T", "T")),
+    alphaName_(coeffs_.lookupOrDefault<word>("alpha", "none")),
     rhoName_(coeffs_.lookupOrDefault<word>("rho", "rho")),
     D_(Function1<scalar>::New("D", coeffs_))
 {}

src/finiteVolume/cfdTools/general/porosityModel/solidification/solidification.H

@@ -31,12 +31,13 @@ Description
     is represented as a porous blockage with the drag-coefficient evaluated from
 
         \f[
-            S = - \rho D(T) U
+            S = - \alpha \rho D(T) U
         \f]
 
     where
     \vartable
-        D(T) | User-defined drag-coefficient as function of temperature
+        \alpha  | Optional phase-fraction of solidifying phase
+        D(T)    | User-defined drag-coefficient as function of temperature
     \endvartable
 
     Note that the latent heat of solidification is not included and the
@@ -55,6 +56,9 @@ Description
                 (330.5     0)     // Liquid above 330.5K
             );
 
+            // Optional phase-fraction of solidifying phase
+            alpha alpha.liquid;
+
             // Solidification porosity is isotropic
             // use the global coordinate system
             coordinateSystem
@@ -103,6 +107,9 @@ class solidification
         //- Name of temperature field, default = "T"
         word TName_;
 
+        //- Name of optional phase-fraction field, default = "none"
+        word alphaName_;
+
         //- Name of density field, default = "rho"
         word rhoName_;
 
@@ -112,6 +119,27 @@ class solidification
 
     // Private Member Functions
 
+        //- Apply resistance
+        template<class AlphaFieldType, class RhoFieldType>
+        void apply
+        (
+            scalarField& Udiag,
+            const scalarField& V,
+            const AlphaFieldType& alpha,
+            const RhoFieldType& rho,
+            const volVectorField& U
+        ) const;
+
+        //- Apply resistance
+        template<class AlphaFieldType, class RhoFieldType>
+        void apply
+        (
+            tensorField& AU,
+            const AlphaFieldType& alpha,
+            const RhoFieldType& rho,
+            const volVectorField& U
+        ) const;
+
         //- Apply resistance
         template<class RhoFieldType>
         void apply

src/finiteVolume/cfdTools/general/porosityModel/solidification/solidificationTemplates.C

@@ -24,14 +24,16 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "volFields.H"
+#include "geometricOneField.H"
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
-template<class RhoFieldType>
+template<class AlphaFieldType, class RhoFieldType>
 void Foam::porosityModels::solidification::apply
 (
     scalarField& Udiag,
     const scalarField& V,
+    const AlphaFieldType& alpha,
     const RhoFieldType& rho,
     const volVectorField& U
 ) const
@@ -48,12 +50,66 @@ void Foam::porosityModels::solidification::apply
         forAll(cells, i)
         {
             const label celli = cells[i];
-            Udiag[celli] += V[celli]*rho[celli]*D_->value(T[celli]);
+            Udiag[celli] +=
+                V[celli]*alpha[celli]*rho[celli]*D_->value(T[celli]);
         }
     }
 }
 
 
+template<class AlphaFieldType, class RhoFieldType>
+void Foam::porosityModels::solidification::apply
+(
+    tensorField& AU,
+    const AlphaFieldType& alpha,
+    const RhoFieldType& rho,
+    const volVectorField& U
+) const
+{
+    const volScalarField& T = mesh_.lookupObject<volScalarField>
+    (
+        IOobject::groupName(TName_, U.group())
+    );
+
+    forAll(cellZoneIDs_, zoneI)
+    {
+        const labelList& cells = mesh_.cellZones()[cellZoneIDs_[zoneI]];
+
+        forAll(cells, i)
+        {
+            const label celli = cells[i];
+            AU[celli] +=
+                tensor::I*alpha[celli]*rho[celli]*D_->value(T[celli]);
+        }
+    }
+}
+
+
+template<class RhoFieldType>
+void Foam::porosityModels::solidification::apply
+(
+    scalarField& Udiag,
+    const scalarField& V,
+    const RhoFieldType& rho,
+    const volVectorField& U
+) const
+{
+    if (alphaName_ == "none")
+    {
+        return apply(Udiag, V, geometricOneField(), rho, U);
+    }
+    else
+    {
+        const volScalarField& alpha = mesh_.lookupObject<volScalarField>
+        (
+            IOobject::groupName(alphaName_, U.group())
+        );
+
+        return apply(Udiag, V, alpha, rho, U);
+    }
+}
+
+
 template<class RhoFieldType>
 void Foam::porosityModels::solidification::apply
 (
@@ -62,20 +118,18 @@ void Foam::porosityModels::solidification::apply
     const volVectorField& U
 ) const
 {
-    const volScalarField& T = mesh_.lookupObject<volScalarField>
+    if (alphaName_ == "none")
-    (
-        IOobject::groupName(TName_, U.group())
-    );
-
-    forAll(cellZoneIDs_, zoneI)
     {
-        const labelList& cells = mesh_.cellZones()[cellZoneIDs_[zoneI]];
+        return apply(AU, geometricOneField(), rho, U);
+    }
+    else
+    {
+        const volScalarField& alpha = mesh_.lookupObject<volScalarField>
+        (
+            IOobject::groupName(alphaName_, U.group())
+        );
 
-        forAll(cells, i)
+        return apply(AU, alpha, rho, U);
-        {
-            const label celli = cells[i];
-            AU[celli] += tensor::I*rho[celli]*D_->value(T[celli]);
-        }
     }
 }