diff --git a/src/fvOptions/sources/derived/solidificationMeltingSource/solidificationMeltingSource.C b/src/fvOptions/sources/derived/solidificationMeltingSource/solidificationMeltingSource.C
new file mode 100644
index 000000000..5f336ca33
--- /dev/null
+++ b/src/fvOptions/sources/derived/solidificationMeltingSource/solidificationMeltingSource.C
@@ -0,0 +1,342 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
+ \\ / O peration |
+ \\ / A nd | Copyright (C) 2014-2018 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 .
+
+\*---------------------------------------------------------------------------*/
+
+#include "solidificationMeltingSource.H"
+#include "fvMatrices.H"
+#include "basicThermo.H"
+#include "uniformDimensionedFields.H"
+#include "zeroGradientFvPatchFields.H"
+#include "extrapolatedCalculatedFvPatchFields.H"
+#include "addToRunTimeSelectionTable.H"
+#include "geometricOneField.H"
+
+// * * * * * * * * * * * * * Static Member Functions * * * * * * * * * * * * //
+
+namespace Foam
+{
+ template<>
+ const char* NamedEnum
+ <
+ fv::solidificationMeltingSource::thermoMode,
+ 2
+ >::names[] =
+ {
+ "thermo",
+ "lookup"
+ };
+
+ namespace fv
+ {
+ defineTypeNameAndDebug(solidificationMeltingSource, 0);
+
+ addToRunTimeSelectionTable
+ (
+ option,
+ solidificationMeltingSource,
+ dictionary
+ );
+ }
+}
+
+const Foam::NamedEnum
+ Foam::fv::solidificationMeltingSource::thermoModeTypeNames_;
+
+
+// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
+
+Foam::tmp
+Foam::fv::solidificationMeltingSource::Cp() const
+{
+ switch (mode_)
+ {
+ case mdThermo:
+ {
+ const basicThermo& thermo =
+ mesh_.lookupObject(basicThermo::dictName);
+
+ return thermo.Cp();
+ break;
+ }
+ case mdLookup:
+ {
+ if (CpName_ == "CpRef")
+ {
+ scalar CpRef = readScalar(coeffs_.lookup("CpRef"));
+
+ return tmp
+ (
+ new volScalarField
+ (
+ IOobject
+ (
+ name_ + ":Cp",
+ mesh_.time().timeName(),
+ mesh_,
+ IOobject::NO_READ,
+ IOobject::NO_WRITE
+ ),
+ mesh_,
+ dimensionedScalar
+ (
+ "Cp",
+ dimEnergy/dimMass/dimTemperature,
+ CpRef
+ ),
+ extrapolatedCalculatedFvPatchScalarField::typeName
+ )
+ );
+ }
+ else
+ {
+ return mesh_.lookupObject(CpName_);
+ }
+
+ break;
+ }
+ default:
+ {
+ FatalErrorInFunction
+ << "Unhandled thermo mode: " << thermoModeTypeNames_[mode_]
+ << abort(FatalError);
+ }
+ }
+
+ return tmp(nullptr);
+}
+
+
+Foam::vector Foam::fv::solidificationMeltingSource::g() const
+{
+ if (mesh_.foundObject("g"))
+ {
+ const uniformDimensionedVectorField& value =
+ mesh_.lookupObject("g");
+ return value.value();
+ }
+ else
+ {
+ return coeffs_.lookup("g");
+ }
+}
+
+
+void Foam::fv::solidificationMeltingSource::update(const volScalarField& Cp)
+{
+ if (curTimeIndex_ == mesh_.time().timeIndex())
+ {
+ return;
+ }
+
+ if (debug)
+ {
+ Info<< type() << ": " << name_ << " - updating phase indicator" << endl;
+ }
+
+ // update old time alpha1 field
+ alpha1_.oldTime();
+
+ const volScalarField& T = mesh_.lookupObject(TName_);
+
+ forAll(cells_, i)
+ {
+ const label celli = cells_[i];
+
+ const scalar Tc = T[celli];
+ const scalar Cpc = Cp[celli];
+ const scalar alpha1New =
+ alpha1_[celli]
+ + relax_*Cpc
+ *(
+ Tc
+ - max
+ (
+ Tsol_,
+ Tsol_
+ + (Tliq_ - Tsol_)*(alpha1_[celli] - alpha1e_)/(1 - alpha1e_)
+ )
+ )/L_;
+
+ alpha1_[celli] = max(0, min(alpha1New, 1));
+ deltaT_[i] =
+ Tc
+ - max
+ (
+ Tsol_,
+ Tsol_
+ + (Tliq_ - Tsol_)*(alpha1_[celli] - alpha1e_)/(1 - alpha1e_)
+ );
+ }
+
+ alpha1_.correctBoundaryConditions();
+
+ curTimeIndex_ = mesh_.time().timeIndex();
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
+
+Foam::fv::solidificationMeltingSource::solidificationMeltingSource
+(
+ const word& sourceName,
+ const word& modelType,
+ const dictionary& dict,
+ const fvMesh& mesh
+)
+:
+ cellSetOption(sourceName, modelType, dict, mesh),
+ Tsol_(readScalar(coeffs_.lookup("Tsol"))),
+ Tliq_(coeffs_.lookupOrDefault("Tliq", Tsol_)),
+ alpha1e_(coeffs_.lookupOrDefault("alpha1e", 0)),
+ L_(readScalar(coeffs_.lookup("L"))),
+ relax_(coeffs_.lookupOrDefault("relax", 0.9)),
+ mode_(thermoModeTypeNames_.read(coeffs_.lookup("thermoMode"))),
+ rhoRef_(readScalar(coeffs_.lookup("rhoRef"))),
+ TName_(coeffs_.lookupOrDefault("T", "T")),
+ CpName_(coeffs_.lookupOrDefault("Cp", "Cp")),
+ UName_(coeffs_.lookupOrDefault("U", "U")),
+ phiName_(coeffs_.lookupOrDefault("phi", "phi")),
+ Cu_(coeffs_.lookupOrDefault("Cu", 100000)),
+ q_(coeffs_.lookupOrDefault("q", 0.001)),
+ beta_(readScalar(coeffs_.lookup("beta"))),
+ alpha1_
+ (
+ IOobject
+ (
+ name_ + ":alpha1",
+ mesh.time().timeName(),
+ mesh,
+ IOobject::READ_IF_PRESENT,
+ IOobject::AUTO_WRITE
+ ),
+ mesh,
+ dimensionedScalar("alpha1", dimless, 0),
+ zeroGradientFvPatchScalarField::typeName
+ ),
+ curTimeIndex_(-1),
+ deltaT_(cells_.size(), 0)
+{
+ fieldNames_.setSize(2);
+ fieldNames_[0] = UName_;
+
+ switch (mode_)
+ {
+ case mdThermo:
+ {
+ const basicThermo& thermo =
+ mesh_.lookupObject(basicThermo::dictName);
+
+ fieldNames_[1] = thermo.he().name();
+ break;
+ }
+ case mdLookup:
+ {
+ fieldNames_[1] = TName_;
+ break;
+ }
+ default:
+ {
+ FatalErrorInFunction
+ << "Unhandled thermo mode: " << thermoModeTypeNames_[mode_]
+ << abort(FatalError);
+ }
+ }
+
+ applied_.setSize(fieldNames_.size(), false);
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
+
+void Foam::fv::solidificationMeltingSource::addSup
+(
+ fvMatrix& eqn,
+ const label fieldi
+)
+{
+ apply(geometricOneField(), eqn);
+}
+
+
+void Foam::fv::solidificationMeltingSource::addSup
+(
+ const volScalarField& rho,
+ fvMatrix& eqn,
+ const label fieldi
+)
+{
+ apply(rho, eqn);
+}
+
+
+void Foam::fv::solidificationMeltingSource::addSup
+(
+ fvMatrix& eqn,
+ const label fieldi
+)
+{
+ if (debug)
+ {
+ Info<< type() << ": applying source to " << eqn.psi().name() << endl;
+ }
+
+ const volScalarField Cp(this->Cp());
+
+ update(Cp);
+
+ vector g = this->g();
+
+ scalarField& Sp = eqn.diag();
+ vectorField& Su = eqn.source();
+ const scalarField& V = mesh_.V();
+
+ forAll(cells_, i)
+ {
+ const label celli = cells_[i];
+
+ const scalar Vc = V[celli];
+ const scalar alpha1c = alpha1_[celli];
+
+ const scalar S = -Cu_*sqr(1.0 - alpha1c)/(pow3(alpha1c) + q_);
+ const vector Sb = rhoRef_*g*beta_*deltaT_[i];
+
+ Sp[celli] += Vc*S;
+ Su[celli] += Vc*Sb;
+ }
+}
+
+
+void Foam::fv::solidificationMeltingSource::addSup
+(
+ const volScalarField& rho,
+ fvMatrix& eqn,
+ const label fieldi
+)
+{
+ // Momentum source uses a Boussinesq approximation - redirect
+ addSup(eqn, fieldi);
+}
+
+
+// ************************************************************************* //
diff --git a/src/fvOptions/sources/derived/solidificationMeltingSource/solidificationMeltingSource.H b/src/fvOptions/sources/derived/solidificationMeltingSource/solidificationMeltingSource.H
new file mode 100644
index 000000000..441ee1a08
--- /dev/null
+++ b/src/fvOptions/sources/derived/solidificationMeltingSource/solidificationMeltingSource.H
@@ -0,0 +1,289 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
+ \\ / O peration |
+ \\ / A nd | Copyright (C) 2014-2018 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 .
+
+Class
+ Foam::fv::solidificationMeltingSource
+
+Description
+ This source is designed to model the effect of solidification and melting
+ processes, e.g. windhield defrosting.
+
+ The isotherm phase change occurs at the melting temperature, \c Tsol (= \c
+ Tliq). The not isotherm phase change occurs between solidus and liquidus
+ temperature, \c Tsol < \c Tliq respectively, as long as the melt fraction is
+ greater than the max eutectic melt fraction, \c alpha1e (0 =
+ pure_substance, 1 = eutectic_mixture is not permitted) , i.e. eutectic to
+ initial solvent concentration difference, where a linear eutectic melt
+ fraction to temperature relation is considered - lever rule.
+
+ 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.
+
+ References:
+ \verbatim
+ Voller, V. R., & Prakash, C. (1987).
+ A fixed grid numerical modelling methodology for convection-diffusion
+ mushy region phase-change problems.
+ International Journal of Heat and Mass Transfer, 30(8), 1709-1719.
+
+ Swaminathan, C. R., & Voller, V. R. (1992).
+ A general enthalpy method for modeling solidification processes.
+ Metallurgical transactions B, 23(5), 651-664.
+ \endverbatim
+
+ The model generates the field \c \:alpha1 which can be visualised to
+ to show the melt distribution as a fraction [0-1].
+
+Usage
+ Example usage:
+ \verbatim
+ solidificationMeltingSource1
+ {
+ type solidificationMeltingSource;
+ active yes;
+
+ selectionMode cellZone;
+ cellZone iceZone;
+
+ Tsol 273;
+ L 334000;
+ thermoMode thermo;
+ beta 50e-6;
+ rhoRef 800;
+ }
+ \endverbatim
+
+ Where:
+ \table
+ Property | Description | Required | Default value
+ Tsol | Solidus temperature [K] | yes |
+ Tliq | Liquidus temperature [K] | no | Tsol
+ alpha1e | Max eutectic melt fraction [0-1[ | no | 0
+ L | Latent heat of fusion [J/kg] | yes |
+ relax | Relaxation coefficient [0-1] | no | 0.9
+ thermoMode | Thermo mode [thermo|lookup] | yes |
+ rhoRef | Reference (solid) density [kg/m^3] | yes |
+ rho | Name of density field | no | rho
+ T | Name of temperature field | no | T
+ Cp | Name of specific heat field | no | Cp
+ U | Name of velocity field | no | U
+ phi | Name of flux field | no | phi
+ Cu | Model coefficient [1/s] | no | 100000
+ q | Model coefficient | no | 0.001
+ beta | Thermal expansion coefficient [1/K] | yes |
+ g | Accelerartion due to gravity | no |
+ \endtable
+
+SourceFiles
+ solidificationMeltingSource.C
+ solidificationMeltingSourceIO.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef solidificationMeltingSource_H
+#define solidificationMeltingSource_H
+
+#include "fvMesh.H"
+#include "volFields.H"
+#include "cellSetOption.H"
+#include "NamedEnum.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace fv
+{
+
+/*---------------------------------------------------------------------------*\
+ Class solidificationMeltingSource Declaration
+\*---------------------------------------------------------------------------*/
+
+class solidificationMeltingSource
+:
+ public cellSetOption
+{
+public:
+
+ enum thermoMode
+ {
+ mdThermo,
+ mdLookup
+ };
+
+ static const NamedEnum thermoModeTypeNames_;
+
+
+private:
+
+ // Private data
+
+ //- Temperature at which isotherm melting occurs
+ // and not isotherm melting starts e.g. eutectic [K]
+ scalar Tsol_;
+
+ //- Temperature at which not isotherm melting stops (Tliq > Tsol) [K]
+ scalar Tliq_;
+
+ //- Max eutectic melt fraction [0-1[
+ // where 0 = pure substance is default value,
+ // 1 = eutectic mixture is not permitted
+ scalar alpha1e_;
+
+ //- Latent heat of fusion [J/kg]
+ scalar L_;
+
+ //- Phase fraction under-relaxation coefficient
+ scalar relax_;
+
+ //- Thermodynamics mode
+ thermoMode mode_;
+
+ //- Reference density - typically the solid density
+ scalar rhoRef_;
+
+ //- Name of temperature field - default = "T" (optional)
+ word TName_;
+
+ //- Name of specific heat capacity field - default = "Cp" (optional)
+ word CpName_;
+
+ //- Name of velocity field - default = "U" (optional)
+ word UName_;
+
+ //- Name of flux field - default = "phi" (optional)
+ word phiName_;
+
+ //- Mushy region momentum sink coefficient [1/s]; default = 10^5
+ scalar Cu_;
+
+ //- Coefficient used in porosity calc - default = 0.001
+ scalar q_;
+
+ //- Thermal expansion coefficient [1/K]
+ scalar beta_;
+
+ //- Phase fraction indicator field
+ volScalarField alpha1_;
+
+ //- Current time index (used for updating)
+ label curTimeIndex_;
+
+ //- Temperature change cached for source calculation when alpha1 updated
+ scalarField deltaT_;
+
+
+ // Private Member Functions
+
+ //- Return the specific heat capacity field
+ tmp Cp() const;
+
+ //- Return the gravity vector
+ vector g() const;
+
+ //- Update the model
+ void update(const volScalarField& Cp);
+
+ //- Helper function to apply to the energy equation
+ template
+ void apply(const RhoFieldType& rho, fvMatrix& eqn);
+
+ //- Disallow default bitwise copy construct
+ solidificationMeltingSource(const solidificationMeltingSource&);
+
+ //- Disallow default bitwise assignment
+ void operator=(const solidificationMeltingSource&);
+
+
+public:
+
+ //- Runtime type information
+ TypeName("solidificationMeltingSource");
+
+
+ // Constructors
+
+ //- Construct from explicit source name and mesh
+ solidificationMeltingSource
+ (
+ 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& eqn, const label fieldi);
+
+ //- Add implicit contribution to momentum equation
+ virtual void addSup(fvMatrix& 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& eqn,
+ const label fieldi
+ );
+
+ //- Add implicit contribution to compressible momentum equation
+ virtual void addSup
+ (
+ const volScalarField& rho,
+ fvMatrix& eqn,
+ const label fieldi
+ );
+
+
+ // IO
+
+ //- Read source dictionary
+ virtual bool read(const dictionary& dict);
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace fv
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+ #include "solidificationMeltingSourceTemplates.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //