From 1876f95e2a238c7af73dc8284f148f45bb401fe7 Mon Sep 17 00:00:00 2001
From: andy <andy>
Date: Wed, 24 Sep 2014 13:28:56 +0100
Subject: [PATCH 1/4] ENH: Surface film modelling - updated thixotropic
 viscosity model

---
 .../thixotropicViscosity.C                    | 42 +++++++++++--------
 .../thixotropicViscosity.H                    | 28 ++-----------
 2 files changed, 28 insertions(+), 42 deletions(-)

diff --git a/src/regionModels/surfaceFilmModels/submodels/thermo/filmViscosityModel/thixotropicViscosity/thixotropicViscosity.C b/src/regionModels/surfaceFilmModels/submodels/thermo/filmViscosityModel/thixotropicViscosity/thixotropicViscosity.C
index 5847ce12..3817aeeb 100644
--- a/src/regionModels/surfaceFilmModels/submodels/thermo/filmViscosityModel/thixotropicViscosity/thixotropicViscosity.C
+++ b/src/regionModels/surfaceFilmModels/submodels/thermo/filmViscosityModel/thixotropicViscosity/thixotropicViscosity.C
@@ -91,12 +91,12 @@ thixotropicViscosity::thixotropicViscosity
 )
 :
     filmViscosityModel(typeName, owner, dict, mu),
-    a_(coeffDict_.lookup("a")),
-    b_(coeffDict_.lookup("b")),
-    c_(coeffDict_.lookup("c")),
-    d_(coeffDict_.lookup("d")),
-    mu0_(coeffDict_.lookup("mu0")),
-    muInf_(coeffDict_.lookup("muInf")),
+    a_("a", dimless/dimTime, coeffDict_.lookup("a")),
+    b_("b", dimless, coeffDict_.lookup("b")),
+    d_("d", dimless, coeffDict_.lookup("d")),
+    c_("c", pow(dimTime, d_.value() - scalar(1)), coeffDict_.lookup("c")),
+    mu0_("mu0", dimPressure*dimTime, coeffDict_.lookup("mu0")),
+    muInf_("muInf", mu0_.dimensions(), coeffDict_.lookup("muInf")),
     K_(1.0 - Foam::sqrt(muInf_/mu0_)),
     lambda_
     (
@@ -144,27 +144,33 @@ void thixotropicViscosity::correct
     const volScalarField& delta = film.delta();
     const volScalarField& deltaRho = film.deltaRho();
     const surfaceScalarField& phi = film.phi();
+    const volScalarField& alpha = film.alpha();
 
-    // gamma-dot (shear rate) raised to the power d
-    volScalarField gDotPowD
-    (
-        "gDotPowD",
-        pow(mag(U - Uw)/(delta + film.deltaSmall()), d_)
-    );
+    // gamma-dot (shear rate)
+    volScalarField gDot("gDot", alpha*mag(U - Uw)/(delta + film.deltaSmall()));
 
-    dimensionedScalar c0("SMALL", dimMass/sqr(dimLength)/dimTime, SMALL);
-    volScalarField coeff(-deltaRho*c_*gDotPowD + c0);
+    if (debug && this->owner().regionMesh().time().outputTime())
+    {
+        gDot.write();
+    }
+
+    dimensionedScalar deltaRho0("deltaRho0", deltaRho.dimensions(), ROOTVSMALL);
+    surfaceScalarField phiU(phi/fvc::interpolate(deltaRho + deltaRho0));
+
+    dimensionedScalar c0("c0", dimless/dimTime, ROOTVSMALL);
+    volScalarField coeff("coeff", -c_*pow(gDot, d_) + c0);
 
     fvScalarMatrix lambdaEqn
     (
-        fvm::ddt(deltaRho, lambda_)
-      + fvm::div(phi, lambda_)
-      - fvm::Sp(fvc::div(phi), lambda_)
+        fvm::ddt(lambda_)
+      + fvm::div(phiU, lambda_)
+      - fvm::Sp(fvc::div(phiU), lambda_)
       ==
-        deltaRho*a_*pow((1.0 - lambda_), b_)
+        a_*pow((1.0 - lambda_), b_)
       + fvm::SuSp(coeff, lambda_)
     );
 
+
     lambdaEqn.relax();
 
     lambdaEqn.solve();
diff --git a/src/regionModels/surfaceFilmModels/submodels/thermo/filmViscosityModel/thixotropicViscosity/thixotropicViscosity.H b/src/regionModels/surfaceFilmModels/submodels/thermo/filmViscosityModel/thixotropicViscosity/thixotropicViscosity.H
index 2c60dc42..ac1a2ac2 100644
--- a/src/regionModels/surfaceFilmModels/submodels/thermo/filmViscosityModel/thixotropicViscosity/thixotropicViscosity.H
+++ b/src/regionModels/surfaceFilmModels/submodels/thermo/filmViscosityModel/thixotropicViscosity/thixotropicViscosity.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2013-2014 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -105,38 +105,18 @@ protected:
 
     // Protected data
 
-/*
-        //- Model `a' coefficient
-        scalar a_;
-
-        //- Model `b' coefficient
-        scalar b_;
-
-        //- Model `c' coefficient
-        scalar c_;
-
-        //- Model `d' coefficient
-        scalar d_;
-
-        //- Limiting viscosity when lambda = 1
-        scalar mu0_;
-
-        //- Limiting viscosity when lambda = 0
-        scalar muInf_;
-*/
-
         //- Model `a' coefficient
         dimensionedScalar a_;
 
         //- Model `b' coefficient
         dimensionedScalar b_;
 
-        //- Model `c' coefficient
-        dimensionedScalar c_;
-
         //- Model `d' coefficient
         dimensionedScalar d_;
 
+        //- Model `c' coefficient (read after d since dims depend on d value)
+        dimensionedScalar c_;
+
         //- Limiting viscosity when lambda = 1
         dimensionedScalar mu0_;
 

From ccc9593194b09fe0a975c6d2de717588d39d77f9 Mon Sep 17 00:00:00 2001
From: andy <andy>
Date: Thu, 16 Oct 2014 10:20:19 +0100
Subject: [PATCH 2/4] ENH: chtMultiRegionFoam solvers - added option to
 calculate using frozen flow

---
 .../fluid/createFluidFields.H                 | 35 ++++++-------------
 .../fluid/setRegionFluidFields.H              |  1 +
 .../fluid/solveFluid.H                        | 22 ++++++++----
 .../fluid/createFluidFields.H                 |  7 +++-
 .../fluid/setRegionFluidFields.H              |  2 ++
 .../chtMultiRegionFoam/fluid/solveFluid.H     | 33 ++++++++++-------
 6 files changed, 54 insertions(+), 46 deletions(-)

diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/createFluidFields.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/createFluidFields.H
index 99f14113..885b1f08 100644
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/createFluidFields.H
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/createFluidFields.H
@@ -11,8 +11,9 @@
     PtrList<radiation::radiationModel> radiation(fluidRegions.size());
 
     List<scalar> initialMassFluid(fluidRegions.size());
-    List<label> pRefCellFluid(fluidRegions.size(),0);
-    List<scalar> pRefValueFluid(fluidRegions.size(),0.0);
+    List<label> pRefCellFluid(fluidRegions.size(), 0);
+    List<scalar> pRefValueFluid(fluidRegions.size(), 0.0);
+    List<bool> frozenFlowFluid(fluidRegions.size(), false);
 
     PtrList<dimensionedScalar> rhoMax(fluidRegions.size());
     PtrList<dimensionedScalar> rhoMin(fluidRegions.size());
@@ -160,38 +161,22 @@
 
         initialMassFluid[i] = fvc::domainIntegrate(rhoFluid[i]).value();
 
+        const dictionary& simpleDict =
+            fluidRegions[i].solutionDict().subDict("SIMPLE");
+
         setRefCell
         (
             thermoFluid[i].p(),
             p_rghFluid[i],
-            fluidRegions[i].solutionDict().subDict("SIMPLE"),
+            simpleDict,
             pRefCellFluid[i],
             pRefValueFluid[i]
         );
 
-        rhoMax.set
-        (
-            i,
-            new dimensionedScalar
-            (
-                fluidRegions[i].solutionDict().subDict("SIMPLE").lookup
-                (
-                    "rhoMax"
-                )
-            )
-        );
+        simpleDict.readIfPresent("frozenFlow", frozenFlowFluid[i]);
 
-        rhoMin.set
-        (
-            i,
-            new dimensionedScalar
-            (
-                fluidRegions[i].solutionDict().subDict("SIMPLE").lookup
-                (
-                    "rhoMin"
-                )
-            )
-        );
+        rhoMax.set(i, new dimensionedScalar(simpleDict.lookup("rhoMax")));
+        rhoMin.set(i, new dimensionedScalar(simpleDict.lookup("rhoMin")));
 
         Info<< "    Adding fvOptions\n" << endl;
         fluidFvOptions.set
diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/setRegionFluidFields.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/setRegionFluidFields.H
index 34aea764..fe5ce1e8 100644
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/setRegionFluidFields.H
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/setRegionFluidFields.H
@@ -25,6 +25,7 @@
 
     const label pRefCell = pRefCellFluid[i];
     const scalar pRefValue = pRefValueFluid[i];
+    const bool frozenFlow = frozenFlowFluid[i];
 
     volScalarField& p_rgh = p_rghFluid[i];
     const volScalarField& gh = ghFluid[i];
diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/solveFluid.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/solveFluid.H
index a17893fb..b43a9490 100644
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/solveFluid.H
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/solveFluid.H
@@ -1,11 +1,19 @@
 //  Pressure-velocity SIMPLE corrector
 
-    p_rgh.storePrevIter();
-    rho.storePrevIter();
     {
-        #include "UEqn.H"
-        #include "EEqn.H"
-        #include "pEqn.H"
-    }
+        if (frozenFlow)
+        {
+            #include "EEqn.H"
+        }
+        else
+        {
+            p_rgh.storePrevIter();
+            rho.storePrevIter();
 
-    turb.correct();
+            #include "UEqn.H"
+            #include "EEqn.H"
+            #include "pEqn.H"
+
+            turb.correct();
+        }
+    }
diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/createFluidFields.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/createFluidFields.H
index 9afa705c..98ef8d7f 100644
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/createFluidFields.H
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/createFluidFields.H
@@ -13,6 +13,7 @@
     PtrList<volScalarField> dpdtFluid(fluidRegions.size());
 
     List<scalar> initialMassFluid(fluidRegions.size());
+    List<bool> frozenFlowFluid(fluidRegions.size(), false);
 
     PtrList<fv::IOoptionList> fluidFvOptions(fluidRegions.size());
 
@@ -102,7 +103,7 @@
             )
         );
 
-        Info<< "    Adding turbulence \n" << endl;
+        Info<< "    Adding to turbulence\n" << endl;
         turbulence.set
         (
             i,
@@ -190,6 +191,10 @@
             )
         );
 
+        const dictionary& pimpleDict =
+            fluidRegions[i].solutionDict().subDict("PIMPLE");
+        pimpleDict.readIfPresent("frozenFlow", frozenFlowFluid[i]);
+
         Info<< "    Adding fvOptions\n" << endl;
         fluidFvOptions.set
         (
diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/setRegionFluidFields.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/setRegionFluidFields.H
index 867585f4..ef3a3df7 100644
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/setRegionFluidFields.H
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/setRegionFluidFields.H
@@ -28,3 +28,5 @@
         dimMass,
         initialMassFluid[i]
     );
+
+    const bool frozenFlow = frozenFlowFluid[i];
\ No newline at end of file
diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/solveFluid.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/solveFluid.H
index 354c4a13..6f511e02 100644
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/solveFluid.H
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/solveFluid.H
@@ -3,24 +3,31 @@ if (finalIter)
     mesh.data::add("finalIteration", true);
 }
 
-if (oCorr == 0)
+if (frozenFlow)
 {
-    #include "rhoEqn.H"
+    #include "EEqn.H"
 }
-
-#include "UEqn.H"
-#include "EEqn.H"
-
-// --- PISO loop
-for (int corr=0; corr<nCorr; corr++)
+else
 {
-    #include "pEqn.H"
+    if (oCorr == 0)
+    {
+        #include "rhoEqn.H"
+    }
+
+    #include "UEqn.H"
+    #include "EEqn.H"
+
+    // --- PISO loop
+    for (int corr=0; corr<nCorr; corr++)
+    {
+        #include "pEqn.H"
+    }
+
+    turb.correct();
+
+    rho = thermo.rho();
 }
 
-turb.correct();
-
-rho = thermo.rho();
-
 if (finalIter)
 {
     mesh.data::remove("finalIteration");

From eeada10f00136703a54b4ec3bd7ab4ac0c4d4822 Mon Sep 17 00:00:00 2001
From: andy <andy>
Date: Thu, 16 Oct 2014 10:21:01 +0100
Subject: [PATCH 3/4] ENH: Added new solidification and melting fvOption

---
 src/fvOptions/Make/files                      |   4 +-
 .../solidificationMeltingSource.C             | 358 ++++++++++++++++++
 .../solidificationMeltingSource.H             | 280 ++++++++++++++
 .../solidificationMeltingSourceIO.C           |  68 ++++
 .../solidificationMeltingSourceTemplates.C    |  70 ++++
 5 files changed, 779 insertions(+), 1 deletion(-)
 create mode 100644 src/fvOptions/sources/derived/solidificationMeltingSource/solidificationMeltingSource.C
 create mode 100644 src/fvOptions/sources/derived/solidificationMeltingSource/solidificationMeltingSource.H
 create mode 100644 src/fvOptions/sources/derived/solidificationMeltingSource/solidificationMeltingSourceIO.C
 create mode 100644 src/fvOptions/sources/derived/solidificationMeltingSource/solidificationMeltingSourceTemplates.C

diff --git a/src/fvOptions/Make/files b/src/fvOptions/Make/files
index d45448f2..b85aac6d 100644
--- a/src/fvOptions/Make/files
+++ b/src/fvOptions/Make/files
@@ -12,6 +12,7 @@ $(generalSources)/semiImplicitSource/semiImplicitSource.C
 
 derivedSources=sources/derived
 $(derivedSources)/actuationDiskSource/actuationDiskSource.C
+$(derivedSources)/effectivenessHeatExchangerSource/effectivenessHeatExchangerSource.C
 $(derivedSources)/explicitPorositySource/explicitPorositySource.C
 $(derivedSources)/MRFSource/MRFSource.C
 $(derivedSources)/pressureGradientExplicitSource/pressureGradientExplicitSource.C
@@ -27,7 +28,8 @@ $(derivedSources)/rotorDiskSource/trimModel/trimModel/trimModel.C
 $(derivedSources)/rotorDiskSource/trimModel/trimModel/trimModelNew.C
 $(derivedSources)/rotorDiskSource/trimModel/fixed/fixedTrim.C
 $(derivedSources)/rotorDiskSource/trimModel/targetCoeff/targetCoeffTrim.C
-$(derivedSources)/effectivenessHeatExchangerSource/effectivenessHeatExchangerSource.C
+$(derivedSources)/solidificationMeltingSource/solidificationMeltingSource.C
+$(derivedSources)/solidificationMeltingSource/solidificationMeltingSourceIO.C
 
 interRegion = sources/interRegion
 $(interRegion)/interRegionHeatTransferModel/constantHeatTransfer/constantHeatTransfer.C
diff --git a/src/fvOptions/sources/derived/solidificationMeltingSource/solidificationMeltingSource.C b/src/fvOptions/sources/derived/solidificationMeltingSource/solidificationMeltingSource.C
new file mode 100644
index 00000000..e17b8c88
--- /dev/null
+++ b/src/fvOptions/sources/derived/solidificationMeltingSource/solidificationMeltingSource.C
@@ -0,0 +1,358 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013-2014 OpenCFD Ltd.
+     \\/     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 "solidificationMeltingSource.H"
+#include "fvMatrices.H"
+#include "basicThermo.H"
+#include "uniformDimensionedFields.H"
+#include "zeroGradientFvPatchFields.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::thermoMode, 2>
+    Foam::fv::solidificationMeltingSource::thermoModeTypeNames_;
+
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+bool Foam::fv::solidificationMeltingSource::solveField
+(
+    const word& fieldName
+) const
+{
+    bool result = true;
+
+    switch (mode_)
+    {
+        case mdThermo:
+        {
+            const basicThermo& thermo =
+                mesh_.lookupObject<basicThermo>("thermophysicalProperties");
+
+            if (fieldName != thermo.he().name())
+            {
+                result = false;
+            }
+            break;
+        }
+        case mdLookup:
+        {
+            if (fieldName != TName_)
+            {
+                result = false;
+            }
+            break;
+        }
+        default:
+        {
+            FatalErrorIn
+            (
+                "bool Foam::fv::solidificationMeltingSource::solveField"
+                "("
+                    "const word&"
+                ") const"
+            )
+                << "Unhandled thermo mode: " << thermoModeTypeNames_[mode_]
+                << abort(FatalError);
+        }
+    }
+    return result;
+}
+
+
+Foam::tmp<Foam::volScalarField>
+Foam::fv::solidificationMeltingSource::Cp() const
+{
+    switch (mode_)
+    {
+        case mdThermo:
+        {
+            const basicThermo& thermo =
+                mesh_.lookupObject<basicThermo>("thermophysicalProperties");
+
+            return thermo.Cp();
+            break;
+        }
+        case mdLookup:
+        {
+            if (CpName_ == "CpRef")
+            {
+                scalar CpRef = readScalar(coeffs_.lookup("CpRef"));
+
+                return tmp<volScalarField>
+                (
+                    new volScalarField
+                    (
+                        IOobject
+                        (
+                            name_ + ":Cp",
+                            mesh_.time().timeName(),
+                            mesh_,
+                            IOobject::NO_READ,
+                            IOobject::NO_WRITE
+                        ),
+                        mesh_,
+                        dimensionedScalar
+                        (
+                            "Cp",
+                            dimEnergy/dimMass/dimTemperature,
+                            CpRef
+                        ),
+                        zeroGradientFvPatchScalarField::typeName
+                    )
+                );
+            }
+            else
+            {
+                return mesh_.lookupObject<volScalarField>(CpName_);
+            }
+
+            break;
+        }
+        default:
+        {
+            FatalErrorIn
+            (
+                "Foam::tmp<Foam::volScalarField> "
+                "Foam::fv::solidificationMeltingSource::Cp() const"
+            )
+                << "Unhandled thermo mode: " << thermoModeTypeNames_[mode_]
+                << abort(FatalError);
+        }
+    }
+
+    return tmp<volScalarField>(NULL);
+}
+
+
+Foam::vector Foam::fv::solidificationMeltingSource::g() const
+{
+    if (mesh_.foundObject<uniformDimensionedVectorField>("g"))
+    {
+        const uniformDimensionedVectorField& value =
+            mesh_.lookupObject<uniformDimensionedVectorField>("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<volScalarField>(TName_);
+
+    forAll(cells_, i)
+    {
+        label cellI = cells_[i];
+
+        scalar Tc = T[cellI];
+        scalar Cpc = Cp[cellI];
+        scalar alpha1New = alpha1_[cellI] + relax_*Cpc*(Tc - Tmelt_)/L_;
+
+        alpha1_[cellI] = max(0, min(alpha1New, 1));
+        deltaT_[i] = Tc - Tmelt_;
+    }
+
+    alpha1_.correctBoundaryConditions();
+
+    curTimeIndex_ = mesh_.time().timeIndex();
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::fv::solidificationMeltingSource::solidificationMeltingSource
+(
+    const word& sourceName,
+    const word& modelType,
+    const dictionary& dict,
+    const fvMesh& mesh
+)
+:
+    option(sourceName, modelType, dict, mesh),
+    Tmelt_(readScalar(coeffs_.lookup("Tmelt"))),
+    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<word>("TName", "T")),
+    CpName_(coeffs_.lookupOrDefault<word>("CpName", "Cp")),
+    UName_(coeffs_.lookupOrDefault<word>("UName", "U")),
+    phiName_(coeffs_.lookupOrDefault<word>("phiName", "phi")),
+    Cu_(coeffs_.lookupOrDefault<scalar>("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.0),
+        zeroGradientFvPatchScalarField::typeName
+    ),
+    curTimeIndex_(-1),
+    deltaT_(cells_.size(), 0)
+{
+    fieldNames_.setSize(1, "source");
+    applied_.setSize(1, false);
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+bool Foam::fv::solidificationMeltingSource::alwaysApply() const
+{
+    return true;
+}
+
+
+void Foam::fv::solidificationMeltingSource::addSup
+(
+    fvMatrix<scalar>& eqn,
+    const label fieldI
+)
+{
+    apply(geometricOneField(), eqn);
+}
+
+
+void Foam::fv::solidificationMeltingSource::addSup
+(
+    const volScalarField& rho,
+    fvMatrix<scalar>& eqn,
+    const label fieldI
+)
+{
+    apply(rho, eqn);
+}
+
+
+void Foam::fv::solidificationMeltingSource::addSup
+(
+    fvMatrix<vector>& eqn,
+    const label fieldI
+)
+{
+    const volVectorField& U = eqn.psi();
+
+    if (U.name() != UName_)
+    {
+        return;
+    }
+
+    if (debug)
+    {
+        Info<< type() << ": applying source to " << UName_ << 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)
+    {
+        label cellI = cells_[i];
+
+        scalar Vc = V[cellI];
+        scalar alpha1c = alpha1_[cellI];
+
+        scalar S = -Cu_*sqr(1.0 - alpha1c)/(pow3(alpha1c) + q_);
+        vector Sb = rhoRef_*g*beta_*deltaT_[i];
+
+        Sp[cellI] += Vc*S;
+        Su[cellI] += Vc*Sb;
+    }
+}
+
+
+void Foam::fv::solidificationMeltingSource::addSup
+(
+    const volScalarField& rho,
+    fvMatrix<vector>& 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 00000000..a8d6f951
--- /dev/null
+++ b/src/fvOptions/sources/derived/solidificationMeltingSource/solidificationMeltingSource.H
@@ -0,0 +1,280 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013-2014 OpenCFD Ltd.
+     \\/     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::solidificationMeltingSource
+
+Description
+    This source is designed to model the effect of solidification and melting
+    processes, e.g. windhield defrosting.  The phase change occurs at the
+    melting temperature, \c Tmelt.
+
+    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.
+
+    Based on the references:
+
+      1. V.R. Voller and C. Prakash, A fixed grid numerical modelling
+         methodology for convection-diffusion mushy phase-change problems,
+         Int. J. Heat Mass Transfer 30(8):17091719, 1987.
+      2. C.R. Swaminathan. and V.R. Voller, A general enthalpy model for
+         modeling solidification processes, Metallurgical Transactions
+         23B:651664, 1992.
+
+    The model generates a field \c \<name\>:alpha1 which can be visualised to
+    to show the melt distribution as a fraction [0-1]
+
+    \heading Source usage
+    Example usage:
+    \verbatim
+    solidificationMeltingSource1
+    {
+        type            solidificationMeltingSource;
+        active          on;
+        selectionMode   cellZone;
+        cellZone        iceZone;
+
+        solidificationMeltingSourceCoeffs
+        {
+            Tmelt           273;
+            L               334000;
+            thermoMode      thermo;
+            beta            50e-6;
+            rhoRef          800;
+        }
+    }
+    \endverbatim
+
+    Where:
+    \table
+        Property     | Description             | Required    | Default value
+        Tmelt        | Melting temperature [K] | yes         |
+        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 | yes       |
+        rhoName      | Name of density field   | no          | rho
+        TName        | Name of temperature field | no        | T
+        CpName       | Name of specific heat capacity field | no | Cp
+        UName        | Name of velocity field  | no          | U
+        phiName      | Name of flux field      | no          | phi
+        Cu           | Model coefficient       | 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 "fvOption.H"
+#include "NamedEnum.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace fv
+{
+
+/*---------------------------------------------------------------------------*\
+                 Class solidificationMeltingSource Declaration
+\*---------------------------------------------------------------------------*/
+
+class solidificationMeltingSource
+:
+    public option
+{
+public:
+
+    enum thermoMode
+    {
+        mdThermo,
+        mdLookup
+    };
+
+    static const NamedEnum<thermoMode, 2> thermoModeTypeNames_;
+
+
+private:
+
+    // Private data
+
+        //- Temperature at which melting occurs [K]
+        scalar Tmelt_;
+
+        //- 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
+
+        //- Flag to indicate whether to solve for given field
+        bool solveField(const word& fieldName) const;
+
+        //- Return the specific heat capacity field
+        tmp<volScalarField> 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<class RhoFieldType>
+        void apply(const RhoFieldType& rho, fvMatrix<scalar>& 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
+
+        //- Flag to bypass the apply flag list checking
+        virtual bool alwaysApply() const;
+
+
+        // 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
+            );
+
+
+        // I-O
+
+            //- Write the source properties
+            virtual void writeData(Ostream&) const;
+
+            //- Read source dictionary
+            virtual bool read(const dictionary& dict);
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace fv
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+    #include "solidificationMeltingSourceTemplates.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/fvOptions/sources/derived/solidificationMeltingSource/solidificationMeltingSourceIO.C b/src/fvOptions/sources/derived/solidificationMeltingSource/solidificationMeltingSourceIO.C
new file mode 100644
index 00000000..d6173d12
--- /dev/null
+++ b/src/fvOptions/sources/derived/solidificationMeltingSource/solidificationMeltingSourceIO.C
@@ -0,0 +1,68 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013-2014 OpenCFD Ltd.
+     \\/     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 "solidificationMeltingSource.H"
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+void Foam::fv::solidificationMeltingSource::writeData(Ostream& os) const
+{
+    os  << indent << name_ << endl;
+    dict_.write(os);
+}
+
+
+bool Foam::fv::solidificationMeltingSource::read(const dictionary& dict)
+{
+    if (option::read(dict))
+    {
+        coeffs_.lookup("Tmelt") >> Tmelt_;
+        coeffs_.lookup("L") >> L_;
+
+        coeffs_.readIfPresent("relax", relax_);
+
+        mode_ = thermoModeTypeNames_.read(coeffs_.lookup("thermoMode"));
+
+        coeffs_.lookup("rhoRef") >> rhoRef_;
+        coeffs_.readIfPresent("TName", TName_);
+        coeffs_.readIfPresent("UName", UName_);
+
+        coeffs_.readIfPresent("Cu", Cu_);
+        coeffs_.readIfPresent("q", q_);
+
+        coeffs_.readIfPresent("beta", beta_);
+
+        return true;
+    }
+    else
+    {
+        return false;
+    }
+
+    return false;
+}
+
+
+// ************************************************************************* //
diff --git a/src/fvOptions/sources/derived/solidificationMeltingSource/solidificationMeltingSourceTemplates.C b/src/fvOptions/sources/derived/solidificationMeltingSource/solidificationMeltingSourceTemplates.C
new file mode 100644
index 00000000..2dfddb76
--- /dev/null
+++ b/src/fvOptions/sources/derived/solidificationMeltingSource/solidificationMeltingSourceTemplates.C
@@ -0,0 +1,70 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2014 OpenCFD Ltd.
+     \\/     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 "fvMatrices.H"
+#include "fvcDdt.H"
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+template<class RhoFieldType>
+void Foam::fv::solidificationMeltingSource::apply
+(
+    const RhoFieldType& rho,
+    fvMatrix<scalar>& eqn
+)
+{
+    if (!solveField(eqn.psi().name()))
+    {
+        return;
+    }
+
+    if (debug)
+    {
+        Info<< type() << ": applying source to " << eqn.psi().name() << endl;
+    }
+
+    const volScalarField Cp(this->Cp());
+
+    update(Cp);
+
+    dimensionedScalar L("L", dimEnergy/dimMass, L_);
+
+    // contributions added to rhs of solver equation
+    if (eqn.psi().dimensions() == dimTemperature)
+    {
+        // isothermal phase change - only include time derivative
+//        eqn -= L/Cp*(fvc::ddt(rho, alpha1_) + fvc::div(phi, alpha1_));
+        eqn -= L/Cp*(fvc::ddt(rho, alpha1_));
+    }
+    else
+    {
+        // isothermal phase change - only include time derivative
+//        eqn -= L*(fvc::ddt(rho, alpha1_) + fvc::div(phi, alpha1_));
+        eqn -= L*(fvc::ddt(rho, alpha1_));
+    }
+}
+
+
+// ************************************************************************* //

From 4eaf1fe74669a7b51211485d3727023cba0b661c Mon Sep 17 00:00:00 2001
From: william <william>
Date: Thu, 16 Oct 2014 16:02:50 +0100
Subject: [PATCH 4/4] BUG: removed unnecessary sed processing from cleanCase

---
 bin/tools/CleanFunctions | 8 +-------
 1 file changed, 1 insertion(+), 7 deletions(-)

diff --git a/bin/tools/CleanFunctions b/bin/tools/CleanFunctions
index bd1659e6..b5edc197 100644
--- a/bin/tools/CleanFunctions
+++ b/bin/tools/CleanFunctions
@@ -2,7 +2,7 @@
 # =========                 |
 # \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
 #  \\    /   O peration     |
-#   \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+#   \\  /    A nd           | Copyright (C) 2011-2014 OpenFOAM Foundation
 #    \\/     M anipulation  |
 #------------------------------------------------------------------------------
 # License
@@ -108,12 +108,6 @@ cleanCase()
     then
         rm -f constant/polyMesh/blockMeshDict > /dev/null 2>&1
     fi
-
-    for f in `find . -name "*Dict"`
-    do
-        sed -e /arguments/d $f > temp.$$
-        mv temp.$$ $f
-    done
 }
 
 removeCase()