From 974624766fedeb1ae199ad982035cc0053b2d96a Mon Sep 17 00:00:00 2001
From: Henry Weller <http://cfd.direct>
Date: Tue, 7 Feb 2017 19:02:30 +0000
Subject: [PATCH] porosityModels::solidification: New porosity model to
 simulate solidification

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 = - \rho D(T) U
        \f]

    where
    \vartable
        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
            );

            // 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/Make/files                   |   1 +
 .../solidification/solidification.C           | 163 +++++++++++++
 .../solidification/solidification.H           | 214 ++++++++++++++++++
 .../solidification/solidificationTemplates.C  |  83 +++++++
 4 files changed, 461 insertions(+)
 create mode 100644 src/finiteVolume/cfdTools/general/porosityModel/solidification/solidification.C
 create mode 100644 src/finiteVolume/cfdTools/general/porosityModel/solidification/solidification.H
 create mode 100644 src/finiteVolume/cfdTools/general/porosityModel/solidification/solidificationTemplates.C

diff --git a/src/finiteVolume/Make/files b/src/finiteVolume/Make/files
index d79bbff4e82..2d7fbce1304 100644
--- a/src/finiteVolume/Make/files
+++ b/src/finiteVolume/Make/files
@@ -417,6 +417,7 @@ $(porosity)/porosityModel/IOporosityModelList.C
 $(porosity)/DarcyForchheimer/DarcyForchheimer.C
 $(porosity)/fixedCoeff/fixedCoeff.C
 $(porosity)/powerLaw/powerLaw.C
+$(porosity)/solidification/solidification.C
 
 MRF = $(general)/MRF
 $(MRF)/MRFZone.C
diff --git a/src/finiteVolume/cfdTools/general/porosityModel/solidification/solidification.C b/src/finiteVolume/cfdTools/general/porosityModel/solidification/solidification.C
new file mode 100644
index 00000000000..ac60cf0d04c
--- /dev/null
+++ b/src/finiteVolume/cfdTools/general/porosityModel/solidification/solidification.C
@@ -0,0 +1,163 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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 "addToRunTimeSelectionTable.H"
+#include "solidification.H"
+#include "geometricOneField.H"
+#include "fvMatrices.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    namespace porosityModels
+    {
+        defineTypeNameAndDebug(solidification, 0);
+        addToRunTimeSelectionTable(porosityModel, solidification, mesh);
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::porosityModels::solidification::solidification
+(
+    const word& name,
+    const word& modelType,
+    const fvMesh& mesh,
+    const dictionary& dict,
+    const word& cellZoneName
+)
+:
+    porosityModel(name, modelType, mesh, dict, cellZoneName),
+    TName_(coeffs_.lookupOrDefault<word>("T", "T")),
+    rhoName_(coeffs_.lookupOrDefault<word>("rho", "rho")),
+    D_(Function1<scalar>::New("D", coeffs_))
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::porosityModels::solidification::~solidification()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void Foam::porosityModels::solidification::calcTransformModelData()
+{}
+
+
+void Foam::porosityModels::solidification::calcForce
+(
+    const volVectorField& U,
+    const volScalarField& rho,
+    const volScalarField& mu,
+    vectorField& force
+) const
+{
+    scalarField Udiag(U.size(), 0.0);
+    const scalarField& V = mesh_.V();
+
+    apply(Udiag, V, rho, U);
+
+    force = Udiag*U;
+}
+
+
+void Foam::porosityModels::solidification::correct
+(
+    fvVectorMatrix& UEqn
+) const
+{
+    const volVectorField& U = UEqn.psi();
+    const scalarField& V = mesh_.V();
+    scalarField& Udiag = UEqn.diag();
+
+    if (UEqn.dimensions() == dimForce)
+    {
+        const volScalarField& rho = mesh_.lookupObject<volScalarField>
+        (
+            IOobject::groupName(rhoName_, U.group())
+        );
+
+        apply(Udiag, V, rho, U);
+    }
+    else
+    {
+        apply(Udiag, V, geometricOneField(), U);
+    }
+}
+
+
+void Foam::porosityModels::solidification::correct
+(
+    fvVectorMatrix& UEqn,
+    const volScalarField& rho,
+    const volScalarField& mu
+) const
+{
+    const volVectorField& U = UEqn.psi();
+    const scalarField& V = mesh_.V();
+    scalarField& Udiag = UEqn.diag();
+
+    apply(Udiag, V, rho, U);
+}
+
+
+void Foam::porosityModels::solidification::correct
+(
+    const fvVectorMatrix& UEqn,
+    volTensorField& AU
+) const
+{
+    const volVectorField& U = UEqn.psi();
+
+    if (UEqn.dimensions() == dimForce)
+    {
+        const volScalarField& rho = mesh_.lookupObject<volScalarField>
+        (
+            IOobject::groupName(rhoName_, U.group())
+        );
+
+        apply(AU, rho, U);
+    }
+    else
+    {
+        apply(AU, geometricOneField(), U);
+    }
+}
+
+
+bool Foam::porosityModels::solidification::writeData(Ostream& os) const
+{
+    os  << indent << name_ << endl;
+    dict_.write(os);
+
+    return true;
+}
+
+
+// ************************************************************************* //
diff --git a/src/finiteVolume/cfdTools/general/porosityModel/solidification/solidification.H b/src/finiteVolume/cfdTools/general/porosityModel/solidification/solidification.H
new file mode 100644
index 00000000000..7f86849f64d
--- /dev/null
+++ b/src/finiteVolume/cfdTools/general/porosityModel/solidification/solidification.H
@@ -0,0 +1,214 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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::porosityModels::solidification
+
+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 = - \rho D(T) U
+        \f]
+
+    where
+    \vartable
+        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
+            );
+
+            // 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
+
+SourceFiles
+    solidification.C
+    solidificationTemplates.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef solidification_H
+#define solidification_H
+
+#include "porosityModel.H"
+#include "Function1.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace porosityModels
+{
+
+/*---------------------------------------------------------------------------*\
+                      Class solidification Declaration
+\*---------------------------------------------------------------------------*/
+
+class solidification
+:
+    public porosityModel
+{
+    // Private data
+
+        //- Name of temperature field, default = "T"
+        word TName_;
+
+        //- Name of density field, default = "rho"
+        word rhoName_;
+
+        //- User-defined drag-coefficient as function of temperature
+        autoPtr<Function1<scalar>> D_;
+
+
+    // Private Member Functions
+
+        //- Apply resistance
+        template<class RhoFieldType>
+        void apply
+        (
+            scalarField& Udiag,
+            const scalarField& V,
+            const RhoFieldType& rho,
+            const volVectorField& U
+        ) const;
+
+        //- Apply resistance
+        template<class RhoFieldType>
+        void apply
+        (
+            tensorField& AU,
+            const RhoFieldType& rho,
+            const volVectorField& U
+        ) const;
+
+        //- Disallow default bitwise copy construct
+        solidification(const solidification&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const solidification&);
+
+
+public:
+
+    //- Runtime type information
+    TypeName("solidification");
+
+    //- Constructor
+    solidification
+    (
+        const word& name,
+        const word& modelType,
+        const fvMesh& mesh,
+        const dictionary& dict,
+        const word& cellZoneName
+    );
+
+    //- Destructor
+    virtual ~solidification();
+
+
+    // Member Functions
+
+        //- Transform the model data wrt mesh changes
+        virtual void calcTransformModelData();
+
+        //- Calculate the porosity force
+        virtual void calcForce
+        (
+            const volVectorField& U,
+            const volScalarField& rho,
+            const volScalarField& mu,
+            vectorField& force
+        ) const;
+
+        //- Add resistance
+        virtual void correct(fvVectorMatrix& UEqn) const;
+
+        //- Add resistance
+        virtual void correct
+        (
+            fvVectorMatrix& UEqn,
+            const volScalarField& rho,
+            const volScalarField& mu
+        ) const;
+
+        //- Add resistance
+        virtual void correct
+        (
+            const fvVectorMatrix& UEqn,
+            volTensorField& AU
+        ) const;
+
+
+    // I-O
+
+        //- Write
+        bool writeData(Ostream& os) const;
+};
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace porosityModels
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+    #include "solidificationTemplates.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/finiteVolume/cfdTools/general/porosityModel/solidification/solidificationTemplates.C b/src/finiteVolume/cfdTools/general/porosityModel/solidification/solidificationTemplates.C
new file mode 100644
index 00000000000..71b290af178
--- /dev/null
+++ b/src/finiteVolume/cfdTools/general/porosityModel/solidification/solidificationTemplates.C
@@ -0,0 +1,83 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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 "volFields.H"
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+template<class RhoFieldType>
+void Foam::porosityModels::solidification::apply
+(
+    scalarField& Udiag,
+    const scalarField& V,
+    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];
+            Udiag[celli] += V[celli]*rho[celli]*D_->value(T[celli]);
+        }
+    }
+}
+
+
+template<class RhoFieldType>
+void Foam::porosityModels::solidification::apply
+(
+    tensorField& AU,
+    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*rho[celli]*D_->value(T[celli]);
+        }
+    }
+}
+
+
+// ************************************************************************* //
-- 
GitLab