diff --git a/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMesh.C b/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMesh.C
index 7ea9e43be5..e1d31a499d 100644
--- a/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMesh.C
+++ b/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMesh.C
@@ -1510,7 +1510,7 @@ int main(int argc, char *argv[])
(
snapDict,
motionDict,
- mergePatchFaces
+ mergePatchFaces,
curvature,
planarAngle,
snapParams
diff --git a/src/dynamicFvMesh/dynamicRefineFvMesh/dynamicRefineFvMesh.C b/src/dynamicFvMesh/dynamicRefineFvMesh/dynamicRefineFvMesh.C
index 3025a5cc91..93559ee08c 100644
--- a/src/dynamicFvMesh/dynamicRefineFvMesh/dynamicRefineFvMesh.C
+++ b/src/dynamicFvMesh/dynamicRefineFvMesh/dynamicRefineFvMesh.C
@@ -70,7 +70,7 @@ Foam::label Foam::dynamicRefineFvMesh::count
return n;
}
-q
+
void Foam::dynamicRefineFvMesh::calculateProtectedCells
(
diff --git a/src/fvOptions/Make/files b/src/fvOptions/Make/files
index d45448f27a..dfe1ac5043 100644
--- a/src/fvOptions/Make/files
+++ b/src/fvOptions/Make/files
@@ -12,6 +12,8 @@ $(generalSources)/semiImplicitSource/semiImplicitSource.C
derivedSources=sources/derived
$(derivedSources)/actuationDiskSource/actuationDiskSource.C
+$(derivedSources)/enthalpyPorositySource/enthalpyPorositySource.C
+$(derivedSources)/enthalpyPorositySource/enthalpyPorositySourceIO.C
$(derivedSources)/explicitPorositySource/explicitPorositySource.C
$(derivedSources)/MRFSource/MRFSource.C
$(derivedSources)/pressureGradientExplicitSource/pressureGradientExplicitSource.C
diff --git a/src/fvOptions/sources/derived/enthalpyPorositySource/enthalpyPorositySource.C b/src/fvOptions/sources/derived/enthalpyPorositySource/enthalpyPorositySource.C
new file mode 100644
index 0000000000..2a2d109941
--- /dev/null
+++ b/src/fvOptions/sources/derived/enthalpyPorositySource/enthalpyPorositySource.C
@@ -0,0 +1,428 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
+ \\ / O peration |
+ \\ / A nd | Copyright (C) 2013-2014 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 "enthalpyPorositySource.H"
+#include "fvMatrices.H"
+#include "specie.H"
+#include "basicThermo.H"
+#include "uniformDimensionedFields.H"
+#include "fixedValueFvPatchFields.H"
+#include "zeroGradientFvPatchFields.H"
+#include "addToRunTimeSelectionTable.H"
+#include "fvcDdt.H"
+#include "fvcDiv.H"
+
+// * * * * * * * * * * * * * Static Member Functions * * * * * * * * * * * * //
+
+namespace Foam
+{
+ template<>
+ const char* NamedEnum
+ <
+ fv::enthalpyPorositySource::thermoMode,
+ 2
+ >::names[] =
+ {
+ "thermo",
+ "lookup"
+ };
+
+ namespace fv
+ {
+ defineTypeNameAndDebug(enthalpyPorositySource, 0);
+
+ addToRunTimeSelectionTable
+ (
+ option,
+ enthalpyPorositySource,
+ dictionary
+ );
+ }
+}
+
+const Foam::NamedEnum
+ Foam::fv::enthalpyPorositySource::thermoModeTypeNames_;
+
+
+// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
+
+Foam::wordList Foam::fv::enthalpyPorositySource::alpha1BoundaryTypes() const
+{
+ const volScalarField& T = mesh_.lookupObject(TName_);
+
+ wordList bTypes(T.boundaryField().size());
+
+ forAll(bTypes, patchI)
+ {
+ const fvPatchField& pf = T.boundaryField()[patchI];
+ if (isA(pf))
+ {
+ bTypes[patchI] = fixedValueFvPatchScalarField::typeName;
+ }
+ else
+ {
+ bTypes[patchI] = zeroGradientFvPatchScalarField::typeName;
+ }
+ }
+
+ return bTypes;
+}
+
+
+bool Foam::fv::enthalpyPorositySource::solveField(const word& fieldName) const
+{
+ bool result = true;
+
+ switch (mode_)
+ {
+ case mdThermo:
+ {
+ const basicThermo& thermo =
+ mesh_.lookupObject("thermophysicalProperties");
+
+ if (fieldName != thermo.he().name())
+ {
+ result = false;
+ }
+ break;
+ }
+ case mdLookup:
+ {
+ if (fieldName != TName_)
+ {
+ result = false;
+ }
+ break;
+ }
+ default:
+ {
+ FatalErrorIn
+ (
+ "bool Foam::fv::enthalpyPorositySource::solveField"
+ "("
+ "const word&"
+ ") const"
+ )
+ << "Unhandled thermo mode: " << thermoModeTypeNames_[mode_]
+ << abort(FatalError);
+ }
+ }
+
+ return result;
+}
+
+
+Foam::tmp Foam::fv::enthalpyPorositySource::rho() const
+{
+ switch (mode_)
+ {
+ case mdThermo:
+ {
+ const basicThermo& thermo =
+ mesh_.lookupObject("thermophysicalProperties");
+
+ return thermo.rho();
+ break;
+ }
+ case mdLookup:
+ {
+ if (rhoName_ == "rhoRef")
+ {
+ scalar rhoRef(readScalar(coeffs_.lookup("rhoRef")));
+
+ return tmp
+ (
+ new volScalarField
+ (
+ IOobject
+ (
+ name_ + ":rho",
+ mesh_.time().timeName(),
+ mesh_,
+ IOobject::NO_READ,
+ IOobject::NO_WRITE
+ ),
+ mesh_,
+ dimensionedScalar("rho", dimDensity, rhoRef),
+ zeroGradientFvPatchScalarField::typeName
+ )
+ );
+ }
+ else
+ {
+ return mesh_.lookupObject(rhoName_);
+ }
+
+ break;
+ }
+ default:
+ {
+ FatalErrorIn
+ (
+ "Foam::tmp "
+ "Foam::fv::enthalpyPorositySource::rho() const"
+ )
+ << "Unhandled thermo mode: " << thermoModeTypeNames_[mode_]
+ << abort(FatalError);
+ }
+ }
+
+ return tmp(NULL);
+}
+
+
+Foam::vector Foam::fv::enthalpyPorositySource::g() const
+{
+ if (mesh_.foundObject("g"))
+ {
+ const uniformDimensionedVectorField& value =
+ mesh_.lookupObject("g");
+ return value.value();
+ }
+ else
+ {
+ return coeffs_.lookup("g");
+ }
+}
+
+
+void Foam::fv::enthalpyPorositySource::update()
+{
+ if (curTimeIndex_ == mesh_.time().timeIndex())
+ {
+ return;
+ }
+
+ const volScalarField& T = mesh_.lookupObject(TName_);
+
+ forAll(cells_, i)
+ {
+ label cellI = cells_[i];
+
+ scalar alpha1New = 0.0;
+ scalar Tc = T[cellI];
+
+ if (Tc > Tliquidus_)
+ {
+ alpha1New = 1.0;
+ }
+ else if (Tc < Tsolidus_)
+ {
+ alpha1New = 0.0;
+ }
+ else
+ {
+ // lever rule
+ alpha1New = (Tc - Tsolidus_)/(Tliquidus_ - Tsolidus_);
+ }
+
+ alpha1New = (1.0 - relax_)*alpha1_[cellI] + relax_*alpha1New;
+
+ dAlpha1_[i] = alpha1New - alpha1_[cellI];
+
+ alpha1_[cellI] = alpha1New;
+
+ curTimeIndex_ = mesh_.time().timeIndex();
+ }
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
+
+Foam::fv::enthalpyPorositySource::enthalpyPorositySource
+(
+ const word& sourceName,
+ const word& modelType,
+ const dictionary& dict,
+ const fvMesh& mesh
+)
+:
+ option(sourceName, modelType, dict, mesh),
+ Tliquidus_(readScalar(coeffs_.lookup("Tliquidus"))),
+ Tsolidus_(readScalar(coeffs_.lookup("Tsolidus"))),
+ L_(readScalar(coeffs_.lookup("L"))),
+ relax_(coeffs_.lookupOrDefault("relax", 0.9)),
+ mode_(thermoModeTypeNames_.read(coeffs_.lookup("thermoMode"))),
+ rhoName_(coeffs_.lookupOrDefault("rhoName", "rho")),
+ TName_(coeffs_.lookupOrDefault("TName", "T")),
+ UName_(coeffs_.lookupOrDefault("UName", "U")),
+ 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.0),
+ alpha1BoundaryTypes()
+ ),
+ dAlpha1_(cells_.size(), 0.0),
+ curTimeIndex_(-1)
+{
+ fieldNames_.setSize(1, "source");
+ applied_.setSize(1, false);
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
+
+bool Foam::fv::enthalpyPorositySource::alwaysApply() const
+{
+ return true;
+}
+
+
+void Foam::fv::enthalpyPorositySource::addSup
+(
+ fvMatrix& eqn,
+ const label fieldI
+)
+{
+ if (!solveField(eqn.psi().name()))
+ {
+ return;
+ }
+
+ if (debug)
+ {
+ Info<< type() << ": applying source to " << eqn.psi().name() << endl;
+ }
+
+ update();
+
+ volScalarField dH
+ (
+ IOobject
+ (
+ name_ + ":dH",
+ mesh_.time().timeName(),
+ mesh_,
+ IOobject::NO_READ,
+ IOobject::NO_WRITE
+ ),
+ mesh_,
+ dimensionedScalar("dH", dimEnergy/dimMass, 0.0)
+ );
+
+ scalarField& dHi = dH.internalField();
+
+ forAll(cells_, i)
+ {
+ label cellI = cells_[i];
+ dHi[cellI] = L_*dAlpha1_[i];
+ }
+
+ const volScalarField rho(this->rho());
+
+ const surfaceScalarField& phi =
+ mesh_.lookupObject("phi");
+
+ dimensionedScalar rhoScale("rhoScale", dimless, 1.0);
+
+ if
+ (
+ (phi.dimensions() == dimVolume/dimTime)
+ && (rho.dimensions() == dimDensity)
+ )
+ {
+ rhoScale.dimensions() /= dimDensity;
+ }
+
+ // contributions added to rhs
+ if (eqn.psi().dimensions() == dimTemperature)
+ {
+ dimensionedScalar Cp
+ (
+ "Cp",
+ dimEnergy/dimMass/dimTemperature,
+ readScalar(coeffs_.lookup("CpRef"))
+ );
+
+ eqn -=
+ fvc::ddt((rho*rhoScale*dH/Cp)())
+ + fvc::div(phi*fvc::interpolate(dH/Cp));
+ }
+ else
+ {
+ eqn -=
+ fvc::ddt((rho*rhoScale*dH)())
+ + fvc::div(phi*fvc::interpolate(dH));
+ }
+}
+
+
+void Foam::fv::enthalpyPorositySource::addSup
+(
+ fvMatrix& eqn,
+ const label fieldI
+)
+{
+ const volVectorField& U = eqn.psi();
+
+ if (U.name() != UName_)
+ {
+ return;
+ }
+
+ if (debug)
+ {
+ Info<< type() << ": applying source to " << UName_ << endl;
+ }
+
+ update();
+
+ scalarField& Sp = eqn.diag();
+ vectorField& Su = eqn.source();
+
+ const scalarField& V = mesh_.V();
+ const volScalarField rho(this->rho());
+ const volScalarField& T = mesh_.lookupObject(TName_);
+
+ vector g = this->g();
+
+ forAll(cells_, i)
+ {
+ label cellI = cells_[i];
+
+ scalar Vc = V[cellI];
+
+ scalar Tc = T[cellI];
+ scalar rhoc = rho[cellI];
+ scalar alpha1c = alpha1_[cellI];
+
+ Sp[cellI] -= Vc*rhoc*Cu_*sqr(1.0 - alpha1c)/(pow3(alpha1c) + q_);
+ Su[cellI] += Vc*rhoc*g*beta_*max(0, (Tc - Tsolidus_));
+ }
+}
+
+
+// ************************************************************************* //
diff --git a/src/fvOptions/sources/derived/enthalpyPorositySource/enthalpyPorositySource.H b/src/fvOptions/sources/derived/enthalpyPorositySource/enthalpyPorositySource.H
new file mode 100644
index 0000000000..8ee254ceb9
--- /dev/null
+++ b/src/fvOptions/sources/derived/enthalpyPorositySource/enthalpyPorositySource.H
@@ -0,0 +1,250 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
+ \\ / O peration |
+ \\ / A nd | Copyright (C) 2013-2014 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::enthalpyPorositySource
+
+Description
+ This source is designed to model the effect of solidification and melting
+ processes, e.g. windhield defrosting. The phase change occurs between the
+ limiting temperatures of Tliquidus and Tsolidus.
+
+ 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.
+
+
+ \heading Source usage
+ Example usage:
+ \verbatim
+ enthalpyPorositySourceCoeffs
+ {
+ type enthalpyPorositySource;
+ active on;
+ selectionMode cellZone;
+ cellZone iceZone;
+
+ enthalpyPorositySourceCoeffs
+ {
+ Tliquidus 288;
+ Tsolidus 268;
+ L 334000;
+ thermoMode thermo;
+ beta 50e-6;
+
+ // only for incompressible solvers:
+ // rhoName rhoRef;
+ // rhoRef 1;
+
+ // only for solvers that do not define gravity:
+ g (0 -9.81 0);
+ }
+ }
+ \endverbatim
+
+ Where:
+ \table
+ Property | Description | Required | Default value
+ Tliquidus | Temperature when liquid [K] | yes |
+ Tsolidus | Temperature when solid [K] | yes |
+ L | Latent heat of fusion [J/kg] | yes |
+ thermoMode | Thermo mode [thermo|lookup] | yes |
+ beta | Thermal expansion coefficient [1/K] | yes |
+ rhoName | Name of density field | no | rho
+ rhoRef | Reference density | no |
+ g | Accelerartion due to gravity | no |
+ \endtable
+
+SourceFiles
+ enthalpyPorositySource.C
+ enthalpyPorositySourceIO.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef enthalpyPorositySource_H
+#define enthalpyPorositySource_H
+
+#include "fvMesh.H"
+#include "volFields.H"
+#include "fvOption.H"
+#include "NamedEnum.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace fv
+{
+
+/*---------------------------------------------------------------------------*\
+ Class enthalpyPorositySource Declaration
+\*---------------------------------------------------------------------------*/
+
+class enthalpyPorositySource
+:
+ public option
+{
+public:
+
+ enum thermoMode
+ {
+ mdThermo,
+ mdLookup
+ };
+
+ static const NamedEnum thermoModeTypeNames_;
+
+
+private:
+
+ // Private data
+
+ //- Solidification initiated when T > Tliquidus_ [K]
+ scalar Tliquidus_;
+
+ //- Solidification temperature [K]
+ scalar Tsolidus_;
+
+ //- Latent heat of fusion [J/kg]
+ scalar L_;
+
+ //- Phase fraction under-relaxation coefficient
+ scalar relax_;
+
+ //- Thermodynamics mode
+ thermoMode mode_;
+
+ //- Name of density field - default = "rho" (optional)
+ word rhoName_;
+
+ //- Name of temperature field - default = "T" (optional)
+ word TName_;
+
+ //- Name of velocity field - default = "U" (optional)
+ word UName_;
+
+ //- 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_;
+
+ //- Change in phase fraction indicator field
+ scalarField dAlpha1_;
+
+ //- Current time index (used for updating)
+ label curTimeIndex_;
+
+
+ // Private Member Functions
+
+ //- Return the list of alpha1 boundary types
+ wordList alpha1BoundaryTypes() const;
+
+ //- Flag to indicate whether to solve for given field
+ bool solveField(const word& fieldName) const;
+
+ //- Return the density field
+ tmp rho() const;
+
+ //- Return the gravity vector
+ vector g() const;
+
+ //- Update the model
+ void update();
+
+ //- Disallow default bitwise copy construct
+ enthalpyPorositySource(const enthalpyPorositySource&);
+
+ //- Disallow default bitwise assignment
+ void operator=(const enthalpyPorositySource&);
+
+
+public:
+
+ //- Runtime type information
+ TypeName("enthalpyPorositySource");
+
+
+ // Constructors
+
+ //- Construct from explicit source name and mesh
+ enthalpyPorositySource
+ (
+ 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;
+
+
+ // Evaluate
+
+ //- 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);
+
+
+ // 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
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/fvOptions/sources/derived/enthalpyPorositySource/enthalpyPorositySourceIO.C b/src/fvOptions/sources/derived/enthalpyPorositySource/enthalpyPorositySourceIO.C
new file mode 100644
index 0000000000..7b5d0c02e1
--- /dev/null
+++ b/src/fvOptions/sources/derived/enthalpyPorositySource/enthalpyPorositySourceIO.C
@@ -0,0 +1,69 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
+ \\ / O peration |
+ \\ / A nd | Copyright (C) 2013-2014 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 "enthalpyPorositySource.H"
+
+// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
+
+void Foam::fv::enthalpyPorositySource::writeData(Ostream& os) const
+{
+ os << indent << name_ << endl;
+ dict_.write(os);
+}
+
+
+bool Foam::fv::enthalpyPorositySource::read(const dictionary& dict)
+{
+ if (option::read(dict))
+ {
+ coeffs_.lookup("Tliquidus") >> Tliquidus_;
+ coeffs_.lookup("Tsolidus") >> Tsolidus_;
+ coeffs_.lookup("L") >> L_;
+
+ coeffs_.readIfPresent("relax", relax_);
+
+ mode_ = thermoModeTypeNames_.read(coeffs_.lookup("thermoMode"));
+
+ coeffs_.readIfPresent("rhoName", rhoName_);
+ 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/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/refinementParameters/wordPairHashTable.H b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/refinementParameters/wordPairHashTable.H
index a2a7ea414d..26ccbbb279 100644
--- a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/refinementParameters/wordPairHashTable.H
+++ b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/refinementParameters/wordPairHashTable.H
@@ -45,7 +45,7 @@ Description
namespace Foam
{
- typedef HashTable, typename FixedList::Hash<> >
+ typedef HashTable, FixedList::Hash<> >
wordPairHashTable;
}
diff --git a/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinement.H b/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinement.H
index 7ae5d5ac29..2a01cd9a87 100644
--- a/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinement.H
+++ b/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinement.H
@@ -571,7 +571,7 @@ private:
const label ownZone,
const label neiZone,
wordPairHashTable& zonesToFaceZone,
- HashTable >&
+ HashTable >&
) const;
//- Remove any loose standing cells
diff --git a/src/thermophysicalModels/basic/psiThermo/psiThermos.C b/src/thermophysicalModels/basic/psiThermo/psiThermos.C
index add4c88cea..f73f2e4ca0 100644
--- a/src/thermophysicalModels/basic/psiThermo/psiThermos.C
+++ b/src/thermophysicalModels/basic/psiThermo/psiThermos.C
@@ -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
@@ -28,6 +28,7 @@ License
#include "specie.H"
#include "perfectGas.H"
+#include "PengRobinsonGas.H"
#include "hConstThermo.H"
#include "eConstThermo.H"
#include "janafThermo.H"
@@ -38,6 +39,10 @@ License
#include "constTransport.H"
#include "sutherlandTransport.H"
+#include "icoPolynomial.H"
+#include "hPolynomialThermo.H"
+#include "polynomialTransport.H"
+
#include "hePsiThermo.H"
#include "pureMixture.H"
@@ -85,6 +90,43 @@ makeThermo
);
+makeThermo
+(
+ psiThermo,
+ hePsiThermo,
+ pureMixture,
+ sutherlandTransport,
+ sensibleEnthalpy,
+ hConstThermo,
+ PengRobinsonGas,
+ specie
+);
+
+makeThermo
+(
+ psiThermo,
+ hePsiThermo,
+ pureMixture,
+ polynomialTransport,
+ sensibleEnthalpy,
+ hPolynomialThermo,
+ PengRobinsonGas,
+ specie
+);
+
+makeThermo
+(
+ psiThermo,
+ hePsiThermo,
+ pureMixture,
+ polynomialTransport,
+ sensibleEnthalpy,
+ janafThermo,
+ PengRobinsonGas,
+ specie
+);
+
+
/* * * * * * * * * * * * * * Internal-energy-based * * * * * * * * * * * * * */
makeThermo
diff --git a/src/thermophysicalModels/basic/rhoThermo/rhoThermos.C b/src/thermophysicalModels/basic/rhoThermo/rhoThermos.C
index 35d49cca22..874cbe50d3 100644
--- a/src/thermophysicalModels/basic/rhoThermo/rhoThermos.C
+++ b/src/thermophysicalModels/basic/rhoThermo/rhoThermos.C
@@ -31,6 +31,7 @@ License
#include "incompressiblePerfectGas.H"
#include "rhoConst.H"
#include "perfectFluid.H"
+#include "PengRobinsonGas.H"
#include "adiabaticPerfectFluid.H"
#include "hConstThermo.H"
#include "janafThermo.H"
@@ -175,6 +176,41 @@ makeThermo
specie
);
+makeThermo
+(
+ rhoThermo,
+ heRhoThermo,
+ pureMixture,
+ sutherlandTransport,
+ sensibleEnthalpy,
+ hConstThermo,
+ PengRobinsonGas,
+ specie
+);
+
+makeThermo
+(
+ rhoThermo,
+ heRhoThermo,
+ pureMixture,
+ polynomialTransport,
+ sensibleEnthalpy,
+ hPolynomialThermo,
+ PengRobinsonGas,
+ specie
+);
+
+makeThermo
+(
+ rhoThermo,
+ heRhoThermo,
+ pureMixture,
+ polynomialTransport,
+ sensibleEnthalpy,
+ janafThermo,
+ PengRobinsonGas,
+ specie
+);
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
diff --git a/src/thermophysicalModels/specie/equationOfState/PengRobinsonGas/PengRobinsonGas.C b/src/thermophysicalModels/specie/equationOfState/PengRobinsonGas/PengRobinsonGas.C
new file mode 100644
index 0000000000..e2b99d08a7
--- /dev/null
+++ b/src/thermophysicalModels/specie/equationOfState/PengRobinsonGas/PengRobinsonGas.C
@@ -0,0 +1,88 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
+ \\ / O peration |
+ \\ / A nd | Copyright (C) 2014 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 "PengRobinsonGas.H"
+#include "IOstreams.H"
+
+// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
+
+template
+Foam::PengRobinsonGas::PengRobinsonGas(Istream& is)
+:
+ Specie(is),
+ Tc_(readScalar(is)),
+ Pc_(readScalar(is)),
+ omega_(readScalar(is))
+{
+ is.check("PengRobinsonGas::PengRobinsonGas(Istream& is)");
+}
+
+
+template
+Foam::PengRobinsonGas::PengRobinsonGas
+(
+ const dictionary& dict
+)
+:
+ Specie(dict),
+ Tc_(readScalar(dict.subDict("equationOfState").lookup("Tc"))),
+ Pc_(readScalar(dict.subDict("equationOfState").lookup("Pc"))),
+ omega_(readScalar(dict.subDict("equationOfState").lookup("omega")))
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
+
+
+template
+void Foam::PengRobinsonGas::write(Ostream& os) const
+{
+ Specie::write(os);
+}
+
+
+// * * * * * * * * * * * * * * * Ostream Operator * * * * * * * * * * * * * //
+
+template
+Foam::Ostream& Foam::operator<<
+(
+ Ostream& os,
+ const PengRobinsonGas& pg
+)
+{
+ os << static_cast(pg)
+ << token::SPACE << pg.Tc_
+ << token::SPACE << pg.Pc_
+ << token::SPACE << pg.omega_;
+
+ os.check
+ (
+ "Ostream& operator<<(Ostream& os, const PengRobinsonGas& st)"
+ );
+ return os;
+}
+
+
+// ************************************************************************* //
diff --git a/src/thermophysicalModels/specie/equationOfState/PengRobinsonGas/PengRobinsonGas.H b/src/thermophysicalModels/specie/equationOfState/PengRobinsonGas/PengRobinsonGas.H
new file mode 100644
index 0000000000..7d1df11245
--- /dev/null
+++ b/src/thermophysicalModels/specie/equationOfState/PengRobinsonGas/PengRobinsonGas.H
@@ -0,0 +1,239 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
+ \\ / O peration |
+ \\ / A nd | Copyright (C) 2014 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::PengRobinsonGas
+
+Description
+ PengRobinsonGas gas equation of state.
+
+SourceFiles
+ PengRobinsonGasI.H
+ PengRobinsonGas.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef PengRobinsonGas_H
+#define PengRobinsonGas_H
+
+#include "autoPtr.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// Forward declaration of friend functions and operators
+
+template class PengRobinsonGas;
+
+template
+inline PengRobinsonGas operator+
+(
+ const PengRobinsonGas&,
+ const PengRobinsonGas&
+);
+
+template
+inline PengRobinsonGas operator-
+(
+ const PengRobinsonGas&,
+ const PengRobinsonGas&
+);
+
+template
+inline PengRobinsonGas operator*
+(
+ const scalar,
+ const PengRobinsonGas&
+);
+
+template
+inline PengRobinsonGas operator==
+(
+ const PengRobinsonGas&,
+ const PengRobinsonGas&
+);
+
+template
+Ostream& operator<<
+(
+ Ostream&,
+ const PengRobinsonGas&
+);
+
+
+
+/*---------------------------------------------------------------------------*\
+ Class PengRobinsonGas Declaration
+\*---------------------------------------------------------------------------*/
+
+template
+class PengRobinsonGas
+:
+ public Specie
+{
+ // Private data
+
+ //- Critical Temperature [K]
+ scalar Tc_;
+
+ //- Critical Pressure [Pa]
+ scalar Pc_;
+
+ //- Accentric factor [-]
+ scalar omega_;
+
+
+public:
+
+ // Constructors
+
+ //- Construct from components
+ inline PengRobinsonGas
+ (
+ const Specie& sp,
+ const scalar& Tc,
+ const scalar& Pc,
+ const scalar& omega
+ );
+
+ //- Construct from Istream
+ PengRobinsonGas(Istream&);
+
+ //- Construct from dictionary
+ PengRobinsonGas(const dictionary& dict);
+
+ //- Construct as named copy
+ inline PengRobinsonGas(const word& name, const PengRobinsonGas&);
+
+ //- Construct and return a clone
+ inline autoPtr clone() const;
+
+ // Selector from Istream
+ inline static autoPtr New(Istream& is);
+
+ // Selector from dictionary
+ inline static autoPtr New
+ (
+ const dictionary& dict
+ );
+
+
+ // Member functions
+
+ //- Return the instantiated type name
+ static word typeName()
+ {
+ return "PengRobinsonGas<" + word(Specie::typeName_()) + '>';
+ }
+
+ // Fundamental properties
+
+
+ //- Is the equation of state is incompressible i.e. rho != f(p)
+ static const bool incompressible = false;
+
+ //- Is the equation of state is isochoric i.e. rho = const
+ static const bool isochoric = false;
+
+ //- Return density [kg/m^3]
+ inline scalar rho(scalar p, scalar T) const;
+
+ //- Return compressibility rho/p [s^2/m^2]
+ inline scalar psi(scalar p, scalar T) const;
+
+ //- Return compression factor []
+ inline scalar Z(scalar p, scalar T) const;
+
+ //- Return (cp - cv) [J/(kmol K]
+ inline scalar cpMcv(scalar p, scalar T) const;
+
+
+ // IO
+
+ //- Write to Ostream
+ void write(Ostream& os) const;
+
+ // Member operators
+
+ inline void operator+=(const PengRobinsonGas&);
+ inline void operator-=(const PengRobinsonGas&);
+
+ inline void operator*=(const scalar);
+
+
+ // Friend operators
+
+ friend PengRobinsonGas operator+
+ (
+ const PengRobinsonGas&,
+ const PengRobinsonGas&
+ );
+
+ friend PengRobinsonGas operator-
+ (
+ const PengRobinsonGas&,
+ const PengRobinsonGas&
+ );
+
+ friend PengRobinsonGas operator*
+ (
+ const scalar s,
+ const PengRobinsonGas&
+ );
+
+ friend PengRobinsonGas operator==
+ (
+ const PengRobinsonGas&,
+ const PengRobinsonGas&
+ );
+
+
+ // Ostream Operator
+
+ friend Ostream& operator<<
+ (
+ Ostream&,
+ const PengRobinsonGas&
+ );
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#include "PengRobinsonGasI.H"
+
+#ifdef NoRepository
+# include "PengRobinsonGas.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/thermophysicalModels/specie/equationOfState/PengRobinsonGas/PengRobinsonGasI.H b/src/thermophysicalModels/specie/equationOfState/PengRobinsonGas/PengRobinsonGasI.H
new file mode 100644
index 0000000000..aae0b92fcb
--- /dev/null
+++ b/src/thermophysicalModels/specie/equationOfState/PengRobinsonGas/PengRobinsonGasI.H
@@ -0,0 +1,323 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
+ \\ / O peration |
+ \\ / A nd | Copyright (C) 2014 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 "PengRobinsonGas.H"
+#include "mathematicalConstants.H"
+
+// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
+
+template
+inline Foam::PengRobinsonGas::PengRobinsonGas
+(
+ const Specie& sp,
+ const scalar& Tc,
+ const scalar& Pc,
+ const scalar& omega
+)
+:
+ Specie(sp),
+ Tc_(Tc),
+ Pc_(Pc),
+ omega_(omega)
+{}
+
+
+// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
+
+template
+inline Foam::PengRobinsonGas::PengRobinsonGas
+(
+ const word& name,
+ const PengRobinsonGas& pg
+)
+:
+ Specie(name, pg),
+ Tc_(pg.Tc),
+ Pc_(pg.Pc),
+ omega_(pg.omega)
+{}
+
+
+template
+inline Foam::autoPtr >
+Foam::PengRobinsonGas::clone() const
+{
+ return autoPtr >
+ (
+ new PengRobinsonGas(*this)
+ );
+}
+
+
+template
+inline Foam::autoPtr >
+Foam::PengRobinsonGas::New
+(
+ Istream& is
+)
+{
+ return autoPtr >(new PengRobinsonGas(is));
+}
+
+
+template
+inline Foam::autoPtr >
+Foam::PengRobinsonGas::New
+(
+ const dictionary& dict
+)
+{
+ return autoPtr >
+ (
+ new PengRobinsonGas(dict)
+ );
+}
+
+// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
+
+template
+inline Foam::scalar Foam::PengRobinsonGas::rho
+(
+ scalar p,
+ scalar T
+) const
+{
+ scalar z = Z(p, T);
+ return p/(z*this->R()*T);
+}
+
+
+template
+inline Foam::scalar Foam::PengRobinsonGas::psi
+(
+ scalar p,
+ scalar T
+) const
+{
+ scalar z = Z(p, T);
+ return 1.0/(z*this->R()*T);
+}
+
+
+template
+inline Foam::scalar Foam::PengRobinsonGas::Z
+(
+ scalar p,
+ scalar T
+) const
+{
+ scalar a = 0.45724*sqr(this->R())*sqr(Tc_)/Pc_;
+ scalar b = 0.07780*this->R()*Tc_/Pc_;
+ scalar Tr = T/Tc_;
+
+ scalar alpha =
+ sqr
+ (
+ 1.0
+ + (0.37464 + 1.54226*omega_- 0.26992*sqr(omega_))
+ * (1.0 - sqrt(Tr))
+ );
+
+ scalar B = b*p/(this->R()*T);
+ scalar A = a*alpha*p/sqr(this->R()*T);
+
+ scalar a2 = B - 1.0;
+ scalar a1 = A - 2.0*B - 3.0*sqr(B);
+ scalar a0 = -A*B + sqr(B) + pow3(B);
+
+ scalar Q = (3.0*a1 - a2*a2)/9.0;
+ scalar Rl = (9.0*a2*a1 - 27.0*a0 - 2.0*a2*a2*a2)/54;
+
+ scalar Q3 = Q*Q*Q;
+ scalar D = Q3 + Rl*Rl;
+
+ scalar root = -1.0;
+
+ if (D <= 0.0)
+ {
+ scalar th = ::acos(Rl/sqrt(-Q3));
+ scalar qm = 2.0*sqrt(-Q);
+ scalar r1 = qm*cos(th/3.0) - a2/3.0;
+ scalar r2 = qm*cos((th + 2.0*constant::mathematical::pi)/3.0) - a2/3.0;
+ scalar r3 = qm*cos((th + 4.0*constant::mathematical::pi)/3.0) - a2/3.0;
+
+ root = max(r1, max(r2, r3));
+ }
+ else
+ {
+ // one root is real
+ scalar D05 = sqrt(D);
+ scalar S = pow(Rl + D05, 1.0/3.0);
+ scalar Tl = 0;
+ if (D05 > Rl)
+ {
+ Tl = -pow(mag(Rl - D05), 1.0/3.0);
+ }
+ else
+ {
+ Tl = pow(Rl - D05, 1.0/3.0);
+ }
+
+ root = S + Tl - a2/3.0;
+ }
+
+ return root;
+}
+
+
+template
+inline Foam::scalar Foam::PengRobinsonGas::cpMcv
+(
+ scalar p,
+ scalar T
+) const
+{
+ return this->RR*Z(p, T);
+}
+
+
+// * * * * * * * * * * * * * * * Member Operators * * * * * * * * * * * * * //
+
+template
+inline void Foam::PengRobinsonGas::operator+=
+(
+ const PengRobinsonGas& pg
+)
+{
+ scalar molr1 = this->nMoles();
+ Specie::operator+=(pg);
+
+ molr1 /= this->nMoles();
+ scalar molr2 = pg.nMoles()/this->nMoles();
+
+ Tc_ = molr1*Tc_ + molr2*pg.Tc_;
+ Pc_ = molr1*Pc_ + molr2*pg.Pc_;
+ omega_ = molr1*omega_ + molr2*pg.omega_;
+}
+
+
+template
+inline void Foam::PengRobinsonGas::operator-=
+(
+ const PengRobinsonGas& pg
+)
+{
+ scalar molr1 = this->nMoles();
+
+ Specie::operator-=(pg);
+
+ molr1 /= this->nMoles();
+ scalar molr2 = pg.nMoles()/this->nMoles();
+
+ Tc_ = molr1*Tc_ - molr2*pg.Tc_;
+ Pc_ = molr1*Pc_ - molr2*pg.Pc_;
+ omega_ = molr1*omega_ - molr2*pg.omega_;
+}
+
+
+template
+inline void Foam::PengRobinsonGas::operator*=(const scalar s)
+{
+ Specie::operator*=(s);
+}
+
+
+// * * * * * * * * * * * * * * * Friend Operators * * * * * * * * * * * * * //
+
+
+template
+Foam::PengRobinsonGas Foam::operator+
+(
+ const PengRobinsonGas& pg1,
+ const PengRobinsonGas& pg2
+)
+{
+ scalar nMoles = pg1.nMoles() + pg2.nMoles();
+ scalar molr1 = pg1.nMoles()/nMoles;
+ scalar molr2 = pg2.nMoles()/nMoles;
+
+ return PengRobinsonGas
+ (
+ static_cast(pg1)
+ + static_cast(pg2),
+ molr1*pg1.Tc_ + molr2*pg2.Tc_,
+ molr1*pg1.Pc_ + molr2*pg2.Pc_,
+ molr1*pg1.omega_ + molr2*pg2.omega_
+ );
+}
+
+
+template
+Foam::PengRobinsonGas Foam::operator-
+(
+ const PengRobinsonGas& pg1,
+ const PengRobinsonGas& pg2
+)
+{
+ scalar nMoles = pg1.nMoles() + pg2.nMoles();
+ scalar molr1 = pg1.nMoles()/nMoles;
+ scalar molr2 = pg2.nMoles()/nMoles;
+
+ return PengRobinsonGas
+ (
+ static_cast(pg1)
+ - static_cast(pg2),
+ molr1*pg1.Tc_ - molr2*pg2.Tc_,
+ molr1*pg1.Pc_ - molr2*pg2.Pc_,
+ molr1*pg1.omega_ - molr2*pg2.omega_
+ );
+}
+
+
+template
+Foam::PengRobinsonGas Foam::operator*
+(
+ const scalar s,
+ const PengRobinsonGas& pg
+)
+{
+ return PengRobinsonGas
+ (
+ s*static_cast(pg),
+ pg.Tc_,
+ pg.Pc_,
+ pg.omega_
+ );
+}
+
+
+template
+Foam::PengRobinsonGas Foam::operator==
+(
+ const PengRobinsonGas& pg1,
+ const PengRobinsonGas& pg2
+)
+{
+ return pg2 - pg1;
+}
+
+
+// ************************************************************************* //