Merge branch 'master' of /home/dm4/OpenFOAM/repositories/OpenFOAM-dev

This commit is contained in:
mattijs
2014-03-05 15:21:30 +00:00
13 changed files with 1482 additions and 5 deletions

View File

@ -1510,7 +1510,7 @@ int main(int argc, char *argv[])
(
snapDict,
motionDict,
mergePatchFaces
mergePatchFaces,
curvature,
planarAngle,
snapParams

View File

@ -70,7 +70,7 @@ Foam::label Foam::dynamicRefineFvMesh::count
return n;
}
q
void Foam::dynamicRefineFvMesh::calculateProtectedCells
(

View File

@ -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

View File

@ -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 <http://www.gnu.org/licenses/>.
\*---------------------------------------------------------------------------*/
#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::thermoMode, 2>
Foam::fv::enthalpyPorositySource::thermoModeTypeNames_;
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
Foam::wordList Foam::fv::enthalpyPorositySource::alpha1BoundaryTypes() const
{
const volScalarField& T = mesh_.lookupObject<volScalarField>(TName_);
wordList bTypes(T.boundaryField().size());
forAll(bTypes, patchI)
{
const fvPatchField<scalar>& pf = T.boundaryField()[patchI];
if (isA<fixedValueFvPatchScalarField>(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<basicThermo>("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::volScalarField> Foam::fv::enthalpyPorositySource::rho() const
{
switch (mode_)
{
case mdThermo:
{
const basicThermo& thermo =
mesh_.lookupObject<basicThermo>("thermophysicalProperties");
return thermo.rho();
break;
}
case mdLookup:
{
if (rhoName_ == "rhoRef")
{
scalar rhoRef(readScalar(coeffs_.lookup("rhoRef")));
return tmp<volScalarField>
(
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<volScalarField>(rhoName_);
}
break;
}
default:
{
FatalErrorIn
(
"Foam::tmp<Foam::volScalarField> "
"Foam::fv::enthalpyPorositySource::rho() const"
)
<< "Unhandled thermo mode: " << thermoModeTypeNames_[mode_]
<< abort(FatalError);
}
}
return tmp<volScalarField>(NULL);
}
Foam::vector Foam::fv::enthalpyPorositySource::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::enthalpyPorositySource::update()
{
if (curTimeIndex_ == mesh_.time().timeIndex())
{
return;
}
const volScalarField& T = mesh_.lookupObject<volScalarField>(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<word>("rhoName", "rho")),
TName_(coeffs_.lookupOrDefault<word>("TName", "T")),
UName_(coeffs_.lookupOrDefault<word>("UName", "U")),
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),
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<scalar>& 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<surfaceScalarField>("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<vector>& 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<volScalarField>(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_));
}
}
// ************************************************************************* //

View File

@ -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 <http://www.gnu.org/licenses/>.
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<thermoMode, 2> 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<volScalarField> 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<scalar>& eqn, const label fieldI);
//- Add implicit contribution to momentum equation
virtual void addSup(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
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

View File

@ -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 <http://www.gnu.org/licenses/>.
\*---------------------------------------------------------------------------*/
#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;
}
// ************************************************************************* //

View File

@ -45,7 +45,7 @@ Description
namespace Foam
{
typedef HashTable<word, Pair<word>, typename FixedList<word, 2>::Hash<> >
typedef HashTable<word, Pair<word>, FixedList<word, 2>::Hash<> >
wordPairHashTable;
}

View File

@ -571,7 +571,7 @@ private:
const label ownZone,
const label neiZone,
wordPairHashTable& zonesToFaceZone,
HashTable<word, labelPair, typename labelPair::Hash<> >&
HashTable<word, labelPair, labelPair::Hash<> >&
) const;
//- Remove any loose standing cells

View File

@ -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

View File

@ -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
);
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

View File

@ -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 <http://www.gnu.org/licenses/>.
\*---------------------------------------------------------------------------*/
#include "PengRobinsonGas.H"
#include "IOstreams.H"
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
template<class Specie>
Foam::PengRobinsonGas<Specie>::PengRobinsonGas(Istream& is)
:
Specie(is),
Tc_(readScalar(is)),
Pc_(readScalar(is)),
omega_(readScalar(is))
{
is.check("PengRobinsonGas<Specie>::PengRobinsonGas(Istream& is)");
}
template<class Specie>
Foam::PengRobinsonGas<Specie>::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<class Specie>
void Foam::PengRobinsonGas<Specie>::write(Ostream& os) const
{
Specie::write(os);
}
// * * * * * * * * * * * * * * * Ostream Operator * * * * * * * * * * * * * //
template<class Specie>
Foam::Ostream& Foam::operator<<
(
Ostream& os,
const PengRobinsonGas<Specie>& pg
)
{
os << static_cast<const Specie&>(pg)
<< token::SPACE << pg.Tc_
<< token::SPACE << pg.Pc_
<< token::SPACE << pg.omega_;
os.check
(
"Ostream& operator<<(Ostream& os, const PengRobinsonGas<Specie>& st)"
);
return os;
}
// ************************************************************************* //

View File

@ -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 <http://www.gnu.org/licenses/>.
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 Specie> class PengRobinsonGas;
template<class Specie>
inline PengRobinsonGas<Specie> operator+
(
const PengRobinsonGas<Specie>&,
const PengRobinsonGas<Specie>&
);
template<class Specie>
inline PengRobinsonGas<Specie> operator-
(
const PengRobinsonGas<Specie>&,
const PengRobinsonGas<Specie>&
);
template<class Specie>
inline PengRobinsonGas<Specie> operator*
(
const scalar,
const PengRobinsonGas<Specie>&
);
template<class Specie>
inline PengRobinsonGas<Specie> operator==
(
const PengRobinsonGas<Specie>&,
const PengRobinsonGas<Specie>&
);
template<class Specie>
Ostream& operator<<
(
Ostream&,
const PengRobinsonGas<Specie>&
);
/*---------------------------------------------------------------------------*\
Class PengRobinsonGas Declaration
\*---------------------------------------------------------------------------*/
template<class Specie>
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<PengRobinsonGas> clone() const;
// Selector from Istream
inline static autoPtr<PengRobinsonGas> New(Istream& is);
// Selector from dictionary
inline static autoPtr<PengRobinsonGas> 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+ <Specie>
(
const PengRobinsonGas&,
const PengRobinsonGas&
);
friend PengRobinsonGas operator- <Specie>
(
const PengRobinsonGas&,
const PengRobinsonGas&
);
friend PengRobinsonGas operator* <Specie>
(
const scalar s,
const PengRobinsonGas&
);
friend PengRobinsonGas operator== <Specie>
(
const PengRobinsonGas&,
const PengRobinsonGas&
);
// Ostream Operator
friend Ostream& operator<< <Specie>
(
Ostream&,
const PengRobinsonGas&
);
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#include "PengRobinsonGasI.H"
#ifdef NoRepository
# include "PengRobinsonGas.C"
#endif
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

View File

@ -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 <http://www.gnu.org/licenses/>.
\*---------------------------------------------------------------------------*/
#include "PengRobinsonGas.H"
#include "mathematicalConstants.H"
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
template<class Specie>
inline Foam::PengRobinsonGas<Specie>::PengRobinsonGas
(
const Specie& sp,
const scalar& Tc,
const scalar& Pc,
const scalar& omega
)
:
Specie(sp),
Tc_(Tc),
Pc_(Pc),
omega_(omega)
{}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
template<class Specie>
inline Foam::PengRobinsonGas<Specie>::PengRobinsonGas
(
const word& name,
const PengRobinsonGas& pg
)
:
Specie(name, pg),
Tc_(pg.Tc),
Pc_(pg.Pc),
omega_(pg.omega)
{}
template<class Specie>
inline Foam::autoPtr<Foam::PengRobinsonGas <Specie> >
Foam::PengRobinsonGas<Specie>::clone() const
{
return autoPtr<PengRobinsonGas<Specie> >
(
new PengRobinsonGas<Specie>(*this)
);
}
template<class Specie>
inline Foam::autoPtr<Foam::PengRobinsonGas<Specie> >
Foam::PengRobinsonGas<Specie>::New
(
Istream& is
)
{
return autoPtr<PengRobinsonGas<Specie> >(new PengRobinsonGas<Specie>(is));
}
template<class Specie>
inline Foam::autoPtr<Foam::PengRobinsonGas<Specie> >
Foam::PengRobinsonGas<Specie>::New
(
const dictionary& dict
)
{
return autoPtr<PengRobinsonGas<Specie> >
(
new PengRobinsonGas<Specie>(dict)
);
}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
template<class Specie>
inline Foam::scalar Foam::PengRobinsonGas<Specie>::rho
(
scalar p,
scalar T
) const
{
scalar z = Z(p, T);
return p/(z*this->R()*T);
}
template<class Specie>
inline Foam::scalar Foam::PengRobinsonGas<Specie>::psi
(
scalar p,
scalar T
) const
{
scalar z = Z(p, T);
return 1.0/(z*this->R()*T);
}
template<class Specie>
inline Foam::scalar Foam::PengRobinsonGas<Specie>::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<class Specie>
inline Foam::scalar Foam::PengRobinsonGas<Specie>::cpMcv
(
scalar p,
scalar T
) const
{
return this->RR*Z(p, T);
}
// * * * * * * * * * * * * * * * Member Operators * * * * * * * * * * * * * //
template<class Specie>
inline void Foam::PengRobinsonGas<Specie>::operator+=
(
const PengRobinsonGas<Specie>& 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<class Specie>
inline void Foam::PengRobinsonGas<Specie>::operator-=
(
const PengRobinsonGas<Specie>& 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<class Specie>
inline void Foam::PengRobinsonGas<Specie>::operator*=(const scalar s)
{
Specie::operator*=(s);
}
// * * * * * * * * * * * * * * * Friend Operators * * * * * * * * * * * * * //
template<class Specie>
Foam::PengRobinsonGas<Specie> Foam::operator+
(
const PengRobinsonGas<Specie>& pg1,
const PengRobinsonGas<Specie>& pg2
)
{
scalar nMoles = pg1.nMoles() + pg2.nMoles();
scalar molr1 = pg1.nMoles()/nMoles;
scalar molr2 = pg2.nMoles()/nMoles;
return PengRobinsonGas<Specie>
(
static_cast<const Specie&>(pg1)
+ static_cast<const Specie&>(pg2),
molr1*pg1.Tc_ + molr2*pg2.Tc_,
molr1*pg1.Pc_ + molr2*pg2.Pc_,
molr1*pg1.omega_ + molr2*pg2.omega_
);
}
template<class Specie>
Foam::PengRobinsonGas<Specie> Foam::operator-
(
const PengRobinsonGas<Specie>& pg1,
const PengRobinsonGas<Specie>& pg2
)
{
scalar nMoles = pg1.nMoles() + pg2.nMoles();
scalar molr1 = pg1.nMoles()/nMoles;
scalar molr2 = pg2.nMoles()/nMoles;
return PengRobinsonGas<Specie>
(
static_cast<const Specie&>(pg1)
- static_cast<const Specie&>(pg2),
molr1*pg1.Tc_ - molr2*pg2.Tc_,
molr1*pg1.Pc_ - molr2*pg2.Pc_,
molr1*pg1.omega_ - molr2*pg2.omega_
);
}
template<class Specie>
Foam::PengRobinsonGas<Specie> Foam::operator*
(
const scalar s,
const PengRobinsonGas<Specie>& pg
)
{
return PengRobinsonGas<Specie>
(
s*static_cast<const Specie&>(pg),
pg.Tc_,
pg.Pc_,
pg.omega_
);
}
template<class Specie>
Foam::PengRobinsonGas<Specie> Foam::operator==
(
const PengRobinsonGas<Specie>& pg1,
const PengRobinsonGas<Specie>& pg2
)
{
return pg2 - pg1;
}
// ************************************************************************* //