compressibleInterFilmFoam::VoFSolidificationMeltingSource: New VoF solidification fvOption

applications/solvers/multiphase/compressibleInterFoam/compressibleInterFilmFoam/Make/files

@@ -1,4 +1,6 @@
 VoFPatchTransfer/VoFPatchTransfer.C
+VoFSolidificationMeltingSource/VoFSolidificationMeltingSource.C
+VoFSolidificationMeltingSource/VoFSolidificationMeltingSourceIO.C
 compressibleInterFilmFoam.C
 
 EXE = $(FOAM_APPBIN)/compressibleInterFilmFoam

applications/solvers/multiphase/compressibleInterFoam/compressibleInterFilmFoam/Make/options

@@ -21,6 +21,7 @@ EXE_INC = \
     -I$(LIB_SRC)/dynamicFvMesh/lnInclude \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/specie/lnInclude \
+    -I$(LIB_SRC)/fvOptions/lnInclude
 
 EXE_LIBS = \
     -ltwoPhaseMixtureThermo \

applications/solvers/multiphase/compressibleInterFoam/compressibleInterFilmFoam/VoFSolidificationMeltingSource/VoFSolidificationMeltingSource.C

@@ -0,0 +1,215 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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 "VoFSolidificationMeltingSource.H"
+#include "twoPhaseMixtureThermo.H"
+#include "zeroGradientFvPatchFields.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * Static Member Functions * * * * * * * * * * * * //
+
+namespace Foam
+{
+    namespace fv
+    {
+        defineTypeNameAndDebug(VoFSolidificationMeltingSource, 0);
+
+        addToRunTimeSelectionTable
+        (
+            option,
+            VoFSolidificationMeltingSource,
+            dictionary
+        );
+    }
+}
+
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+void Foam::fv::VoFSolidificationMeltingSource::update()
+{
+    if (curTimeIndex_ == mesh_.time().timeIndex())
+    {
+        return;
+    }
+
+    if (debug)
+    {
+        Info<< type() << ": " << name_
+            << " - updating solid phase fraction" << endl;
+    }
+
+    alphaSolid_.oldTime();
+
+    const twoPhaseMixtureThermo& thermo
+    (
+        mesh_.lookupObject<twoPhaseMixtureThermo>
+        (
+            twoPhaseMixtureThermo::dictName
+        )
+    );
+
+    const volScalarField& TVoF = thermo.thermo1().T();
+    const volScalarField CpVoF(thermo.thermo1().Cp());
+    const volScalarField& alphaVoF = thermo.alpha1();
+
+    forAll(cells_, i)
+    {
+        const label celli = cells_[i];
+
+        alphaSolid_[celli] = min
+        (
+            relax_*alphaVoF[celli]*alphaSolidT_->value(TVoF[celli])
+          + (1 - relax_)*alphaSolid_[celli],
+            alphaVoF[celli]
+        );
+    }
+
+    alphaSolid_.correctBoundaryConditions();
+
+    curTimeIndex_ = mesh_.time().timeIndex();
+}
+
+
+Foam::word Foam::fv::VoFSolidificationMeltingSource::alphaSolidName() const
+{
+    const twoPhaseMixtureThermo& thermo
+    (
+        mesh_.lookupObject<twoPhaseMixtureThermo>
+        (
+            twoPhaseMixtureThermo::dictName
+        )
+    );
+
+    const volScalarField& alphaVoF = thermo.alpha1();
+
+    return IOobject::groupName(alphaVoF.name(), "solid");
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::fv::VoFSolidificationMeltingSource::VoFSolidificationMeltingSource
+(
+    const word& sourceName,
+    const word& modelType,
+    const dictionary& dict,
+    const fvMesh& mesh
+)
+:
+    cellSetOption(sourceName, modelType, dict, mesh),
+    alphaSolidT_(Function1<scalar>::New("alphaSolidT", coeffs_)),
+    L_("L", dimEnergy/dimMass, coeffs_),
+    relax_(coeffs_.lookupOrDefault("relax", 0.9)),
+    Cu_(coeffs_.lookupOrDefault<scalar>("Cu", 100000)),
+    q_(coeffs_.lookupOrDefault("q", 0.001)),
+    alphaSolid_
+    (
+        IOobject
+        (
+            alphaSolidName(),
+            mesh.time().timeName(),
+            mesh,
+            IOobject::READ_IF_PRESENT,
+            IOobject::AUTO_WRITE
+        ),
+        mesh,
+        dimensionedScalar("alpha1", dimless, 0.0),
+        zeroGradientFvPatchScalarField::typeName
+    ),
+    curTimeIndex_(-1)
+{
+    fieldNames_.setSize(2);
+    fieldNames_[0] = "U";
+    fieldNames_[1] = "T";
+    applied_.setSize(fieldNames_.size(), false);
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void Foam::fv::VoFSolidificationMeltingSource::addSup
+(
+    fvMatrix<scalar>& eqn,
+    const label fieldi
+)
+{
+    apply(geometricOneField(), eqn);
+}
+
+
+void Foam::fv::VoFSolidificationMeltingSource::addSup
+(
+    const volScalarField& rho,
+    fvMatrix<scalar>& eqn,
+    const label fieldi
+)
+{
+    apply(rho, eqn);
+}
+
+
+void Foam::fv::VoFSolidificationMeltingSource::addSup
+(
+    fvMatrix<vector>& eqn,
+    const label fieldi
+)
+{
+    if (debug)
+    {
+        Info<< type() << ": applying source to " << eqn.psi().name() << endl;
+    }
+
+    update();
+
+    scalarField& Sp = eqn.diag();
+    const scalarField& V = mesh_.V();
+
+    forAll(cells_, i)
+    {
+        const label celli = cells_[i];
+        const scalar Vc = V[celli];
+        const scalar alphaFluid = 1 - alphaSolid_[celli];
+
+        const scalar S = Cu_*sqr(1 - alphaFluid)/(pow3(alphaFluid) + q_);
+
+        Sp[celli] -= Vc*S;
+    }
+}
+
+
+void Foam::fv::VoFSolidificationMeltingSource::addSup
+(
+    const volScalarField& rho,
+    fvMatrix<vector>& eqn,
+    const label fieldi
+)
+{
+    // Momentum source uses a Boussinesq approximation - redirect
+    addSup(eqn, fieldi);
+}
+
+
+// ************************************************************************* //

applications/solvers/multiphase/compressibleInterFoam/compressibleInterFilmFoam/VoFSolidificationMeltingSource/VoFSolidificationMeltingSource.H

@@ -0,0 +1,215 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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::VoFSolidificationMeltingSource
+
+Description
+    Solidification and melting model for VoF simulations.
+
+    The presence of the solid phase in the flow field is incorporated into the
+    model as a momentum porosity contribution; the energy associated with the
+    phase change is added as an enthalpy contribution.  The solid fraction as a
+    function of temperature \c alphaSolidT is specified as a Foam::Function1.
+
+    The model writes the field \c alpha[01].solid which can be visualised to to
+    show the solid distribution.
+
+Usage
+    Example usage:
+    \verbatim
+    VoFSolidificationMeltingSource1
+    {
+        type            VoFSolidificationMeltingSource;
+        active          yes;
+
+        selectionMode   cellZone;
+        cellZone        solidZone;
+
+        alphaSolidT     table
+        (
+            (330 1)
+            (335 0)
+        );
+
+        L               334000;
+    }
+    \endverbatim
+
+    Where:
+    \table
+        Property     | Description             | Required    | Default value
+        alphaSolidT  | Solid fraction as function of temperature | yes |
+        L            | Latent heat of fusion [J/kg] | yes    |
+        relax        | Relaxation coefficient [0-1] | no     | 0.9
+        Cu           | Model coefficient       | no          | 100000
+        q            | Model coefficient       | no          | 0.001
+    \endtable
+
+See also
+    Foam::fv::solidificationMeltingSource
+    Foam::Function1
+
+SourceFiles
+    VoFSolidificationMeltingSource.C
+    VoFSolidificationMeltingSourceIO.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef VoFSolidificationMeltingSource_H
+#define VoFSolidificationMeltingSource_H
+
+#include "fvMesh.H"
+#include "volFields.H"
+#include "cellSetOption.H"
+#include "Function1.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace fv
+{
+
+/*---------------------------------------------------------------------------*\
+                 Class VoFSolidificationMeltingSource Declaration
+\*---------------------------------------------------------------------------*/
+
+class VoFSolidificationMeltingSource
+:
+    public cellSetOption
+{
+    // Private data
+
+        //- Solid fraction as a function of temperature
+        autoPtr<Function1<scalar>> alphaSolidT_;
+
+        //- Latent heat of fusion [J/kg]
+        dimensionedScalar L_;
+
+        //- Phase fraction under-relaxation coefficient
+        scalar relax_;
+
+        //- Mushy region momentum sink coefficient [1/s]; default = 10^5
+        scalar Cu_;
+
+        //- Coefficient used in porosity calc - default = 0.001
+        scalar q_;
+
+        //- Solid phase fraction
+        volScalarField alphaSolid_;
+
+        //- Current time index (used for updating)
+        label curTimeIndex_;
+
+
+    // Private Member Functions
+
+        //- Return the name of the solid phase fraction
+        word alphaSolidName() const;
+
+        //- Update the model
+        void update();
+
+        //- Helper function to apply to the energy equation
+        template<class RhoFieldType>
+        void apply(const RhoFieldType& rho, fvMatrix<scalar>& eqn);
+
+        //- Disallow default bitwise copy construct
+        VoFSolidificationMeltingSource(const VoFSolidificationMeltingSource&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const VoFSolidificationMeltingSource&);
+
+
+public:
+
+    //- Runtime type information
+    TypeName("VoFSolidificationMeltingSource");
+
+
+    // Constructors
+
+        //- Construct from explicit source name and mesh
+        VoFSolidificationMeltingSource
+        (
+            const word& sourceName,
+            const word& modelType,
+            const dictionary& dict,
+            const fvMesh& mesh
+        );
+
+
+    // Member Functions
+
+        // Add explicit and implicit contributions
+
+            //- Add explicit contribution to enthalpy equation
+            virtual void addSup(fvMatrix<scalar>& eqn, const label fieldi);
+
+            //- Add implicit contribution to momentum equation
+            virtual void addSup(fvMatrix<vector>& eqn, const label fieldi);
+
+
+        // Add explicit and implicit contributions to compressible equation
+
+            //- Add explicit contribution to compressible enthalpy equation
+            virtual void addSup
+            (
+                const volScalarField& rho,
+                fvMatrix<scalar>& eqn,
+                const label fieldi
+            );
+
+            //- Add implicit contribution to compressible momentum equation
+            virtual void addSup
+            (
+                const volScalarField& rho,
+                fvMatrix<vector>& eqn,
+                const label fieldi
+            );
+
+
+        // IO
+
+            //- Read source dictionary
+            virtual bool read(const dictionary& dict);
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace fv
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+    #include "VoFSolidificationMeltingSourceTemplates.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

applications/solvers/multiphase/compressibleInterFoam/compressibleInterFilmFoam/VoFSolidificationMeltingSource/VoFSolidificationMeltingSourceIO.C

@@ -0,0 +1,51 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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 "VoFSolidificationMeltingSource.H"
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+bool Foam::fv::VoFSolidificationMeltingSource::read(const dictionary& dict)
+{
+    if (cellSetOption::read(dict))
+    {
+        alphaSolidT_ = Function1<scalar>::New("alphaSolidT", coeffs_);
+        coeffs_.lookup("L") >> L_;
+        coeffs_.readIfPresent("relax", relax_);
+        coeffs_.readIfPresent("Cu", Cu_);
+        coeffs_.readIfPresent("q", q_);
+
+        return true;
+    }
+    else
+    {
+        return false;
+    }
+
+    return false;
+}
+
+
+// ************************************************************************* //

applications/solvers/multiphase/compressibleInterFoam/compressibleInterFilmFoam/VoFSolidificationMeltingSource/VoFSolidificationMeltingSourceTemplates.C

@@ -0,0 +1,66 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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 "fvcDdt.H"
+#include "twoPhaseMixtureThermo.H"
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+template<class RhoFieldType>
+void Foam::fv::VoFSolidificationMeltingSource::apply
+(
+    const RhoFieldType& rho,
+    fvMatrix<scalar>& eqn
+)
+{
+    if (debug)
+    {
+        Info<< type() << ": applying source to " << eqn.psi().name() << endl;
+    }
+
+    update();
+
+    const twoPhaseMixtureThermo& thermo
+    (
+        mesh_.lookupObject<twoPhaseMixtureThermo>
+        (
+            twoPhaseMixtureThermo::dictName
+        )
+    );
+
+    const volScalarField CpVoF(thermo.thermo1().Cp());
+
+    if (eqn.psi().dimensions() == dimTemperature)
+    {
+        eqn += L_/CpVoF*(fvc::ddt(rho, alphaSolid_));
+    }
+    else
+    {
+        eqn += L_*(fvc::ddt(rho, alphaSolid_));
+    }
+}
+
+
+// ************************************************************************* //