exaFOAM/foam-extend4.1-coherent-io
Archived
1
1
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity

Merged Christian's update

Browse source
This commit is contained in:
Henrik Rusche 2016-06-06 14:10:12 +02:00
commit 2477805d86
121 changed files with 4601 additions and 4261 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

4
applications/solvers/compressible/realFluidPisoFoam/createFields.H
View file

@ -11,8 +11,6 @@
const volScalarField& psi = thermo.psi();
const volScalarField& drhodh = thermo.drhodh();
bool realFluid=mesh.solutionDict().subDict("PISO").lookupOrDefault<bool>("realFluid",false);
volScalarField rho
(
IOobject
@ -42,7 +40,6 @@
# include "compressibleCreatePhi.H"
Info<< "Creating turbulence model\n" << endl;
autoPtr<compressible::turbulenceModel> turbulence
(
@ -55,7 +52,6 @@
)
);
Info<< "Creating field DpDt\n" << endl;
volScalarField DpDt =
fvc::DDt(surfaceScalarField("phiU", phi/fvc::interpolate(rho)), p);

61
applications/solvers/compressible/realFluidPisoFoam/pEqn.H
View file

@ -3,8 +3,6 @@ rho = thermo.rho();
volScalarField rUA = 1.0/UEqn.A();
U = rUA*UEqn.H();
if (realFluid)
{
phi =
fvc::interpolate(rho)*
(
@ -29,65 +27,6 @@ if (realFluid)
phi += pEqn.flux();
}
}
}
else
{
if (piso.transonic())
{
surfaceScalarField phid
(
"phid",
fvc::interpolate(psi)
*(
(fvc::interpolate(U) & mesh.Sf())
+ fvc::ddtPhiCorr(rUA, rho, U, phi)
)
);
while (piso.correctNonOrthogonal())
{
fvScalarMatrix pEqn
(
fvm::ddt(psi, p)
+ fvm::div(phid, p)
- fvm::laplacian(rho*rUA, p)
);
pEqn.solve();
if (piso.finalNonOrthogonalIter())
{
phi == pEqn.flux();
}
}
}
else
{
phi =
fvc::interpolate(rho)*
(
(fvc::interpolate(U) & mesh.Sf())
+ fvc::ddtPhiCorr(rUA, rho, U, phi)
);
while (piso.correctNonOrthogonal())
{
fvScalarMatrix pEqn
(
fvm::ddt(psi, p)
+ fvc::div(phi)
- fvm::laplacian(rho*rUA, p)
);
pEqn.solve();
if (piso.finalNonOrthogonalIter())
{
phi += pEqn.flux();
}
}
}
}
#include "rhoEqn.H"
#include "compressibleContinuityErrs.H"

8
applications/solvers/compressible/realFluidPisoFoam/realFluidPisoFoam.C
View file

@ -23,13 +23,13 @@ License
Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
Application
rhoPisoFoam
realFluidPisoFoam
Description
Transient PISO solver for compressible, laminar or turbulent flow.
Transient PISO solver for compressible, laminar or turbulent flow
of real fluids e.g. real gases (cubic equations of state)
CL: rhoPisoFoam with a changed pressure equation for non-perfect gas fluids
CL: see realFluid flag
Solver cannot be used with perfect gas library
\*---------------------------------------------------------------------------*/

15
src/thermophysicalModels/basic/IAPWS_Waterproperties/IAPWSThermo/IAPWS-IF97.C
View file

@ -27,8 +27,8 @@ Christian Lucas
Institut für Thermodynamik
Technische Universität Braunschweig
Germany
\*---------------------------------------------------------------------------*/
\*---------------------------------------------------------------------------*/
#include "IAPWS-IF97.H"
@ -55,6 +55,7 @@ void Foam::calculateProperties_ph
calculateProperties_h(S,p,h,T,rho,psi,drhodh,mu,alpha,x);
}
//CL: calculated all (minimal) needed properties + the vapor mass fraction for a given pressure and enthalpy
void Foam::calculateProperties_ph
(
@ -75,6 +76,7 @@ void Foam::calculateProperties_ph
calculateProperties_h(S,p,h,T,rho,psi,drhodh,mu,alpha,x);
}
//CL: calculated all (minimal) needed properties for a given pressure and temperature
void Foam::calculateProperties_pT
(
@ -98,6 +100,7 @@ void Foam::calculateProperties_pT
calculateProperties_h(S,p,h,T,rho,psi,drhodh,mu,alpha,x);
}
//CL: calculated all (minimal) needed properties + the vapor mass fraction for a given pressure and temperature
void Foam::calculateProperties_pT
(
@ -305,24 +308,28 @@ void Foam::calculateProperties_h
}
}
//CL: returns density for given pressure and temperature
Foam::scalar Foam::rho_pT(scalar p,scalar T)
{
return 1/freesteam_v(freesteam_set_pT(p,T));
}
//CL: returns density for given pressure and enthalpy
Foam::scalar Foam::rho_ph(scalar p,scalar h)
{
return 1/freesteam_v(freesteam_set_ph(p,h));
}
//CL: returns Cp(heat capacity @ contant pressure) for given pressure and temperature
Foam::scalar Foam::cp_pT(scalar p,scalar T)
{
return freesteam_cp(freesteam_set_pT(p,T));
}
//CL: returns Cp(heat capacity @ contant pressure) for given pressure and enthalpy
Foam::scalar Foam::cp_ph(scalar p,scalar h)
{
@ -335,18 +342,21 @@ Foam::scalar Foam::cv_pT(scalar p,scalar T)
return freesteam_cv(freesteam_set_pT(p,T));
}
//CL: returns Cv (heat capacity @ contant volume) for given pressure and enthalpy
Foam::scalar Foam::cv_ph(scalar p,scalar h)
{
return freesteam_cv(freesteam_set_ph(p,h));
}
//CL: returns enthalpy for given pressure and temperature
Foam::scalar Foam::h_pT(scalar p,scalar T)
{
return freesteam_h(freesteam_set_pT(p,T));
}
//CL: returns temperature for given pressure and enthalpy
Foam::scalar Foam::T_ph(scalar p,scalar h)
{
@ -430,7 +440,6 @@ Foam::scalar Foam::psiH(SteamState S)
scalar dvldp,dvvdp,dhldp,dhvdp;
scalar dpdT,dvdp,dxdp;
SteamState Sl,Sv;
rho=1/freesteam_region4_v_Tx(S.R4.T,S.R4.x);
@ -539,7 +548,6 @@ Foam::scalar Foam::drhodh(SteamState S)
}
else if (region==3)
{
scalar gamma,cv,p;
p=freesteam_region3_p_rhoT(S.R3.rho,S.R3.T);
@ -557,7 +565,6 @@ Foam::scalar Foam::drhodh(SteamState S)
}
else if (region==4)
{
scalar vv,vl,hl,hv,p;
SteamState Sl,Sv;

1
src/thermophysicalModels/basic/IAPWS_Waterproperties/IAPWSThermo/IAPWS-IF97.H
View file

@ -37,6 +37,7 @@ Christian Lucas
Institut für Thermodynamik
Technische Universität Braunschweig
Germany
\*---------------------------------------------------------------------------*/
#ifndef IAPWSIF97_H
#define IAPWSIF97_H

2
src/thermophysicalModels/basic/IAPWS_Waterproperties/IAPWSThermo/IAPWSThermo.C
View file

@ -211,14 +211,12 @@ Foam::IAPWSThermo::IAPWSThermo
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
Foam::IAPWSThermo::~IAPWSThermo()
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
void Foam::IAPWSThermo::correct()
{
if (debug)

1
src/thermophysicalModels/basic/IAPWS_Waterproperties/IAPWSThermo/IAPWSThermo.H
View file

@ -43,7 +43,6 @@ Institut für Thermodynamik
Technische Universität Braunschweig
Germany
\*---------------------------------------------------------------------------*/
#ifndef IAPWSThermo_H

20
src/thermophysicalModels/basic/IAPWS_Waterproperties/IAPWSThermo/IAPWSThermos.C
View file

@ -30,17 +30,9 @@ Germany
\*---------------------------------------------------------------------------*/
#include "makeBasicPsiThermo.H"
#include "IAPWSThermo.H"
// including dummy classes --> this classes do nothing
// except satisfy the template structure
//#include "dummyEqnOfState.H"
//#include "dummyThermo.H"
//#include "dummyH.H"
//#include "dummyTransport.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
@ -48,18 +40,6 @@ namespace Foam
/* * * * * * * * * * * * * * * private static data * * * * * * * * * * * * * */
/*
makeBasicRealGasThermo
(
IAPWSThermo,
pureMixture,
dummyTransport,
dummyThermo,
dummyH,
dummyEqnOfState
);
*/
makeBasicExternalLibraryBasedThermo
(
IAPWSThermo

93
src/thermophysicalModels/basic/mixtures/basicMixture/basicMixtures.C
View file

@ -33,12 +33,19 @@ Description
#include "makeBasicMixture.H"
#include "perfectGas.H"
#include "redlichKwong.H"
#include "pengRobinson.H"
#include "aungierRedlichKwong.H"
#include "soaveRedlichKwong.H"
#include "eConstThermo.H"
#include "hConstThermo.H"
#include "janafThermo.H"
#include "nasaHeatCapacityPolynomial.H"
#include "constantHeatCapacity.H"
#include "specieThermo.H"
#include "realGasSpecieThermo.H"
#include "constTransport.H"
#include "sutherlandTransport.H"
@ -49,13 +56,6 @@ Description
#include "addToRunTimeSelectionTable.H"
#include "redlichKwong.H"
#include "pengRobinson.H"
#include "aungierRedlichKwong.H"
#include "soaveRedlichKwong.H"
#include "realGasSpecieThermo.H"
#include "nasaHeatCapacityPolynomial.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
@ -118,7 +118,6 @@ makeBasicRealFluidMixture
redlichKwong
);
makeBasicRealFluidMixture
(
pureMixture,
@ -137,7 +136,6 @@ makeBasicRealFluidMixture
aungierRedlichKwong
);
makeBasicRealFluidMixture
(
pureMixture,
@ -147,8 +145,6 @@ makeBasicRealFluidMixture
soaveRedlichKwong
);
makeBasicRealFluidMixture
(
pureMixture,
@ -158,7 +154,6 @@ makeBasicRealFluidMixture
redlichKwong
);
makeBasicRealFluidMixture
(
pureMixture,
@ -168,7 +163,6 @@ makeBasicRealFluidMixture
pengRobinson
);
makeBasicRealFluidMixture
(
pureMixture,
@ -178,7 +172,6 @@ makeBasicRealFluidMixture
aungierRedlichKwong
);
makeBasicRealFluidMixture
(
pureMixture,
@ -188,6 +181,78 @@ makeBasicRealFluidMixture
soaveRedlichKwong
);
makeBasicRealFluidMixture
(
pureMixture,
sutherlandTransport,
realGasSpecieThermo,
constantHeatCapacity,
redlichKwong
);
makeBasicRealFluidMixture
(
pureMixture,
sutherlandTransport,
realGasSpecieThermo,
constantHeatCapacity,
pengRobinson
);
makeBasicRealFluidMixture
(
pureMixture,
sutherlandTransport,
realGasSpecieThermo,
constantHeatCapacity,
aungierRedlichKwong
);
makeBasicRealFluidMixture
(
pureMixture,
sutherlandTransport,
realGasSpecieThermo,
constantHeatCapacity,
soaveRedlichKwong
);
makeBasicRealFluidMixture
(
pureMixture,
constTransport,
realGasSpecieThermo,
constantHeatCapacity,
redlichKwong
);
makeBasicRealFluidMixture
(
pureMixture,
constTransport,
realGasSpecieThermo,
constantHeatCapacity,
pengRobinson
);
makeBasicRealFluidMixture
(
pureMixture,
constTransport,
realGasSpecieThermo,
constantHeatCapacity,
aungierRedlichKwong
);
makeBasicRealFluidMixture
(
pureMixture,
constTransport,
realGasSpecieThermo,
constantHeatCapacity,
soaveRedlichKwong
);
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

11
src/thermophysicalModels/basic/mixtures/basicMixture/makeBasicMixture.H
View file

@ -68,6 +68,17 @@ defineTemplateTypeNameAndDebugWithName \
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#define makeBasicRealFluidMixture(Mixture,Transport,SpecieThermo,Thermo,EqnOfState) \
\
typedef Mixture<Transport<SpecieThermo<Thermo<EqnOfState> > > > \
Mixture##Transport##SpecieThermo##Thermo##EqnOfState; \
\
defineTemplateTypeNameAndDebugWithName \
(Mixture##Transport##SpecieThermo##Thermo##EqnOfState, \
#Mixture"<"#Transport"<"#SpecieThermo"<"#Thermo"<"#EqnOfState">>>>", 0)
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

7
src/thermophysicalModels/basic/psiThermo/basicPsiThermo/basicPsiThermo.C
View file

@ -1,9 +1,9 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | foam-extend: Open Source CFD
\\ / O peration | Version: 3.2
\\ / A nd | Web: http://www.foam-extend.org
\\/ M anipulation | For copyright notice see file Copyright
\\ / O peration |
\\ / A nd | For copyright notice see file Copyright
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of foam-extend.
@ -62,4 +62,5 @@ const Foam::volScalarField& Foam::basicPsiThermo::drhode() const
return const_cast<volScalarField&>(volScalarField::null());
}
// ************************************************************************* //

2
src/thermophysicalModels/basic/psiThermo/basicPsiThermo/makeBasicPsiThermo.H
View file

@ -58,8 +58,6 @@ addToRunTimeSelectionTable \
)
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#define makeBasicRealGasThermo(Cthermo,Mixture,Transport,SpecieThermo,Thermo,EqnOfState) \

6
src/thermophysicalModels/basic/psiThermo/realGasEThermo/realGasEThermo.C
View file

@ -221,6 +221,7 @@ void Foam::realGasEThermo<MixtureType>::correct()
}
}
template<class MixtureType>
Foam::tmp<Foam::scalarField> Foam::realGasEThermo<MixtureType>::e
(
@ -265,6 +266,7 @@ Foam::tmp<Foam::scalarField> Foam::realGasEThermo<MixtureType>::e
return te;
}
template<class MixtureType>
Foam::tmp<Foam::scalarField> Foam::realGasEThermo<MixtureType>::rho
(
@ -287,6 +289,7 @@ Foam::tmp<Foam::scalarField> Foam::realGasEThermo<MixtureType>::rho
return trho;
}
template<class MixtureType>
Foam::tmp<Foam::scalarField> Foam::realGasEThermo<MixtureType>::Cp
(
@ -309,6 +312,7 @@ Foam::tmp<Foam::scalarField> Foam::realGasEThermo<MixtureType>::Cp
return tCp;
}
template<class MixtureType>
Foam::tmp<Foam::volScalarField> Foam::realGasEThermo<MixtureType>::Cp() const
{
@ -376,6 +380,7 @@ Foam::tmp<Foam::scalarField> Foam::realGasEThermo<MixtureType>::Cv
return tCv;
}
// CL: Maybe this function should be changed so that it is not "const" function anymore
template<class MixtureType>
Foam::tmp<Foam::volScalarField> Foam::realGasEThermo<MixtureType>::rho() const
@ -442,6 +447,7 @@ Foam::tmp<Foam::volScalarField> Foam::realGasEThermo<MixtureType>::rho() const
return trho;
}
template<class MixtureType>
Foam::tmp<Foam::volScalarField> Foam::realGasEThermo<MixtureType>::Cv() const
{

88
src/thermophysicalModels/basic/psiThermo/realGasEThermo/realGasEThermos.C
View file

@ -28,14 +28,10 @@ Institut für Thermodynamik
Technische Universität Braunschweig
Germany
\*---------------------------------------------------------------------------*/
#include "makeBasicPsiThermo.H"
#include "redlichKwong.H"
#include "pengRobinson.H"
#include "aungierRedlichKwong.H"
@ -44,7 +40,7 @@ Germany
#include "realGasSpecieThermo.H"
#include "constTransport.H"
#include "sutherlandTransport.H"
#include "constantHeatCapacity.H"
#include "pureMixture.H"
#include "realGasEThermo.H"
@ -136,6 +132,88 @@ makeBasicRealGasThermo
soaveRedlichKwong
);
makeBasicRealGasThermo
(
realGasEThermo,
pureMixture,
sutherlandTransport,
realGasSpecieThermo,
constantHeatCapacity,
pengRobinson
);
makeBasicRealGasThermo
(
realGasEThermo,
pureMixture,
sutherlandTransport,
realGasSpecieThermo,
constantHeatCapacity,
aungierRedlichKwong
);
makeBasicRealGasThermo
(
realGasEThermo,
pureMixture,
sutherlandTransport,
realGasSpecieThermo,
constantHeatCapacity,
redlichKwong
);
makeBasicRealGasThermo
(
realGasEThermo,
pureMixture,
sutherlandTransport,
realGasSpecieThermo,
constantHeatCapacity,
soaveRedlichKwong
);
makeBasicRealGasThermo
(
realGasEThermo,
pureMixture,
constTransport,
realGasSpecieThermo,
constantHeatCapacity,
pengRobinson
);
makeBasicRealGasThermo
(
realGasEThermo,
pureMixture,
constTransport,
realGasSpecieThermo,
constantHeatCapacity,
aungierRedlichKwong
);
makeBasicRealGasThermo
(
realGasEThermo,
pureMixture,
constTransport,
realGasSpecieThermo,
constantHeatCapacity,
redlichKwong
);
makeBasicRealGasThermo
(
realGasEThermo,
pureMixture,
constTransport,
realGasSpecieThermo,
constantHeatCapacity,
soaveRedlichKwong
);
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam

7
src/thermophysicalModels/basic/psiThermo/realGasHThermo/realGasHThermo.C
View file

@ -220,6 +220,7 @@ void Foam::realGasHThermo<MixtureType>::correct()
}
}
template<class MixtureType>
Foam::tmp<Foam::scalarField> Foam::realGasHThermo<MixtureType>::h
(
@ -242,6 +243,7 @@ Foam::tmp<Foam::scalarField> Foam::realGasHThermo<MixtureType>::h
return th;
}
template<class MixtureType>
Foam::tmp<Foam::scalarField> Foam::realGasHThermo<MixtureType>::h
(
@ -264,6 +266,7 @@ Foam::tmp<Foam::scalarField> Foam::realGasHThermo<MixtureType>::h
return th;
}
template<class MixtureType>
Foam::tmp<Foam::scalarField> Foam::realGasHThermo<MixtureType>::rho
(
@ -286,6 +289,7 @@ Foam::tmp<Foam::scalarField> Foam::realGasHThermo<MixtureType>::rho
return trho;
}
template<class MixtureType>
Foam::tmp<Foam::scalarField> Foam::realGasHThermo<MixtureType>::Cp
(
@ -308,6 +312,7 @@ Foam::tmp<Foam::scalarField> Foam::realGasHThermo<MixtureType>::Cp
return tCp;
}
template<class MixtureType>
Foam::tmp<Foam::volScalarField> Foam::realGasHThermo<MixtureType>::Cp() const
{
@ -375,6 +380,7 @@ Foam::tmp<Foam::scalarField> Foam::realGasHThermo<MixtureType>::Cv
return tCv;
}
// CL: Maybe this function should be changed so that it is not "const" fucntion anymore
template<class MixtureType>
Foam::tmp<Foam::volScalarField> Foam::realGasHThermo<MixtureType>::rho() const
@ -441,6 +447,7 @@ Foam::tmp<Foam::volScalarField> Foam::realGasHThermo<MixtureType>::rho() const
return trho;
}
template<class MixtureType>
Foam::tmp<Foam::volScalarField> Foam::realGasHThermo<MixtureType>::Cv() const
{

91
src/thermophysicalModels/basic/psiThermo/realGasHThermo/realGasHThermos.C
View file

@ -28,14 +28,10 @@ Institut für Thermodynamik
Technische Universität Braunschweig
Germany
\*---------------------------------------------------------------------------*/
#include "makeBasicPsiThermo.H"
#include "redlichKwong.H"
#include "pengRobinson.H"
#include "aungierRedlichKwong.H"
@ -44,7 +40,7 @@ Germany
#include "realGasSpecieThermo.H"
#include "constTransport.H"
#include "sutherlandTransport.H"
#include "constantHeatCapacity.H"
#include "pureMixture.H"
#include "realGasHThermo.H"
@ -106,9 +102,6 @@ makeBasicRealGasThermo
pengRobinson
);
makeBasicRealGasThermo
(
realGasHThermo,
@ -139,6 +132,88 @@ makeBasicRealGasThermo
soaveRedlichKwong
);
makeBasicRealGasThermo
(
realGasHThermo,
pureMixture,
sutherlandTransport,
realGasSpecieThermo,
constantHeatCapacity,
pengRobinson
);
makeBasicRealGasThermo
(
realGasHThermo,
pureMixture,
sutherlandTransport,
realGasSpecieThermo,
constantHeatCapacity,
aungierRedlichKwong
);
makeBasicRealGasThermo
(
realGasHThermo,
pureMixture,
sutherlandTransport,
realGasSpecieThermo,
constantHeatCapacity,
redlichKwong
);
makeBasicRealGasThermo
(
realGasHThermo,
pureMixture,
sutherlandTransport,
realGasSpecieThermo,
constantHeatCapacity,
soaveRedlichKwong
);
makeBasicRealGasThermo
(
realGasHThermo,
pureMixture,
constTransport,
realGasSpecieThermo,
constantHeatCapacity,
pengRobinson
);
makeBasicRealGasThermo
(
realGasHThermo,
pureMixture,
constTransport,
realGasSpecieThermo,
constantHeatCapacity,
aungierRedlichKwong
);
makeBasicRealGasThermo
(
realGasHThermo,
pureMixture,
constTransport,
realGasSpecieThermo,
constantHeatCapacity,
redlichKwong
);
makeBasicRealGasThermo
(
realGasHThermo,
pureMixture,
constTransport,
realGasSpecieThermo,
constantHeatCapacity,
soaveRedlichKwong
);
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam

13
src/thermophysicalModels/specie/Make/files
View file

@ -8,17 +8,12 @@ $(atomicWeights)/atomicWeights.C
$(specie)/specie.C
$(speciesTable)/speciesTable.C
$(equationOfState)/perfectGas/perfectGas.C
$(equationOfState)/cubicEquationOfState/redlichKwong/redlichKwong.C
$(equationOfState)/cubicEquationOfState/aungierRedlichKwong/aungierRedlichKwong.C
$(equationOfState)/cubicEquationOfState/pengRobinson/pengRobinson.C
$(equationOfState)/cubicEquationOfState/soaveRedlichKwong/soaveRedlichKwong.C
$(reactions)/makeChemkinReactions.C
$(reactions)/makeReactionThermoReactions.C
$(reactions)/makeLangmuirHinshelwoodReactions.C
$(equationOfState)/cubicEquationOfState/aungierRedlichKwong/aungierRedlichKwong.C
$(equationOfState)/cubicEquationOfState/pengRobinson/pengRobinson.C
$(equationOfState)/cubicEquationOfState/redlichKwong/redlichKwong.C
$(equationOfState)/cubicEquationOfState/soaveRedlichKwong/soaveRedlichKwong.C
$(equationOfState)/cubicEquationOfState/mixtures/mixtureRedlichKwong/mixtureRedlichKwong.C
$(equationOfState)/cubicEquationOfState/mixtures/mixturePengRobinson/mixturePengRobinson.C
$(equationOfState)/cubicEquationOfState/mixtures/mixtureSoaveRedlichKwong/mixtureSoaveRedlichKwong.C
LIB = $(FOAM_LIBBIN)/libspecie

71
src/thermophysicalModels/specie/equationOfState/cubicEquationOfState/aungierRedlichKwong/aungierRedlichKwong.C
View file

@ -36,67 +36,72 @@ Germany
#include "aungierRedlichKwong.H"
#include "IOstreams.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
aungierRedlichKwong::aungierRedlichKwong(Istream& is)
Foam::aungierRedlichKwong::aungierRedlichKwong(Istream& is)
:
specie(is),
pcrit_(readScalar(is)),
Tcrit_(readScalar(is)),
azentricFactor_(readScalar(is)),
rhocrit_(readScalar(is)),
azentricFactor_(readScalar(is)),
a0_(0.42747*pow(this->RR(), 2)*pow(Tcrit_, 2)/pcrit_),
b_(0.08664*this->RR()*Tcrit_/pcrit_),
c_(this->RR()*Tcrit_/(pcrit_+(a0_/(this->W()/rhocrit_*(this->W()/rhocrit_ + b_)))) + b_ - this->W()/rhocrit_),
n_(0.4986 + 1.2735*azentricFactor_ + 0.4754*pow(azentricFactor_, 2)),
b2_(b_*b_),
b3_(b2_*b_),
b4_(b3_*b_),
b5_(b4_*b_),
c2_(c_*c_),
//CL: Only uses the default values
b2_(pow(b_,2)),
b3_(pow(b_,3)),
b4_(pow(b_,4)),
b5_(pow(b_,5)),
c2_(pow(c_,2)),
rhoMax_(1500),
rhoMin_(1e-3),
rhoMax_(1500),
// Starting GUESS for the density by ideal gas law
TSave(0.0),
rhostd_(this->rho(this->Pstd(),this->Tstd(),this->Pstd()/(this->Tstd()*this->R())))
rhostd_(this->rho(this->Pstd(), this->Tstd(), this->Pstd()/(this->Tstd()*this->R()))),
aSave(0.0),
daSave(0.0),
d2aSave(0.0),
TSave(0.0)
{
is.check("aungierRedlichKwong::aungierRedlichKwong(Istream& is)");
}
//CL: Constructed needed in OpenFOAM 2.x.x
//CL: Code works fine, but compiling problem in OpenFOAM 1.6.ext
//CL: because specie has no constructor using dict
/*
aungierRedlichKwong::aungierRedlichKwong(const dictionary& dict)
Foam::aungierRedlichKwong::aungierRedlichKwong(const dictionary& dict)
:
specie(dict),
pcrit_(readScalar(dict.subDict("equationOfState").lookup("pCritical"))),
Tcrit_(readScalar(dict.subDict("equationOfState").lookup("TCritical"))),
azentricFactor_(readScalar(dict.subDict("equationOfState").lookup("azentricFactor"))),
rhocrit_(readScalar(dict.subDict("equationOfState").lookup("rhoCritical"))),
azentricFactor_(readScalar(dict.subDict("equationOfState").lookup("azentricFactor"))),
a0_(0.42747*pow(this->RR(), 2)*pow(Tcrit_, 2)/pcrit_),
b_(0.08664*this->RR()*Tcrit_/pcrit_),
c_(this->RR()*Tcrit_/(pcrit_+(a0_/(this->W()/rhocrit_*(this->W()/rhocrit_+b_))))+b_-this->W()/rhocrit_),
n_(0.4986+1.2735*azentricFactor_+0.4754*pow(azentricFactor_,2)),
b2_(b_*b_),
b3_(b2_*b_),
b4_(b3_*b_),
b5_(b4_*b_),
c2_(pow(c_,2)),
//CL: rhoMin and rhoMax are only used as boundaries for the bisection method (see rho function)
//CL: important: rhoMin and rhoMax are not used as boundary for the newton solver
//CL: therefore, rho can be larger than rhoMax and smaller than rhoMin
rhoMin_(dict.subDict("equationOfState").lookupOrDefault("rhoMin",1e-3)),
rhoMax_(dict.subDict("equationOfState").lookupOrDefault("rhoMax",1500)),
a0_(0.42747*pow(this->RR(),2)*pow(Tcrit_,2)/pcrit_),
b_(0.08664*this->RR()*Tcrit_/pcrit_),
c_(this->RR()*Tcrit_/(pcrit_+(a0_/(this->W()/rhocrit_*(this->W()/rhocrit_+b_))))+b_-this->W()/rhocrit_),
n_(0.4986+1.2735*azentricFactor_+0.4754*pow(azentricFactor_,2)),
TSave(0.0),
b2_(pow(b_,2)),
b3_(pow(b_,3)),
b4_(pow(b_,4)),
b5_(pow(b_,5)),
c2_(pow(c_,2)),
// Starting GUESS for the density by ideal gas law
rhostd_(this->rho(Pstd,Tstd,Pstd/(Tstd*this->R())))
{}
rhostd_(this->rho(this->Pstd(), this->Tstd(), this->Pstd()/(this->Tstd()*this->R()))),
aSave(0.0),
daSave(0.0),
d2aSave(0.0),
TSave(0.0)
{
is.check("aungierRedlichKwong::aungierRedlichKwong(Istream& is)");
}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
@ -115,9 +120,11 @@ void Foam::aungierRedlichKwong::write(Ostream& os) const
os << indent << dict.dictName() << dict;
}
*/
// * * * * * * * * * * * * * * * Ostream Operator * * * * * * * * * * * * * //
Ostream& operator<<(Ostream& os, const aungierRedlichKwong& ark)
Foam::Ostream& Foam::operator<<(Ostream& os, const aungierRedlichKwong& ark)
{
os << static_cast<const specie&>(ark)<< token::SPACE
<< ark.pcrit_ << tab<< ark.Tcrit_<< tab<<ark.azentricFactor_<< tab<<ark.rhocrit_;
@ -127,8 +134,4 @@ Ostream& operator<<(Ostream& os, const aungierRedlichKwong& ark)
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// ************************************************************************* //

87
src/thermophysicalModels/specie/equationOfState/cubicEquationOfState/aungierRedlichKwong/aungierRedlichKwong.H
View file

@ -63,33 +63,33 @@ namespace Foam
class aungierRedlichKwong
:
public specie
{
protected:
// Private data
// private data
//CL: data at critical point
scalar pcrit_;
scalar Tcrit_;
scalar azentricFactor_;
scalar rhocrit_;
scalar azentricFactor_;
//Aungier Redlich Kwong factors
mutable scalar a0_;
mutable scalar b_;
mutable scalar c_;
mutable scalar n_;
scalar a0_;
scalar b_;
scalar c_;
scalar n_;
//CL: pow of constants (b_, c_) used in the code e.g. b2_=b*b;
mutable scalar b2_;
mutable scalar b3_;
mutable scalar b4_;
mutable scalar b5_;
mutable scalar c2_;
scalar b2_;
scalar b3_;
scalar b4_;
scalar b5_;
scalar c2_;
//CL: rhoMin and rhoMax are only used as boundaries for the bisection methode (see rho function)
scalar rhoMax_;
scalar rhoMin_;
scalar rhoMax_;
//- Density @STD, initialise after a0, b, c, n!
scalar rhostd_;
//CL: Variables to save the values of a, dadT and d2adT2 of the mixture
//CL: Variables must corrected for changing temperatures
@ -100,16 +100,12 @@ protected:
//CL: save the temperature for which the save coefficients (amix,dadTmix,d2adT2mix) are correct
mutable scalar TSave;
//Density @STD, initialise after a, b!
scalar rhostd_;
//Protected functions
//CL: function updates the coefficients (aSave, daSave, d2aSave)
inline void updateModelCoefficients(const scalar T) const;
public:
// Constructors
//- Construct from components
@ -135,17 +131,6 @@ public:
// Member functions
inline scalar rhostd() const;
//CL: Model coefficient a(T)
inline scalar a(const scalar T)const;
//CL: temperature deriviative of model coefficient a(T)
inline scalar dadT(const scalar T)const;
//CL: second order temperature deriviative of model coefficient a(T)
inline scalar d2adT2(const scalar T)const;
//Return Aungier Redlich Kwong factors
inline scalar a0() const;
@ -155,17 +140,28 @@ public:
inline scalar n() const;
//CL: return power of constants (b_, c_)
inline scalar b2()const;
inline scalar rhostd() const;
inline scalar b3()const;
inline scalar rhoMin() const;
inline scalar b4()const;
inline scalar rhoMax() const;
inline scalar b5()const;
inline scalar Tcrit() const;
inline scalar c2()const;
inline scalar pcrit() const;
inline scalar rhocrit() const;
inline scalar azentricFactor() const;
//CL: Model coefficient a(T)
inline scalar a(const scalar T) const;
//CL: temperature deriviative of model coefficient a(T)
inline scalar dadT(const scalar T) const;
//CL: second order temperature deriviative of model coefficient a(T)
inline scalar d2adT2(const scalar T) const;
//CL: Equation of state
inline scalar p(const scalar rho, const scalar T) const;
@ -197,17 +193,25 @@ public:
//CL: second order derivatives, not Used At The Moment
inline scalar d2pdv2(const scalar rho, const scalar T) const;
inline scalar d2pdT2(const scalar rho, scalar T) const;
inline scalar d2pdT2(const scalar rho, const scalar T) const;
inline scalar d2pdvdT(const scalar rho, const scalar T) const;
inline scalar d2vdT2(const scalar rho, const scalar T) const;
//CL: Used for internal Energy
inline scalar integral_p_dv(const scalar rho, const scalar T) const;
inline scalar integral_p_dv
(
const scalar rho,
const scalar T
) const;
//CL: Used for Entropy
inline scalar integral_dpdT_dv(const scalar rho, const scalar T) const;
inline scalar integral_dpdT_dv
(
const scalar rho,
const scalar T
) const;
//- Return density [kg/m^3]
// rho0 is the starting point of the newton solver used to calculate rho
@ -242,6 +246,8 @@ public:
inline void operator+=(const aungierRedlichKwong&);
inline void operator*=(const scalar);
// Friend operators
@ -257,6 +263,7 @@ public:
const aungierRedlichKwong&
);
// Ostream Operator
friend Ostream& operator<<(Ostream&, const aungierRedlichKwong&);

218
src/thermophysicalModels/specie/equationOfState/cubicEquationOfState/aungierRedlichKwong/aungierRedlichKwongI.H
View file

@ -45,7 +45,8 @@ inline aungierRedlichKwong::aungierRedlichKwong
const specie& sp
)
:
specie(sp)
specie(sp),
TSave(0)
{}
@ -55,8 +56,8 @@ inline aungierRedlichKwong::aungierRedlichKwong(const word& name, const aungierR
specie(name, pg),
pcrit_(pg.pcrit_),
Tcrit_(pg.Tcrit_),
azentricFactor_(pg.azentricFactor_),
rhocrit_(pg.rhocrit_),
azentricFactor_(pg.azentricFactor_),
a0_(pg.a0_),
b_(pg.b_),
c_(pg.c_),
@ -66,11 +67,16 @@ inline aungierRedlichKwong::aungierRedlichKwong(const word& name, const aungierR
b4_(pg.b4_),
b5_(pg.b5_),
c2_(pg.c2_),
rhostd_(pg.rhostd_)
rhoMin_(pg.rhoMin_),
rhoMax_(pg.rhoMax_),
rhostd_(pg.rhostd_),
aSave(0.0),
daSave(0.0),
d2aSave(0.0),
TSave(0.0)
{}
// Construct and return a clone
inline autoPtr<aungierRedlichKwong> aungierRedlichKwong::clone() const
{
@ -85,18 +91,55 @@ inline autoPtr<aungierRedlichKwong> aungierRedlichKwong::New(Istream& is)
}
// * * * * * * * * * * * * * Member Functions * * * * * * * * * * * //
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
inline scalar aungierRedlichKwong::rhostd() const
{
return rhostd_;
}
inline scalar aungierRedlichKwong::rhoMin() const
{
return rhoMin_;
}
inline scalar aungierRedlichKwong::rhoMax() const
{
return rhoMax_;
}
inline scalar aungierRedlichKwong::Tcrit() const
{
return Tcrit_;
}
inline scalar aungierRedlichKwong::pcrit() const
{
return pcrit_;
}
inline scalar aungierRedlichKwong::rhocrit() const
{
return rhocrit_;
}
inline scalar aungierRedlichKwong::azentricFactor() const
{
return azentricFactor_;
}
inline void aungierRedlichKwong::updateModelCoefficients(const scalar T)const
{
aSave=a0_*pow(T/Tcrit_,-n());
daSave=-a0_*n()*pow(T/Tcrit_,-n())/T;
d2aSave=a0_*(n() *n()+n())/(T*T)*pow(T/Tcrit_,-n());
aSave=a0_*pow(T/Tcrit_, -n_);
daSave=-a0_*n_*pow(T/Tcrit_, -n_)/T;
d2aSave=a0_*(n_*n_ + n_)/(T*T)*pow(T/Tcrit_, -n_);
//CL: saving the temperature at which the coefficients are valid
TSave=T;
@ -153,6 +196,7 @@ inline scalar aungierRedlichKwong::d2adT2(const scalar T)const
}
}
//Aungier Redlich Kwong factors
inline scalar aungierRedlichKwong::a0() const
{
@ -165,6 +209,7 @@ inline scalar aungierRedlichKwong::b()const
return b_;
}
inline scalar aungierRedlichKwong::c() const
{
return c_;
@ -175,39 +220,12 @@ inline scalar aungierRedlichKwong::n()const
return n_;
}
//CL: pow of constants (b(), c()) used in the code e.g. b2_=b*b;
inline scalar aungierRedlichKwong::b2()const
{
return b2_;
}
inline scalar aungierRedlichKwong::b3()const
{
return b3_;
}
inline scalar aungierRedlichKwong::b4()const
{
return b4_;
}
inline scalar aungierRedlichKwong::b5()const
{
return b5_;
}
inline scalar aungierRedlichKwong::c2()const
{
return c2_;
}
//returns the pressure for a given density and temperature
inline scalar aungierRedlichKwong::p(const scalar rho, const scalar T) const
{
scalar Vm = this->W()/rho;
return this->RR()*T/(Vm-b()+c())
-a(T)/(Vm*(Vm+b()));
return this->RR()*T/(Vm - b_ + c_) - a(T)/(Vm*(Vm + b_));
}
@ -218,9 +236,12 @@ inline scalar aungierRedlichKwong::dpdv(const scalar rho, const scalar T) const
scalar Vm = this->W()/rho;
scalar Vm2 = Vm*Vm;
return (a(T)*(b3()-2*b2()*c()+b()*(c()+Vm)*(c()-3*Vm)+2*Vm*pow(c()+Vm,2))
-this->RR()*T*Vm2*(b2()+2*b()*Vm+Vm2))
/(Vm2*pow(b()-c()-Vm,2)*pow(b()+Vm,2));
return
(
a(T)*(b3_ - 2*b2_*c_ + b_*(c_ + Vm)*(c_ - 3*Vm) + 2*Vm*pow(c_ + Vm, 2))
- this->RR()*T*Vm2*(b2_ + 2*b_*Vm + Vm2)
)
/(Vm2*pow(b_ - c_ - Vm, 2)*pow(b_ + Vm, 2));
}
@ -229,7 +250,7 @@ inline scalar aungierRedlichKwong::dpdv(const scalar rho, const scalar T) const
inline scalar aungierRedlichKwong::dpdT(const scalar rho, const scalar T) const
{
scalar Vm = this->W()/rho;
return this->RR()/(Vm-b()+c())-dadT(T)/(Vm*(Vm+b()));
return this->RR()/(Vm - b_ + c_)-dadT(T)/(Vm*(Vm + b_));
}
@ -238,7 +259,7 @@ inline scalar aungierRedlichKwong::dpdT(const scalar rho, const scalar T) const
//(molar values)
inline scalar aungierRedlichKwong::dvdT(const scalar rho, const scalar T) const
{
return (-1)*this->dpdT(rho,T)/this->dpdv(rho,T);
return -this->dpdT(rho, T)/this->dpdv(rho, T);
}
@ -259,7 +280,8 @@ inline scalar aungierRedlichKwong::integral_p_dv
) const
{
scalar Vm = this->W()/rho;
return this->RR()*T*log(-b()+c()+Vm)+a(T)*log(b()+Vm)/b()-a(T)*log(Vm)/b();
return this->RR()*T*log(-b_ + c_ + Vm) + a(T)/b_*(log(b_ + Vm) - log(Vm));
//return this->RR()*T*log(-b_ + c_ + Vm) + a(T)*log(b_ + Vm)/b_ - a(T)*log(Vm)/b_;
}
@ -272,7 +294,8 @@ inline scalar aungierRedlichKwong::integral_dpdT_dv
) const
{
scalar Vm = this->W()/rho;
return this->RR()*log(-b()+c()+Vm)+dadT(T)*log(b()+Vm)/b()-dadT(T)*log(Vm)/b();
return this->RR()*log(-b_ + c_ + Vm) + dadT(T)/b_*(log(b_ + Vm) - log(Vm));
//return this->RR()*log(-b_ + c_ + Vm) + dadT(T)*log(b_ + Vm)/b_ - dadT(T)*log(Vm)/b_;
}
@ -280,7 +303,7 @@ inline scalar aungierRedlichKwong::integral_dpdT_dv
inline scalar aungierRedlichKwong::d2pdT2(const scalar rho,const scalar T) const
{
scalar Vm = this->W()/rho;
return -d2adT2(T)/(Vm*(Vm+b()));
return -d2adT2(T)/(Vm*(Vm + b_));
}
@ -289,39 +312,61 @@ inline scalar aungierRedlichKwong::d2pdv2(const scalar rho,const scalar T) const
{
scalar Vm = this->W()/rho;
scalar Vm2 = Vm*Vm;
scalar Vm3 = Vm*Vm*Vm;
scalar Vm3 = Vm2*Vm;
return -2*(a(T)*(b5()-3*b4()*c()+3*b3()*(c2()-c()*Vm-Vm2)
-b2()*(c()+Vm)*(c2()-7*c()*Vm+Vm2)+3*b()*Vm*pow(c()+Vm,2)*(2*Vm-c())
-3*Vm2*pow(c()+Vm,3))+this->RR()*T*Vm3*(b3()+3*b2()*Vm
+3*b()*Vm2+Vm3))
/(Vm3*pow(b()-c()-Vm,3)*pow(b()+Vm,3));
return -2*
(
a(T)*
(
b5_ - 3*b4_*c_ + 3*b3_*(c2_ - c_*Vm - Vm2)
- b2_*(c_ + Vm)*(c2_ - 7*c_*Vm+Vm2)
+ 3*b_*Vm*pow(c_ + Vm, 2)*(2*Vm - c_)
- 3*Vm2*pow(c_ + Vm, 3)
)
+ this->RR()*T*Vm3*(b3_ + 3*b2_*Vm + 3*b_*Vm2 + Vm3)
)
/(Vm3*pow(b_ - c_ - Vm, 3)*pow(b_ + Vm, 3));
}
//(molar values)
//using second order implicit differentiation
inline scalar aungierRedlichKwong::d2vdT2(const scalar rho, const scalar T) const
inline scalar aungierRedlichKwong::d2vdT2
(
const scalar rho,
const scalar T
) const
{
scalar dpdT2=this->dpdT(rho, T)*this->dpdT(rho, T);
scalar dpdv2=this->dpdv(rho, T)*this->dpdv(rho, T);
scalar dpdv3=dpdv2*this->dpdv(rho, T);
return
-(
pow(this->dpdT(rho,T),2)*this->d2pdv2(rho,T)
+ pow(this->dpdv(rho,T),2)*this->d2pdT2(rho,T)
dpdT2*this->d2pdv2(rho, T)
+ dpdv2*this->d2pdT2(rho, T)
- 2*this->dpdv(rho, T)*this->dpdT(rho, T)*this->d2pdvdT(rho, T)
)
/(pow(this->dpdv(rho,T),3));
/dpdv3;
}
//(molar values)
inline scalar aungierRedlichKwong::d2pdvdT(const scalar rho, const scalar T) const
inline scalar aungierRedlichKwong::d2pdvdT
(
const scalar rho,
const scalar T
) const
{
scalar Vm = this->W()/rho;
scalar Vm2 = Vm*Vm;
return (dadT(T)*(b3()-2*b2()*c()+b()*(c()+Vm)*(c()-3*Vm)+2*Vm*pow(c()+Vm,2))
-this->RR()*Vm2*(b2()+2*b()*Vm+Vm2))
/(Vm2*pow(b()-c()-Vm,2)*pow(b()+Vm,2));
return
(
dadT(T)*(b3_ - 2*b2_*c_ + b_*(c_ + Vm)*(c_ - 3*Vm)+2*Vm*pow(c_ + Vm, 2))
- this->RR()*Vm2*(b2_ + 2*b_*Vm + Vm2)
)
/(Vm2*pow(b_ - c_ - Vm, 2)*pow(b_ + Vm, 2));
}
@ -334,13 +379,17 @@ inline scalar aungierRedlichKwong::integral_d2pdT2_dv
) const
{
scalar Vm = this->W()/rho;
return d2adT2(T)*log(b()+Vm)/b()-d2adT2(T)*log(Vm)/b();
return d2adT2(T)*log(b_ + Vm)/b_ - d2adT2(T)*log(Vm)/b_;
}
//Isobar expansion Coefficent beta = 1/v (dv/dt) at constant p
//(molar values)
inline scalar aungierRedlichKwong::isobarExpCoef(const scalar rho,const scalar T) const
inline scalar aungierRedlichKwong::isobarExpCoef
(
const scalar rho,
const scalar T
) const
{
return this->dvdT(rho, T)*rho/this->W();
}
@ -348,7 +397,11 @@ inline scalar aungierRedlichKwong::isobarExpCoef(const scalar rho,const scalar T
//isothemal compressiblity kappa (not Thermal conductivity)
//(molar values)
inline scalar aungierRedlichKwong::isothermalCompressiblity(const scalar rho,const scalar T) const
inline scalar aungierRedlichKwong::isothermalCompressiblity
(
const scalar rho,
const scalar T
) const
{
return this->isobarExpCoef(rho, T)/this->dpdT(rho, T);
//also possible : return -this->dvdp(rho,T)*rho/this->W();
@ -380,6 +433,7 @@ inline scalar aungierRedlichKwong::rho
label i=0;
do
{
//CL: modified Newton solver
molarVolume=molarVolumePrevIteration
-(
(this->p((this->W()/molarVolumePrevIteration),T) - p)
@ -390,8 +444,7 @@ inline scalar aungierRedlichKwong::rho
if (i>8)
{
//CL: using bisection methode as backup,
//CL: solution must be between rho=0.001 to rho=1500;
//CL: if not, change rhoMax_ and rhoMin_
//CL: solution must be between rhoMin_ to rhoMax
for(i=0; i<200; i++)
{
scalar f1 = this->p(rho1, T) - p;
@ -435,7 +488,8 @@ inline scalar aungierRedlichKwong::rho
{
FatalErrorIn
(
"inline scalar aungierRedlichKwong::rho(const scalar p, const scalar T, const scalar rho0) const "
"inline scalar aungierRedlichKwong::rho"
"(const scalar p, const scalar T, const scalar rho0) const "
) << "Maximum number of iterations exceeded"
<< abort(FatalError);
}
@ -449,7 +503,7 @@ inline scalar aungierRedlichKwong::rho
//- Return density [kg/m^3]
inline scalar aungierRedlichKwong::rho(const scalar p, const scalar T) const
{
//using perfect gas equation as starting point
//CL: using perfect gas equation as starting point
return rho(p, T, p/(this->R()*T));
}
@ -457,12 +511,17 @@ inline scalar aungierRedlichKwong::rho(const scalar p,const scalar T) const
//- Return compressibility drho/dp at T=constant [s^2/m^2]
inline scalar aungierRedlichKwong::psi(const scalar rho, const scalar T) const
{
return -this->dvdp(rho,T)*pow(rho,2)/this->W();
return -this->dvdp(rho, T)*rho*rho/this->W();
}
//- Return compression factor []
inline scalar aungierRedlichKwong::Z( const scalar p, const scalar T,const scalar rho0) const
inline scalar aungierRedlichKwong::Z
(
const scalar p,
const scalar T,
const scalar rho0
) const
{
return p/(this->R()*T*this->rho(p, T, rho0));
}
@ -470,10 +529,15 @@ inline scalar aungierRedlichKwong::Z( const scalar p, const scalar T,const scala
// * * * * * * * * * * * * * * * Member Operators * * * * * * * * * * * * * //
inline void aungierRedlichKwong::operator+=(const aungierRedlichKwong& pr)
inline void aungierRedlichKwong::operator+=(const aungierRedlichKwong& ark)
{
specie::operator+=(pr);
specie::operator+=(ark);
}
inline void aungierRedlichKwong::operator*=(const scalar s)
{
specie::operator*=(s);
}
@ -481,26 +545,28 @@ inline void aungierRedlichKwong::operator+=(const aungierRedlichKwong& pr)
inline aungierRedlichKwong operator+
(
const aungierRedlichKwong& pr1,
const aungierRedlichKwong& pr2
const aungierRedlichKwong& ark1,
const aungierRedlichKwong& ark2
)
{
return aungierRedlichKwong
(
static_cast<const specie&>(pr1)
+ static_cast<const specie&>(pr2)
static_cast<const specie&>(ark1)
+ static_cast<const specie&>(ark2)
);
}
inline aungierRedlichKwong operator*
(
const scalar s,
const aungierRedlichKwong& pr
const aungierRedlichKwong& ark
)
{
return aungierRedlichKwong(s*static_cast<const specie&>(pr));
return aungierRedlichKwong(s*static_cast<const specie&>(ark));
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam

11
src/thermophysicalModels/specie/equationOfState/cubicEquationOfState/mixtures/README
View file

@ -1,11 +0,0 @@
CL: Real gas mixture classes using mixture models.
CL: other models are possible
CL: Mixtures based on a pseudo critical point approach
--> mixtureAungierRedlichKwong
CL: Mixtures using the van der waals mixing rule
--> mixtureSoaveRedlichKwong
--> mixturePengRobinson
--> mixtureRedlichKwong

87
src/thermophysicalModels/specie/equationOfState/cubicEquationOfState/mixtures/mixtureAungierRedlichKwong/mixtureAungierRedlichKwong.C
View file

@ -1,87 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright held by original author
\\/ 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 2 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, write to the Free Software Foundation,
Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
Description
Mixture Aungier Redlich Kwong equation of state.
Author
Christian Lucas
Institut für Thermodynamik
Technische Universität Braunschweig
Germany
\*---------------------------------------------------------------------------*/
#include "mixtureAungierRedlichKwong.H"
#include "IOstreams.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
mixtureAungierRedlichKwong::mixtureAungierRedlichKwong(Istream& is)
:
aungierRedlichKwong(is),
Vcrit_(this->W()/rhocrit_),
Zcrit_(pcrit_/((this->R()*Tcrit_*rhocrit_)))
{
is.check("mixtureAungierRedlichKwong::mixtureAungierRedlichKwong(Istream& is)");
}
//CL: Constructed needed in OpenFOAM 2.x.x
//CL: Code works fine, but compiling problem in OpenFOAM 1.6.ext
//CL: because specie has no constructor using dict
/*
mixtureAungierRedlichKwong::mixtureAungierRedlichKwong(const dictionary& dict)
:
aungierRedlichKwong(dict),
Vcrit_(this->W()/rhocrit_),
Zcrit_(pcrit_/((this->R()*Tcrit_*rhocrit_)))
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
void Foam::mixtureAungierRedlichKwong::write(Ostream& os) const
{
aungierRedlichKwong::write(os);
}
*/
// * * * * * * * * * * * * * * * Ostream Operator * * * * * * * * * * * * * //
Ostream& operator<<(Ostream& os, const mixtureAungierRedlichKwong& ark)
{
os << static_cast<const specie&>(ark)<< token::SPACE;
os.check("Ostream& operator<<(Ostream& os, const mixtureAungierRedlichKwong& st)");
return os;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// ************************************************************************* //

164
src/thermophysicalModels/specie/equationOfState/cubicEquationOfState/mixtures/mixtureAungierRedlichKwong/mixtureAungierRedlichKwong.H
View file

@ -1,164 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright held by original author
\\/ 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 2 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, write to the Free Software Foundation,
Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
Class
Foam::mixtureAungierRedlichKwong
Description
Mixture Aungier Redlich Kwong equation of state.
Mixture based on a pseudo critical point approach.
For further information, see:
BOOK Title: The Properties of Gases And Liquids, Fifth Edition, McGraw-Hill, Chapter 5
Authors: B. Poling; J. Prausnitz; J. O'Connell
SourceFiles
mixtureAungierRedlichKwongI.H
mixtureAungierRedlichKwong.C
Author
Christian Lucas
Institut für Thermodynamik
Technische Universität Braunschweig
Germany
\*---------------------------------------------------------------------------*/
#ifndef mixtureAungierRedlichKwong_H
#define mixtureAungierRedlichKwong_H
#include "specie.H"
#include "aungierRedlichKwong.H"
#include "autoPtr.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
/*---------------------------------------------------------------------------*\
Class mixtureAungierRedlichKwong Declaration
\*---------------------------------------------------------------------------*/
class mixtureAungierRedlichKwong
:
public aungierRedlichKwong
{
protected:
// data at the critical point
scalar Zcrit_;
//CL. molar volume at the critical point
scalar Vcrit_;
//CL: save the concentrations of each component of the mixture
//CL: needs to be multiplied by this->W() to get the molar fractions
//mutable DynamicList<scalar> weigths;
//CL: saves a pointer to the pure component classes of the mixture
// mutable DynamicList<mixtureAungierRedlichKwong*> mixtureComponents;
public:
// Constructors
//- Construct from components
inline mixtureAungierRedlichKwong
(
const aungierRedlichKwong& sp,
scalar Tcrit,
scalar azentricFactor,
scalar Vcrit,
scalar Zcrit,
scalar rhoMin,
scalar rhoMax
);
//- Construct from Istream
mixtureAungierRedlichKwong(Istream&);
//- Construct from dictionary
//mixtureAungierRedlichKwong(const dictionary& dict);
//- Construct as named copy
inline mixtureAungierRedlichKwong(const word& name, const mixtureAungierRedlichKwong&);
//- Construct and return a clone
inline autoPtr<mixtureAungierRedlichKwong> clone() const;
// Selector from Istream
inline static autoPtr<mixtureAungierRedlichKwong> New(Istream& is);
// Member functions
// I-O
//- Write to Ostream
//void write(Ostream& os) const;
// Member operators
inline void operator+=(const mixtureAungierRedlichKwong&);
// Friend operators
inline friend mixtureAungierRedlichKwong operator+
(
const mixtureAungierRedlichKwong&,
const mixtureAungierRedlichKwong&
);
inline friend mixtureAungierRedlichKwong operator*
(
const scalar s,
const mixtureAungierRedlichKwong&
);
// Ostream Operator
friend Ostream& operator<<(Ostream&, const mixtureAungierRedlichKwong&);
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#include "mixtureAungierRedlichKwongI.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

212
src/thermophysicalModels/specie/equationOfState/cubicEquationOfState/mixtures/mixtureAungierRedlichKwong/mixtureAungierRedlichKwongI.H
View file

@ -1,212 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright held by original author
\\/ 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 2 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, write to the Free Software Foundation,
Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
Author
Christian Lucas
Institut für Thermodynamik
Technische Universität Braunschweig
Germany
\*---------------------------------------------------------------------------*/
#include "mixtureAungierRedlichKwong.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
// Construct from components
inline mixtureAungierRedlichKwong::mixtureAungierRedlichKwong
(
const aungierRedlichKwong& sp,
scalar Tcrit,
scalar azentricFactor,
scalar Vcrit,
scalar Zcrit,
scalar rhoMin,
scalar rhoMax
)
:
aungierRedlichKwong(sp),
Vcrit_(Vcrit),
Zcrit_(Zcrit)
{
//CL: Saving critical data
Tcrit_=Tcrit;
pcrit_=Zcrit*this->RR*Tcrit/Vcrit;
rhocrit_=this->W()/Vcrit_;
azentricFactor_=azentricFactor;
//CL: calculating the aungier redlich kwong coefficience
a0_=0.42747*pow(this->RR,2)*pow(Tcrit_,2)/pcrit_;
b_=0.08664*this->RR*Tcrit_/pcrit_;
c_=this->RR*Tcrit_/(pcrit_+(a0_/(this->W()/rhocrit_*(this->W()/rhocrit_+b_))))+b_-this->W()/rhocrit_;
n_=0.4986+1.2735*azentricFactor_+0.4754*pow(azentricFactor_,2);
//CL: set rhoMax and rhoMin
rhoMin_=rhoMin;
rhoMax_=rhoMax;
//CL:
b2_=pow(b_,2);
b3_=pow(b_,3);
b4_=pow(b_,4);
b5_=pow(b_,5);
c2_=pow(c_,2);
//CL: Starting GUESS for the density by ideal gas law
rhostd_=this->rho(Pstd,Tstd,Pstd/(Tstd*this->R()));
}
// Construct as named copy
inline mixtureAungierRedlichKwong::mixtureAungierRedlichKwong(const word& name, const mixtureAungierRedlichKwong& pg)
:
aungierRedlichKwong(name, pg),
Vcrit_(pg.Vcrit_),
Zcrit_(pg.Zcrit_)
{}
// Construct and return a clone
inline autoPtr<mixtureAungierRedlichKwong> mixtureAungierRedlichKwong::clone() const
{
return autoPtr<mixtureAungierRedlichKwong>(new mixtureAungierRedlichKwong(*this));
}
// Selector from Istream
inline autoPtr<mixtureAungierRedlichKwong> mixtureAungierRedlichKwong::New(Istream& is)
{
return autoPtr<mixtureAungierRedlichKwong>(new mixtureAungierRedlichKwong(is));
}
// * * * * * * * * * * * * * * * Member Operators * * * * * * * * * * * * * //
inline void mixtureAungierRedlichKwong::operator+=(const mixtureAungierRedlichKwong& pr)
{
scalar molr1 = this->nMoles();
aungierRedlichKwong::operator+=(pr);
molr1 /= this->nMoles();
scalar molr2 = pr.nMoles()/this->nMoles();
//CL: calculating new critical point data
Tcrit_=molr1*Tcrit_+molr2*pr.Tcrit_;
Zcrit_=molr1*Zcrit_+molr2*pr.Zcrit_;
Vcrit_=molr1*Vcrit_+molr2*pr.Vcrit_;
pcrit_=Zcrit_*this->RR*Tcrit_/Vcrit_;
rhocrit_=this->W()/Vcrit_;
//CL: calculating new azentric factor
azentricFactor_=molr1*azentricFactor_+molr2*pr.azentricFactor_;
//CL: set rhoMax and rhoMin
rhoMin_=min(rhoMin_, pr.rhoMin_);
rhoMax_=max(rhoMax_, pr.rhoMax_);
//CL: calculating the aungier redlich kwong coefficience
a0_=0.42747*pow(this->RR,2)*pow(Tcrit_,2)/pcrit_;
b_=0.08664*this->RR*Tcrit_/pcrit_;
c_=this->RR*Tcrit_/(pcrit_+(a0_/(this->W()/rhocrit_*(this->W()/rhocrit_+b_))))+b_-this->W()/rhocrit_;
n_=0.4986+1.2735*azentricFactor_+0.4754*pow(azentricFactor_,2);
//CL:
b2_=pow(b_,2);
b3_=pow(b_,3);
b4_=pow(b_,4);
b5_=pow(b_,5);
c2_=pow(c_,2);
//CL: Starting GUESS for the density by ideal gas law
rhostd_=this->rho(Pstd,Tstd,Pstd/(Tstd*this->R()));
}
// * * * * * * * * * * * * * * * Friend Operators * * * * * * * * * * * * * //
inline mixtureAungierRedlichKwong operator+
(
const mixtureAungierRedlichKwong& pr1,
const mixtureAungierRedlichKwong& pr2
)
{
aungierRedlichKwong ark
(
static_cast<const aungierRedlichKwong&>(pr1)
+static_cast<const aungierRedlichKwong&>(pr2)
);
return mixtureAungierRedlichKwong
(
ark,
pr1.nMoles()/ark.nMoles()*pr1.Tcrit_
+pr2.nMoles()/ark.nMoles()*pr2.Tcrit_,
pr1.nMoles()/ark.nMoles()*pr1.azentricFactor_
+pr2.nMoles()/ark.nMoles()*pr2.azentricFactor_,
pr1.nMoles()/ark.nMoles()*pr1.Vcrit_
+pr2.nMoles()/ark.nMoles()*pr2.Vcrit_,
pr1.nMoles()/ark.nMoles()*pr1.Zcrit_
+pr2.nMoles()/ark.nMoles()*pr2.Zcrit_,
min(pr1.rhoMin_,pr2.rhoMin_),
max(pr1.rhoMax_,pr2.rhoMax_)
);
}
inline mixtureAungierRedlichKwong operator*
(
const scalar s,
const mixtureAungierRedlichKwong& pr
)
{
return mixtureAungierRedlichKwong
(
s*static_cast<const aungierRedlichKwong&>(pr),
pr.Tcrit_,
pr.azentricFactor_,
pr.Vcrit_,
pr.Zcrit_,
pr.rhoMin_,
pr.rhoMax_
);
}
} // End namespace Foam
// ************************************************************************* //

157
src/thermophysicalModels/specie/equationOfState/cubicEquationOfState/mixtures/mixturePengRobinson/mixturePengRobinson.C
View file

@ -1,157 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright held by original author
\\/ 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 2 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, write to the Free Software Foundation,
Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
Description
Peng Robionson equation of state for mixtures.
Author
Christian Lucas
Institut für Thermodynamik
Technische Universität Braunschweig
Germany
\*---------------------------------------------------------------------------*/
#include "mixturePengRobinson.H"
#include "IOstreams.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
mixturePengRobinson::mixturePengRobinson(Istream& is)
:
pengRobinson(is),
singleComponent(1),
numOfComp(1),
//CL: no real gas mixture correction when stream constructor is used
realMixtureCorr_(false)
{
//CL: set size of weigths, mixtureComponents ... to 10,
//CL: when more mixture components are used
//CL: size of the DynamicLis increases automatically
weigths.setSize(10);
mixtureComponents.setSize(10);
//Save a pointer of this object in the mixtureComponents array
mixtureComponents[0]=this;
is.check("mixturePengRobinson::mixturePengRobinson(Istream& is)");
}
//CL: Constructed needed in OpenFOAM 2.x.x
//CL: Code works fine, but compiling problem in OpenFOAM 1.6.ext
//CL: because specie has no constructor using dict
/*
mixturePengRobinson::mixturePengRobinson(const dictionary& dict)
:
pengRobinson(dict),
numOfComp(1),
singleComponent(1),
realMixtureCorr_(dict.subDict("equationOfState").lookupOrDefault("realMixtureCorrection",false))
{
//CL: set size of weigths, mixtureComponents ... to 10,
//CL: when more mixture components are used
//CL: size of the DynamicLis increases automatically
weigths.setSize(10);
mixtureComponents.setSize(10);
//Save a pointer of this object in the mixtureComponents array
mixtureComponents[0]=this;
if(realMixtureCorr_==true)
{
//CL: reads number of mixture components
//CL: this number is needed to calculate the number of correction coefficients
nCom_=dict.subDict("equationOfState").lookupOrDefault("numberOfComponents",1);
if (nCom_<1||nCom_==1)
{
FatalErrorIn
(
"mixturePengRobinson::mixturePengRobinson(const dictionary& dict)"
) << "mixture has only one component or number is not defined correctly, "
<< "recheck numberOfComponents in the thermophysicalProperties dict of your case"
<< abort(FatalError);
}
label i;
label j;
label k=0;
// CL: saves the number of mixtureCorrectionCoefficients needed for this mixture
// CL: need to be set to 1
label sizeOfVector_=1;
//CL: size of the vector depends on the number of components
//CL: determine the size of the vector
for (i=3; i<=nCom_;i++)
{
sizeOfVector_=sizeOfVector_+(i-1);
}
//CL: setting the size
realMixtureCorrCoef_.setSize(sizeOfVector_);
//CL: Reading the real mixture correction coefficients
//CL: Naming convention e.g for 3 component mixture:
//CL: mixtureCorrectionCoefficient_12, mixtureCorrectionCoefficient_13, mixtureCorrectionCoefficient_23
for(i=1;i<=nCom_-1;i++)
{
for(j=i+1;j<=nCom_;j++)
{
realMixtureCorrCoef_[k]=dict.subDict("equationOfState")
.lookupOrDefault("mixtureCorrectionCoefficient_" + Foam::name(i) + Foam::name(j),0.0);
k++;
}
}
}
}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
void Foam::mixturePengRobinson::write(Ostream& os) const
{
pengRobinson::write(os);
}
*/
// * * * * * * * * * * * * * * * Ostream Operator * * * * * * * * * * * * * //
Ostream& operator<<(Ostream& os, const mixturePengRobinson& pr)
{
os << static_cast<const pengRobinson&>(pr)<< token::SPACE ;
os.check("Ostream& operator<<(Ostream& os, const mixturePengRobinson& st)");
return os;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// ************************************************************************* //

274
src/thermophysicalModels/specie/equationOfState/cubicEquationOfState/mixtures/mixturePengRobinson/mixturePengRobinson.H
View file

@ -1,274 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright held by original author
\\/ 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 2 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, write to the Free Software Foundation,
Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
Class
Foam::mixturePengRobinson
Description
Peng Robinson equation of state for multiple component mixtures using the van der Waals mixing rule
Mixing Rule --> see Paper
Title: Van der waals mixing rules for cublic equations of state. Applications for supercritical fluids extraction modelling
Authors: T.Y. Kwak and G.A. Mansoori
Journal: Chemical Engineering Science, Vol 41, No. 5, pp. 1303-1309, 1986
SourceFiles
mixturePengRobinsonI.H
mixturePengRobinson.C
Author
Christian Lucas
Institut für Thermodynamik
Technische Universität Braunschweig
Germany
\*---------------------------------------------------------------------------*/
#ifndef mixturePengRobinson_H
#define mixturePengRobinson_H
#include "specie.H"
#include "autoPtr.H"
#include "pengRobinson.H"
#include "scalarList.H"
#include "DynamicList.H"
#include "label.H"
#include "List.H"
#include "dictionary.H"
#include "scalar.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
/*---------------------------------------------------------------------------*\
Class mixturePengRobinson Declaration
\*---------------------------------------------------------------------------*/
class mixturePengRobinson
:
public pengRobinson
{
protected:
//CL: bool used to make sure the model behaves like a single component model when needed
//CL: this is needed during the construction of the mixture
mutable bool singleComponent;
//CL: save the concentrations of each component of the mixture
//CL: needs to be multiplied by this->W() to get the molar fractions
mutable DynamicList<scalar> weigths;
//CL: saves a pointer to the pure component classes of the mixture
mutable DynamicList<mixturePengRobinson*> mixtureComponents;
//CL: counts the number of components
mutable label numOfComp;
//Protected functions
//CL: function updates the coefficients (aSave, daSave, d2aSave, b) of the mixture
//CL: this is the function with the mixing rule
inline void updateModelCoefficients(const scalar T) const;
//CL: Variables used in real gas mixture correction
//CL: If true, the real gas mixture correction is used
mutable bool realMixtureCorr_;
//CL: number of mixture components, needed to calculate the mixture correction factors needed
//CL: do not mistake this variable with numOfComp,
//CL: numOfComp is a counter to counts the number of components while the mixture is constructed
mutable label nCom_;
//CL: stores real mixture correction coefficients
mutable DynamicList<scalar> realMixtureCorrCoef_;
public:
// Constructors
//- Construct from components
//CL: needed for operator*
inline mixturePengRobinson
(
const pengRobinson& pr,
label numOfComp,
DynamicList<scalar> weigths,
DynamicList<mixturePengRobinson*> mixtureComponents,
scalar a0,
scalar b,
scalar Tcrit,
scalar n,
scalar rhostd,
scalar rhoMin,
scalar rhoMax
);
//- Construct from components
//CL: needed for operator+
inline mixturePengRobinson
(
const pengRobinson& pr,
label numOfComp,
DynamicList<scalar> weigths,
DynamicList<mixturePengRobinson*> mixtureComponents,
scalar rhoMin,
scalar rhoMax
);
//- Construct from Istream
mixturePengRobinson(Istream&);
//- Construct from dictionary
//mixturePengRobinson(const dictionary& dict);
//- Construct as named copy
inline mixturePengRobinson(const word& name, const mixturePengRobinson&);
//- Construct and return a clone
inline autoPtr<mixturePengRobinson> clone() const;
// Selector from Istream
inline static autoPtr<mixturePengRobinson> New(Istream& is);
// Member functions
//CL: mixture coefficience
inline scalar amix(const scalar T, const label i, const label j) const;
inline scalar dadTmix(const scalar T, const label i, const label j) const;
inline scalar d2adT2mix(const scalar T, const label i, const label j) const;
inline scalar p(const scalar rho, const scalar T) const;
//first order derivatives
inline scalar dpdv(const scalar rho, const scalar T) const;
inline scalar dpdT(const scalar rho, const scalar T) const;
inline scalar dvdT(const scalar rho, const scalar T) const;
inline scalar dvdp(const scalar rho, const scalar T) const;
inline scalar isobarExpCoef(const scalar rho, const scalar T) const;
inline scalar isothermalCompressiblity
(
const scalar rho,
const scalar T
) const;
// Used for cv
inline scalar integral_d2pdT2_dv
(
const scalar rho,
const scalar T
) const;
//Used for internal Energy
inline scalar integral_p_dv(const scalar rho, const scalar T) const;
// Used for Entropy
inline scalar integral_dpdT_dv(const scalar rho, const scalar T) const;
// second order derivatives, not Used At The Moment
inline scalar d2pdv2(const scalar rho, const scalar T) const;
inline scalar d2pdT2(const scalar rho, const scalar T) const;
inline scalar d2pdvdT(const scalar rho, const scalar T) const;
inline scalar d2vdT2(const scalar rho, const scalar T) const;
//- Return density [kg/m^3]
// rho0 is the starting point of the newton solver used to calculate rho
inline scalar rho
(
const scalar p,
const scalar T,
const scalar rho0
) const;
inline scalar rho(const scalar p, const scalar T) const;
//- Return compressibility drho/dp at T=constant [s^2/m^2]
inline scalar psi(const scalar rho, const scalar T) const;
//- Return compression factor []
inline scalar Z
(
const scalar p,
const scalar T,
const scalar rho0
) const;
// I-O
//- Write to Ostream
//void write(Ostream& os) const;
// Member operators
inline void operator+=(const mixturePengRobinson&);
// Friend operators
inline friend mixturePengRobinson operator+
(
const mixturePengRobinson&,
const mixturePengRobinson&
);
inline friend mixturePengRobinson operator*
(
const scalar s,
const mixturePengRobinson&
);
// Ostream Operator
friend Ostream& operator<<(Ostream&, const mixturePengRobinson&);
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#include "mixturePengRobinsonI.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

497
src/thermophysicalModels/specie/equationOfState/cubicEquationOfState/mixtures/mixturePengRobinson/mixturePengRobinsonI.H
View file

@ -1,497 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright held by original author
\\/ 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 2 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, write to the Free Software Foundation,
Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
Author
Christian Lucas
Institut für Thermodynamik
Technische Universität Braunschweig
Germany
\*---------------------------------------------------------------------------*/
#include "mixturePengRobinson.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
// Construct from components
// CL: needed for operator*
inline mixturePengRobinson::mixturePengRobinson
(
const pengRobinson& pr,
label numOfComp,
DynamicList<scalar> weigths,
DynamicList<mixturePengRobinson*> mixtureComponents,
scalar a0,
scalar b,
scalar Tcrit,
scalar n,
scalar rhostd,
scalar rhoMin,
scalar rhoMax
)
:
pengRobinson(pr),
singleComponent(1),
weigths(weigths),
mixtureComponents(mixtureComponents),
numOfComp(numOfComp)
{
a0_=a0;
b_=b;
Tcrit_=Tcrit;
n_=n;
rhostd_=rhostd;
rhoMin_=rhoMin;
rhoMax_=rhoMax;
}
// Construct from components
// CL: needed for operator+
inline mixturePengRobinson::mixturePengRobinson
(
const pengRobinson& pr,
label numOfComp,
DynamicList<scalar> weigths,
DynamicList<mixturePengRobinson*> mixtureComponents,
scalar rhoMin,
scalar rhoMax
)
:
pengRobinson(pr),
singleComponent(0),
weigths(weigths),
mixtureComponents(mixtureComponents),
numOfComp(numOfComp)
{
TSave=0.0;
rhoMin_=rhoMin;
rhoMax_=rhoMax;
rhostd_=this->rho(this->Pstd(), this->Tstd(), this->Pstd()*this->W()/(this->Tstd()*this->R()));
}
// Construct as named copy
inline mixturePengRobinson::mixturePengRobinson(const word& name, const mixturePengRobinson& pr)
:
pengRobinson(name, pr)
{}
// Construct and return a clone
inline autoPtr<mixturePengRobinson> mixturePengRobinson::clone() const
{
return autoPtr<mixturePengRobinson>(new mixturePengRobinson(*this));
}
// Selector from Istream
inline autoPtr<mixturePengRobinson> mixturePengRobinson::New(Istream& is)
{
return autoPtr<mixturePengRobinson>(new mixturePengRobinson(is));
}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
//CL: updates the coefficents of the model after the construction of the mixture
//CL: uses the van der waals mixing rule
inline void mixturePengRobinson::updateModelCoefficients(const scalar T) const
{
// CL: to make sure that the coefficents are only updated if the mixture has more than 1 component
if (singleComponent==0)
{
// Checking if the mixture coefficient were already calculated for this temperature
if(TSave!=T)
{
label i,j,k;
aSave=0;
daSave=0;
d2aSave=0;
b_=0;
for (i=0;i<numOfComp;i++)
{
//CL: getting a(T), dadT(T) and d2adT2(T) for the mixture
//CL: using van der waals mixing rule
for (j=0;j<numOfComp;j++)
{
//CL: use mixture correction cofficients k_ij (see paper reference)
if(mixtureComponents[0]->realMixtureCorr_==true)
{
// first and second order temperature derivative of the van der waals mixing rule for a(T)
if(i==j)
{
aSave=aSave+amix(T,i,j);
daSave=daSave+dadTmix(T,i,j);
d2aSave=d2aSave+d2adT2mix(T,i,j);
}
else
{
//CL: gives the position of the correction factor in the vector realMixtureCorrCoef_
k=i+j-1;
aSave=aSave+amix(T,i,j)*(1-mixtureComponents[0]->realMixtureCorrCoef_[k]);
daSave=daSave+dadTmix(T,i,j)*(1-mixtureComponents[0]->realMixtureCorrCoef_[k]);
d2aSave=d2aSave+d2adT2mix(T,i,j)*(1-mixtureComponents[0]->realMixtureCorrCoef_[k]);
}
}
else
{
aSave=aSave+amix(T,i,j);
daSave=daSave+dadTmix(T,i,j);
d2aSave=d2aSave+d2adT2mix(T,i,j);
}
}
//CL: getting b for the mixture
//CL: using van der waals mixing rule
b_=b_+weigths[i]*mixtureComponents[i]->b()*this->W();
//CL: saving the temperature at which the mixture coefficients are valid
TSave=T;
}
b2_=b_*b_;
b3_=pow(b_,3);
b4_=pow(b_,4);
b5_=pow(b_,5);
}
}
}
//CL: returns part of "a" for the mixture
//CL: Important: a is a function of T
//CL: TODO: RECHECK
inline scalar mixturePengRobinson::amix(const scalar T, const label i, const label j) const
{
return weigths[i]*weigths[j]*this->W()*this->W()*
pow(mixtureComponents[i]->a(T)*mixtureComponents[j]->a(T),0.5);
}
//CL: returns part of the first derivative of "a" for the mixture
//CL: Important: a is a function of T
//CL: TODO: RECHECK
inline scalar mixturePengRobinson::dadTmix(const scalar T, const label i, const label j) const
{
return 0.5*weigths[i]*weigths[j]*this->W()*this->W()
*(pow(mixtureComponents[i]->a(T)/mixtureComponents[j]->a(T),0.5)*mixtureComponents[j]->dadT(T)
+pow(mixtureComponents[j]->a(T)/mixtureComponents[i]->a(T),0.5)*mixtureComponents[i]->dadT(T));
}
//CL: returns part of the second derivative of "a" for the mixture
//CL: Important: a is a function of T
//CL: TODO: RECHECK
inline scalar mixturePengRobinson::d2adT2mix(const scalar T, const label i, const label j) const
{
return 0.5*weigths[i]*weigths[j]*this->W()*this->W()*
(
pow(mixtureComponents[i]->a(T)/mixtureComponents[j]->a(T),0.5)*mixtureComponents[j]->d2adT2(T)
+pow(mixtureComponents[j]->a(T)/mixtureComponents[i]->a(T),0.5)*mixtureComponents[i]->d2adT2(T)
+pow(mixtureComponents[i]->a(T)*mixtureComponents[j]->a(T),-0.5)
*mixtureComponents[i]->dadT(T)*mixtureComponents[j]->dadT(T)
-0.5*pow(mixtureComponents[i]->a(T)/mixtureComponents[j]->a(T),0.5)*1/mixtureComponents[j]->a(T)
*pow(mixtureComponents[j]->dadT(T),2)
-0.5*pow(mixtureComponents[j]->a(T)/mixtureComponents[i]->a(T),0.5)
*1/mixtureComponents[i]->a(T)*pow(mixtureComponents[i]->dadT(T),2)
);
}
//returns the pressure for a given density and temperature
inline scalar mixturePengRobinson::p(const scalar rho, const scalar T) const
{
updateModelCoefficients(T);
return pengRobinson::p(rho,T);
}
//Real deviative dp/dv at constant temperature
//(molar values)
inline scalar mixturePengRobinson::dpdv(const scalar rho, const scalar T) const
{
updateModelCoefficients(T);
return pengRobinson::dpdv(rho,T);
}
//Real deviative dp/dT at constant molar volume
//(molar values)
inline scalar mixturePengRobinson::dpdT(const scalar rho, const scalar T) const
{
updateModelCoefficients(T);
return pengRobinson::dpdT(rho,T);
}
//Real deviative dv/dT at constant pressure
//(molar values)
inline scalar mixturePengRobinson::dvdT(const scalar rho, const scalar T) const
{
updateModelCoefficients(T);
return pengRobinson::dvdT(rho,T);
}
//Real deviative dv/dp at constant temperature
//(molar values)
inline scalar mixturePengRobinson::dvdp(const scalar rho, const scalar T) const
{
updateModelCoefficients(T);
return pengRobinson::dvdp(rho,T);
}
//needed to calculate the internal energy
//(molar values)
inline scalar mixturePengRobinson::integral_p_dv
(
const scalar rho,
const scalar T
) const
{
updateModelCoefficients(T);
return pengRobinson::integral_p_dv(rho,T);
}
//needed to calculate the entropy
//(molar values)
inline scalar mixturePengRobinson::integral_dpdT_dv
(
const scalar rho,
const scalar T
) const
{
updateModelCoefficients(T);
return pengRobinson::integral_dpdT_dv(rho,T);
}
//(molar values)
inline scalar mixturePengRobinson::d2pdT2(const scalar rho, const scalar T) const
{
updateModelCoefficients(T);
return pengRobinson::d2pdT2(rho,T);
}
//(molar values)
inline scalar mixturePengRobinson::d2pdv2(const scalar rho, const scalar T) const
{
updateModelCoefficients(T);
return pengRobinson::d2pdv2(rho,T);
}
//(molar values)
inline scalar mixturePengRobinson::d2vdT2(const scalar rho, const scalar T) const
{
updateModelCoefficients(T);
return pengRobinson::d2vdT2(rho,T);
}
//(molar values)
inline scalar mixturePengRobinson::d2pdvdT(const scalar rho, const scalar T) const
{
updateModelCoefficients(T);
return pengRobinson::d2pdvdT(rho,T);
}
// the result of this intergal is needed for the nasa based cp polynomial
//(molar values)
inline scalar mixturePengRobinson::integral_d2pdT2_dv
(
const scalar rho,
const scalar T
) const
{
updateModelCoefficients(T);
return pengRobinson::integral_d2pdT2_dv(rho,T);
}
//Isobar expansion Coefficent beta = 1/v (dv/dt) at constant p
//(molar values)
inline scalar mixturePengRobinson::isobarExpCoef
(
const scalar rho,
const scalar T
) const
{
updateModelCoefficients(T);
return pengRobinson::isobarExpCoef(rho,T);
}
//isothemal compressiblity kappa (not Thermal conductivity)
//(molar values)
inline scalar mixturePengRobinson::isothermalCompressiblity
(
const scalar rho,
const scalar T
) const
{
updateModelCoefficients(T);
return pengRobinson::isothermalCompressiblity(rho,T);
}
//- Return density [kg/m^3]on
inline scalar mixturePengRobinson::rho
(
const scalar p,
const scalar T,
const scalar rho0
) const
{
updateModelCoefficients(T);
return pengRobinson::rho(p,T,rho0);
}
//- Return density [kg/m^3]on
inline scalar mixturePengRobinson::rho(const scalar p, const scalar T) const
{
updateModelCoefficients(T);
return pengRobinson::rho(p,T);
}
//- Return compressibility drho/dp at T=constant [s^2/m^2]
inline scalar mixturePengRobinson::psi(const scalar rho, const scalar T) const
{
this->updateModelCoefficients(T);
return pengRobinson::psi(rho,T);
}
//- Return compression factor []
inline scalar mixturePengRobinson::Z
(
const scalar p,
const scalar T,
const scalar rho0
) const
{
updateModelCoefficients(T);
return pengRobinson::Z(p,T,rho0);
}
// * * * * * * * * * * * * * * * Member Operators * * * * * * * * * * * * * //
inline void mixturePengRobinson::operator+=(const mixturePengRobinson& pr)
{
//CL: Resizes the DynamicList if number of compontens is larger than the DynamicList
if((weigths.size()<numOfComp+1)||(mixtureComponents.size()<numOfComp+1))
{
weigths.setSize(2*numOfComp);
mixtureComponents.setSize(2*numOfComp);
}
//Cl: oldClass+=pr
//CL: Saving the object pointer and weigths of pr (which is @ mixtureComponents[0] and weigths[0]) at the numOfComp's value of the oldClass
weigths[numOfComp]=pr.weigths[0];
mixtureComponents[numOfComp]=pr.mixtureComponents[0];
pengRobinson::operator+=(pr);
//CL: increase number of Components by 1
numOfComp=numOfComp+1;
//CL: setting TSave=0--> results that all coefficients (a, da, d2a) are calculated in updateModelCoefficients
TSave=0.0;
singleComponent=0;
//CL:setting rho boundaries
rhoMin_=min(rhoMin_,pr.rhoMin_);
rhoMax_=max(rhoMax_,pr.rhoMax_);
//CL: calculating rho @ std
rhostd_=this->rho(this->Pstd(), this->Tstd(), this->Pstd()*this->W()/(this->Tstd()*this->R()));
}
// * * * * * * * * * * * * * * * Friend Operators * * * * * * * * * * * * * //
inline mixturePengRobinson operator+
(
const mixturePengRobinson& pr1,
const mixturePengRobinson& pr2
)
{
//CL:save both a pointer to the objector pr1 and pr2 and the weights in the new object
DynamicList<scalar> weigths=pr1.weigths;
DynamicList<mixturePengRobinson*> mixtureComponents=pr1.mixtureComponents;
//CL: Resizes the DynamicList if number of compontens is larger than the DynamicList
if((weigths.size()<pr1.numOfComp+1)||(mixtureComponents.size()<pr1.numOfComp+1))
{
weigths.setSize(2*pr1.numOfComp);
mixtureComponents.setSize(2*pr1.numOfComp);
}
//CL: Getting the new weigths and mixtureComponents lists,
//CL: Saving the object pointer and weigths of pr2 (which is @ mixtureComponents[0] and weigths[0]) at the numOfComp's value of the new object
weigths[pr1.numOfComp]=pr2.weigths[0];
mixtureComponents[pr1.numOfComp]=pr2.mixtureComponents[0];
return mixturePengRobinson(static_cast<const pengRobinson&>(pr1)+static_cast<const pengRobinson&>(pr2),
pr1.numOfComp+1, weigths, mixtureComponents, min(pr1.rhoMin_,pr2.rhoMin_), max(pr1.rhoMax_,pr2.rhoMax_));
}
inline mixturePengRobinson operator*
(
const scalar s,
const mixturePengRobinson& pr
)
{
//CL: saving the "concentraction" of the component of the mixture in the vector weights
//CL: saved at the Position "numOfComp-1"
DynamicList<scalar> weigths=pr.weigths;
weigths[pr.numOfComp-1]=s*pr.nMoles();
return mixturePengRobinson(s*static_cast<const pengRobinson&>(pr), pr.numOfComp, weigths,
pr.mixtureComponents,pr.a0_,pr.b_,pr.Tcrit_,pr.n_,pr.rhostd_, pr.rhoMin_, pr.rhoMax_);
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// ************************************************************************* //

154
src/thermophysicalModels/specie/equationOfState/cubicEquationOfState/mixtures/mixtureRedlichKwong/mixtureRedlichKwong.C
View file

@ -1,154 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright held by original author
\\/ 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 2 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, write to the Free Software Foundation,
Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
Description
Redlich Kwong equation of state for mixtures.
Author
Christian Lucas
Institut für Thermodynamik
Technische Universität Braunschweig
Germany
\*---------------------------------------------------------------------------*/
#include "mixtureRedlichKwong.H"
#include "IOstreams.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
mixtureRedlichKwong::mixtureRedlichKwong(Istream& is)
:
redlichKwong(is),
numOfComp(1),
//CL: no real gas mixture correction when stream constructor is used
realMixtureCorr_(false)
{
//CL: set size of weigths and mixtureComponents to 10
//CL: when more mixture components are used
//CL: size of the DynamicLis increases automatically
weigths.setSize(10);
mixtureComponents.setSize(10);
//Save a pointer of this object in the mixtureComponents array
mixtureComponents[0]=this;
is.check("mixtureRedlichKwong::mixtureRedlichKwong(Istream& is)");
}
//CL: Constructed needed in OpenFOAM 2.x.x
//CL: Code works fine, but compiling problem in OpenFOAM 1.6.ext
//CL: because specie has no constructor using dict
/*
mixtureRedlichKwong::mixtureRedlichKwong(const dictionary& dict)
:
redlichKwong(dict),
numOfComp(1),
realMixtureCorr_(dict.subDict("equationOfState").lookupOrDefault("realMixtureCorrection",false))
{
//CL: set size of weigths and mixtureComponents to 10
//CL: when more mixture components are used
//CL: size of the DynamicLis increases automatically
weigths.setSize(10);
mixtureComponents.setSize(10);
//Save a pointer of this object in the mixtureComponents array
mixtureComponents[0]=this;
if(realMixtureCorr_==true)
{
//CL: reads number of mixture components
//CL: this number is needed to calculate the number of correction coefficients
nCom_=dict.subDict("equationOfState").lookupOrDefault("numberOfComponents",1);
if (nCom_<1||nCom_==1)
{
FatalErrorIn
(
"mixtureRedlichKwong::mixtureRedlichKwong(const dictionary& dict)"
) << "mixture has only one component or number is not defined correctly, "
<< "recheck numberOfComponents in the thermophysicalProperties dict of your case"
<< abort(FatalError);
}
label i;
label j;
label k=0;
// CL: saves the number of mixtureCorrectionCoefficients needed for this mixture
// CL: need to be set to 1
label sizeOfVector_=1;
//CL: size of the vector depends on the number of components
//CL: determine the size of the vector
for (i=3; i<=nCom_;i++)
{
sizeOfVector_=sizeOfVector_+(i-1);
}
//CL: setting the size
realMixtureCorrCoef_.setSize(sizeOfVector_);
//CL: Reading the real mixture correction coefficients
//CL: Naming convention e.g for 3 component mixture:
//CL: mixtureCorrectionCoefficient_12, mixtureCorrectionCoefficient_13, mixtureCorrectionCoefficient_23
for(i=1;i<=nCom_-1;i++)
{
for(j=i+1;j<=nCom_;j++)
{
realMixtureCorrCoef_[k]=dict.subDict("equationOfState")
.lookupOrDefault("mixtureCorrectionCoefficient_" + Foam::name(i) + Foam::name(j),0.0);
k++;
}
}
}
}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
void Foam::mixtureRedlichKwong::write(Ostream& os) const
{
redlichKwong::write(os);
}
*/
// * * * * * * * * * * * * * * * Ostream Operator * * * * * * * * * * * * * //
Ostream& operator<<(Ostream& os, const mixtureRedlichKwong& rk)
{
os << static_cast<const redlichKwong&>(rk)<< token::SPACE ;
os.check("Ostream& operator<<(Ostream& os, const mixtureRedlichKwong& st)");
return os;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// ************************************************************************* //

198
src/thermophysicalModels/specie/equationOfState/cubicEquationOfState/mixtures/mixtureRedlichKwong/mixtureRedlichKwong.H
View file

@ -1,198 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright held by original author
\\/ 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 2 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, write to the Free Software Foundation,
Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
Class
Foam::mixtureRedlichKwong
Description
Redlich Kwong equation of state for multiple component mixtures using the van der Waals mixing rule
Mixing Rule --> see Paper
Title: Van der waals mixing rules for cublic equations of state. Applications for supercritical fluids extraction modelling
Authors: T.Y. Kwak and G.A. Mansoori
Journal: Chemical Engineering Science, Vol 41, No. 5, pp. 1303-1309, 1986
SourceFiles
mixtureRedlichKwongI.H
mixtureRedlichKwong.C
Author
Christian Lucas
Institut für Thermodynamik
Technische Universität Braunschweig
Germany
\*---------------------------------------------------------------------------*/
#ifndef mixtureRedlichKwong_H
#define mixtureRedlichKwong_H
#include "specie.H"
#include "autoPtr.H"
#include "redlichKwong.H"
#include "scalarList.H"
#include "DynamicList.H"
#include "label.H"
#include "List.H"
#include "dictionary.H"
#include "scalar.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
/*---------------------------------------------------------------------------*\
Class mixtureRedlichKwong Declaration
\*---------------------------------------------------------------------------*/
class mixtureRedlichKwong
:
public redlichKwong
{
// Private data
// Private functions
//CL: function updates the model coefficients (a,b) of the mixture
//CL: this is the function with the mixing rule
inline void updateModelCoefficients() const;
protected:
//CL: save the concentrations of each component of the mixture
//CL: needs to be multiplied by this->W() to get the molar fractions
mutable DynamicList<scalar> weigths;
//CL: saves a pointer to the pure component classes of the mixture
mutable DynamicList<mixtureRedlichKwong*> mixtureComponents;
//CL: counts the number of components
mutable label numOfComp;
//CL: Variables used in real gas mixture correction
//CL: If true, the real gas mixture correction is used
mutable bool realMixtureCorr_;
//CL: number of mixture components, needed to calculate the mixture correction factors needed
//CL: do not mistake this variable with numOfComp,
//CL: numOfComp is a counter to counts the number of components while the mixture is constructed
mutable label nCom_;
//CL: stores real mixture correction coefficients
mutable DynamicList<scalar> realMixtureCorrCoef_;
public:
// Constructors
//- Construct from components
//CL: needed for operator*
inline mixtureRedlichKwong
(
const redlichKwong& rK,
label numOfComp,
DynamicList<scalar> weigths,
DynamicList<mixtureRedlichKwong*> mixtureComponents,
scalar a,
scalar b,
scalar rhostd,
scalar rhoMin,
scalar rhoMax
);
//- Construct from components
//CL: needed for operator+
inline mixtureRedlichKwong
(
const redlichKwong& rK,
label numOfComp,
DynamicList<scalar> weigths,
DynamicList<mixtureRedlichKwong*> mixtureComponents,
scalar rhoMin,
scalar rhoMax
);
//- Construct from Istream
mixtureRedlichKwong(Istream&);
//- Construct from dictionary
//mixtureRedlichKwong(const dictionary& dict);
//- Construct as named copy
inline mixtureRedlichKwong(const word& name, const mixtureRedlichKwong&);
//- Construct and return a clone
inline autoPtr<mixtureRedlichKwong> clone() const;
// Selector from Istream
inline static autoPtr<mixtureRedlichKwong> New(Istream& is);
// I-O
//- Write to Ostream
//void write(Ostream& os) const;
// Member operators
inline void operator+=(const mixtureRedlichKwong&);
// Friend operators
inline friend mixtureRedlichKwong operator+
(
const mixtureRedlichKwong&,
const mixtureRedlichKwong&
);
inline friend mixtureRedlichKwong operator*
(
const scalar s,
const mixtureRedlichKwong&
);
// Ostream Operator
friend Ostream& operator<<(Ostream&, const mixtureRedlichKwong&);
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#include "mixtureRedlichKwongI.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

254
src/thermophysicalModels/specie/equationOfState/cubicEquationOfState/mixtures/mixtureRedlichKwong/mixtureRedlichKwongI.H
View file

@ -1,254 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright held by original author
\\/ 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 2 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, write to the Free Software Foundation,
Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
Author
Christian Lucas
Institut für Thermodynamik
Technische Universität Braunschweig
Germany
\*---------------------------------------------------------------------------*/
#include "mixtureRedlichKwong.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
// Construct from components
// CL: needed for operator*
inline mixtureRedlichKwong::mixtureRedlichKwong
(
const redlichKwong& rK,
label numOfComp,
DynamicList<scalar> weigths,
DynamicList<mixtureRedlichKwong*> mixtureComponents,
scalar a,
scalar b,
scalar rhostd,
scalar rhoMin,
scalar rhoMax
)
:
redlichKwong(rK),
weigths(weigths),
mixtureComponents(mixtureComponents),
numOfComp(numOfComp)
{
a_=a;
b_=b;
b2_=b*b;
b3_=pow(b,3);
b5_=pow(b,4);
rhostd_=rhostd;
//CL:setting rho boundaries
rhoMin_=rhoMin;
rhoMax_=rhoMax;
}
// Construct from components
// CL: needed for operator+
inline mixtureRedlichKwong::mixtureRedlichKwong
(
const redlichKwong& rK,
label numOfComp,
DynamicList<scalar> weigths,
DynamicList<mixtureRedlichKwong*> mixtureComponents,
scalar rhoMin,
scalar rhoMax
)
:
redlichKwong(rK),
weigths(weigths),
mixtureComponents(mixtureComponents),
numOfComp(numOfComp)
{
//CL: update model coefficients
updateModelCoefficients();
//CL:setting rho boundaries
rhoMin_=rhoMin;
rhoMax_=rhoMax;
//CL: calculating rho @ std
rhostd_=this->rho(this->Pstd(), this->Tstd(), this->Pstd()/(this->Tstd()*this->R()));
}
// Construct as named copy
inline mixtureRedlichKwong::mixtureRedlichKwong(const word& name, const mixtureRedlichKwong& rK)
:
redlichKwong(name, rK)
{}
// Construct and return a clone
inline autoPtr<mixtureRedlichKwong> mixtureRedlichKwong::clone() const
{
return autoPtr<mixtureRedlichKwong>(new mixtureRedlichKwong(*this));
}
// Selector from Istream
inline autoPtr<mixtureRedlichKwong> mixtureRedlichKwong::New(Istream& is)
{
return autoPtr<mixtureRedlichKwong>(new mixtureRedlichKwong(is));
}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
//CL: updates the coefficents of the model after the final construction of the mixture
//CL: uses the van der waals mixing rule
inline void mixtureRedlichKwong::updateModelCoefficients() const
{
label i,j,k;
a_=0;
b_=0;
for (i=0;i<numOfComp;i++)
{
//CL: getting a for the mixture
//CL: using van der waals mixing rule
for (j=0;j<numOfComp;j++)
{
//CL: Important: a is not a function of T
if(mixtureComponents[0]->realMixtureCorr_==true)
{
if(j!=i)
{
//CL: gives the position of the correction factor in the vector realMixtureCorrCoef_
k=i+j-1;
a_=a_+weigths[i]*weigths[j]*pow(mixtureComponents[i]->a()*mixtureComponents[j]->a(),0.5)
*this->W()*this->W()*(1-mixtureComponents[0]->realMixtureCorrCoef_[k]);
}
else
{
a_=a_+weigths[i]*weigths[j]*pow(mixtureComponents[i]->a()*mixtureComponents[j]->a(),0.5)*this->W()*this->W();
}
}
//CL: don't use mixture correction cofficients k_ij (see paper reference)
else
{
a_=a_+weigths[i]*weigths[j]*pow(mixtureComponents[i]->a()*mixtureComponents[j]->a(),0.5)*this->W()*this->W();
}
//CL: getting b for the mixture
//CL: using van der waals mixing rule
b_=b_+weigths[i]*mixtureComponents[i]->b()*this->W();
}
}
b2_=b_*b_;
b3_=pow(b_,3);
b5_=pow(b_,4);
}
// * * * * * * * * * * * * * * * Member Operators * * * * * * * * * * * * * //
inline void mixtureRedlichKwong::operator+=(const mixtureRedlichKwong& rK)
{
//CL: Resizes the DynamicList if number of compontens is larger than the DynamicList
if((weigths.size()<numOfComp+1)||(mixtureComponents.size()<numOfComp+1))
{
weigths.setSize(2*numOfComp);
mixtureComponents.setSize(2*numOfComp);
}
//Cl: oldClass+=rK
//CL: Saving the object pointer and weigths of rK (which is @ mixtureComponents[0] and weigths[0]) at the numOfComp's value of the oldClass
weigths[numOfComp]=rK.weigths[0];
mixtureComponents[numOfComp]=rK.mixtureComponents[0];
redlichKwong::operator+=(rK);
//CL: increase number of Components by 1
numOfComp=numOfComp+1;
//CL: update model coefficients
updateModelCoefficients();
//CL:setting rho boundaries
rhoMin_=min(rhoMin_,rK.rhoMin_);
rhoMax_=max(rhoMax_,rK.rhoMax_);
//CL: calculating rho @ std
rhostd_=this->rho(this->Pstd(), this->Tstd(), this->Pstd()*this->W()/(this->Tstd()*this->R()));
}
// * * * * * * * * * * * * * * * Friend Operators * * * * * * * * * * * * * //
inline mixtureRedlichKwong operator+
(
const mixtureRedlichKwong& rK1,
const mixtureRedlichKwong& rK2
)
{
//CL:save both a pointer to the objector rK1 and rK2 and the weights in the new object
DynamicList<scalar> weigths=rK1.weigths;
DynamicList<mixtureRedlichKwong*> mixtureComponents=rK1.mixtureComponents;
//CL: Resizes the DynamicList if number of compontens is larger than the DynamicList
if((weigths.size()<rK1.numOfComp+1)||(mixtureComponents.size()<rK1.numOfComp+1))
{
weigths.setSize(2*rK1.numOfComp);
mixtureComponents.setSize(2*rK1.numOfComp);
}
//CL: Getting the new weigths and mixtureComponents lists,
//CL: Saving the object pointer and weigths of rK2 (which is @ mixtureComponents[0] and weigths[0]) at the numOfComp's value of the new object
weigths[rK1.numOfComp]=rK2.weigths[0];
mixtureComponents[rK1.numOfComp]=rK2.mixtureComponents[0];
return mixtureRedlichKwong(static_cast<const redlichKwong&>(rK1)+static_cast<const redlichKwong&>(rK2),
rK1.numOfComp+1, weigths, mixtureComponents, min(rK1.rhoMin_,rK2.rhoMin_), max(rK1.rhoMax_,rK2.rhoMax_));
}
inline mixtureRedlichKwong operator*
(
const scalar s,
const mixtureRedlichKwong& rK
)
{
//CL: saving the "concentraction" of the component of the mixture in the vector weights
DynamicList<scalar> weigths=rK.weigths;
weigths[rK.numOfComp-1]=s*rK.nMoles();
return mixtureRedlichKwong(s*static_cast<const redlichKwong&>(rK), rK.numOfComp, weigths,
rK.mixtureComponents,rK.a_,rK.b_,rK.rhostd_, rK.rhoMin_, rK.rhoMax_);
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// ************************************************************************* //

157
src/thermophysicalModels/specie/equationOfState/cubicEquationOfState/mixtures/mixtureSoaveRedlichKwong/mixtureSoaveRedlichKwong.C
View file

@ -1,157 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright held by original author
\\/ 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 2 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, write to the Free Software Foundation,
Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
Description
Soave Redlich Kwong equation of state for mixtures.
Author
Christian Lucas
Institut für Thermodynamik
Technische Universität Braunschweig
Germany
\*---------------------------------------------------------------------------*/
#include "mixtureSoaveRedlichKwong.H"
#include "IOstreams.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
mixtureSoaveRedlichKwong::mixtureSoaveRedlichKwong(Istream& is)
:
soaveRedlichKwong(is),
singleComponent(1),
numOfComp(1),
//CL: no real gas mixture correction when stream constructor is used
realMixtureCorr_(false)
{
//CL: set size of weigths, mixtureComponents ... to 10,
//CL: when more mixture componentents are used
//CL: size of the DynamicLis increases automatically
weigths.setSize(10);
mixtureComponents.setSize(10);
//Save a pointer of this object in the mixtureComponents array
mixtureComponents[0]=this;
is.check("mixtureSoaveRedlichKwong::mixtureSoaveRedlichKwong(Istream& is)");
}
//CL: Constructed needed in OpenFOAM 2.x.x
//CL: Code works fine, but compiling problem in OpenFOAM 1.6.ext
//CL: because specie has no constructor using dict
/*
mixtureSoaveRedlichKwong::mixtureSoaveRedlichKwong(const dictionary& dict)
:
soaveRedlichKwong(dict),
numOfComp(1),
singleComponent(1),
realMixtureCorr_(dict.subDict("equationOfState").lookupOrDefault("realMixtureCorrection",false))
{
//CL: set size of weigths, mixtureComponents ... to 10,
//CL: when more mixture componentents are used
//CL: size of the DynamicLis increases automatically
weigths.setSize(10);
mixtureComponents.setSize(10);
//Save a pointer of this object in the mixtureComponents array
mixtureComponents[0]=this;
if(realMixtureCorr_==true)
{
//CL: reads number of mixture components
//CL: this number is needed to calculate the number of correction coefficients
nCom_=dict.subDict("equationOfState").lookupOrDefault("numberOfComponents",1);
if (nCom_<1||nCom_==1)
{
FatalErrorIn
(
"mixturePengRobinson::mixturePengRobinson(const dictionary& dict)"
) << "mixture has only one component or number is not defined correctly, "
<< "recheck numberOfComponents in the thermophysicalProperties dict of your case"
<< abort(FatalError);
}
label i;
label j;
label k=0;
// CL: saves the number of mixtureCorrectionCoefficients needed for this mixture
// CL: need to be set to 1
label sizeOfVector_=1;
//CL: size of the vector depends on the number of components
//CL: determine the size of the vector
for (i=3; i<=nCom_;i++)
{
sizeOfVector_=sizeOfVector_+(i-1);
}
//CL: setting the size
realMixtureCorrCoef_.setSize(sizeOfVector_);
//CL: Reading the real mixture correction coefficients
//CL: Naming convention e.g for 3 component mixture:
//CL: mixtureCorrectionCoefficient_12, mixtureCorrectionCoefficient_13, mixtureCorrectionCoefficient_23
for(i=1;i<=nCom_-1;i++)
{
for(j=i+1;j<=nCom_;j++)
{
realMixtureCorrCoef_[k]=dict.subDict("equationOfState")
.lookupOrDefault("mixtureCorrectionCoefficient_" + Foam::name(i) + Foam::name(j),0.0);
k++;
}
}
}
}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
void Foam::mixtureSoaveRedlichKwong::write(Ostream& os) const
{
soaveRedlichKwong::write(os);
}
*/
// * * * * * * * * * * * * * * * Ostream Operator * * * * * * * * * * * * * //
Ostream& operator<<(Ostream& os, const mixtureSoaveRedlichKwong& srk)
{
os << static_cast<const soaveRedlichKwong&>(srk)<< token::SPACE ;
os.check("Ostream& operator<<(Ostream& os, const mixtureSoaveRedlichKwong& st)");
return os;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// ************************************************************************* //

273
src/thermophysicalModels/specie/equationOfState/cubicEquationOfState/mixtures/mixtureSoaveRedlichKwong/mixtureSoaveRedlichKwong.H
View file

@ -1,273 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright held by original author
\\/ 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 2 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, write to the Free Software Foundation,
Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
Class
Foam::mixtureSoaveRedlichKwong
Description
Soave Redlich Kwong equation of state for multiple component mixtures using the van der Waals mixing rule
Mixing Rule --> see Paper
Title: Van der eaals mixting rules for cublic equations of state. Applications for supercritical fluids extraction modelling
Authors: T.Y. Kwak and G.A. Mansoori
Journal: Chemical Engineering Science, Vol 41, No. 5, pp. 1303-1309, 1986
SourceFiles
mixtureSoaveRedlichKwongI.H
mixtureSoaveRedlichKwong.C
Author
Christian Lucas
Institut für Thermodynamik
Technische Universität Braunschweig
Germany
\*---------------------------------------------------------------------------*/
#ifndef mixtureSoaveRedlichKwong_H
#define mixtureSoaveRedlichKwong_H
#include "specie.H"
#include "autoPtr.H"
#include "soaveRedlichKwong.H"
#include "DynamicList.H"
#include "scalarList.H"
#include "label.H"
#include "List.H"
#include "dictionary.H"
#include "scalar.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
/*---------------------------------------------------------------------------*\
Class mixtureSoaveRedlichKwong Declaration
\*---------------------------------------------------------------------------*/
class mixtureSoaveRedlichKwong
:
public soaveRedlichKwong
{
protected:
//CL: bool used to make sure the model behaves like a single component model when needed
//CL: this is needed during the construction of the mixture
mutable bool singleComponent;
//CL: save the concentrations of each component of the mixture
//CL: needs to be multiplied by this->W() to get the molar fractions
mutable DynamicList<scalar> weigths;
//CL: saves a pointer to the pure component classes of the mixture
mutable DynamicList<mixtureSoaveRedlichKwong*> mixtureComponents;
//CL: counts the number of components
mutable label numOfComp;
//Protected functions
//CL: function updates the coefficients (aSave, daSave, d2aSave, b) of the mixture
inline void updateModelCoefficients(const scalar T) const;
//CL: Variables used in real gas mixture correction
//CL: If true, the real gas mixture correction is used
mutable bool realMixtureCorr_;
//CL: number of mixture components, needed to calculate the mixture correction factors needed
//CL: do not mistake this variable with numOfComp,
//CL: numOfComp is a counter to counts the number of components while the mixture is constructed
mutable label nCom_;
//CL: stores real mixture correction coefficients
mutable DynamicList<scalar> realMixtureCorrCoef_;
public:
// Constructors
//- Construct from components
//CL: needed for operator*
inline mixtureSoaveRedlichKwong
(
const soaveRedlichKwong& srk,
label numOfComp,
DynamicList<scalar> weigths,
DynamicList<mixtureSoaveRedlichKwong*> mixtureComponents,
scalar a0,
scalar b,
scalar Tcrit,
scalar n,
scalar rhostd,
scalar rhoMin,
scalar rhoMax
);
//- Construct from components
//CL: needed for operator+
inline mixtureSoaveRedlichKwong
(
const soaveRedlichKwong& srk,
label numOfComp,
DynamicList<scalar> weigths,
DynamicList<mixtureSoaveRedlichKwong*> mixtureComponents,
scalar rhoMin,
scalar rhoMax
);
//- Construct from Istream
mixtureSoaveRedlichKwong(Istream&);
//- Construct from dictionary
//mixtureSoaveRedlichKwong(const dictionary& dict);
//- Construct as named copy
inline mixtureSoaveRedlichKwong(const word& name, const mixtureSoaveRedlichKwong&);
//- Construct and return a clone
inline autoPtr<mixtureSoaveRedlichKwong> clone() const;
// Selector from Istream
inline static autoPtr<mixtureSoaveRedlichKwong> New(Istream& is);
// Member functions
//CL: mixture coefficience
inline scalar amix(const scalar T, const label i, const label j) const;
inline scalar dadTmix(const scalar T, const label i, const label j) const;
inline scalar d2adT2mix(const scalar T, const label i, const label j) const;
inline scalar p(const scalar rho, const scalar T) const;
//first order derivatives
inline scalar dpdv(const scalar rho, const scalar T) const;
inline scalar dpdT(const scalar rho, const scalar T) const;
inline scalar dvdT(const scalar rho, const scalar T) const;
inline scalar dvdp(const scalar rho, const scalar T) const;
inline scalar isobarExpCoef(const scalar rho, const scalar T) const;
inline scalar isothermalCompressiblity
(
const scalar rho,
const scalar T
) const;
// Used for cv
inline scalar integral_d2pdT2_dv
(
const scalar rho,
const scalar T
) const;
//Used for internal Energy
inline scalar integral_p_dv(const scalar rho, const scalar T) const;
// Used for Entropy
inline scalar integral_dpdT_dv(const scalar rho, const scalar T) const;
// second order derivatives, not Used At The Moment
inline scalar d2pdv2(const scalar rho, const scalar T) const;
inline scalar d2pdT2(const scalar rho, const scalar T) const;
inline scalar d2pdvdT(const scalar rho, const scalar T) const;
inline scalar d2vdT2(const scalar rho, const scalar T) const;
//- Return density [kg/m^3]
// rho0 is the starting point of the newton solver used to calculate rho
inline scalar rho
(
const scalar p,
const scalar T,
const scalar rho0
) const;
inline scalar rho(const scalar p, const scalar T) const;
//- Return compressibility drho/dp at T=constant [s^2/m^2]
inline scalar psi(const scalar rho, const scalar T) const;
//- Return compression factor []
inline scalar Z
(
const scalar p,
const scalar T,
const scalar rho0
) const;
// I-O
//- Write to Ostream
//void write(Ostream& os) const;
// Member operators
inline void operator+=(const mixtureSoaveRedlichKwong&);
// Friend operators
inline friend mixtureSoaveRedlichKwong operator+
(
const mixtureSoaveRedlichKwong&,
const mixtureSoaveRedlichKwong&
);
inline friend mixtureSoaveRedlichKwong operator*
(
const scalar s,
const mixtureSoaveRedlichKwong&
);
// Ostream Operator
friend Ostream& operator<<(Ostream&, const mixtureSoaveRedlichKwong&);
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#include "mixtureSoaveRedlichKwongI.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

498
src/thermophysicalModels/specie/equationOfState/cubicEquationOfState/mixtures/mixtureSoaveRedlichKwong/mixtureSoaveRedlichKwongI.H
View file

@ -1,498 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright held by original author
\\/ 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 2 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, write to the Free Software Foundation,
Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
Author
Christian Lucas
Institut für Thermodynamik
Technische Universität Braunschweig
Germany
\*---------------------------------------------------------------------------*/
#include "mixtureSoaveRedlichKwong.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
// Construct from components
// CL: needed for operator*
inline mixtureSoaveRedlichKwong::mixtureSoaveRedlichKwong
(
const soaveRedlichKwong& srk,
label numOfComp,
DynamicList<scalar> weigths,
DynamicList<mixtureSoaveRedlichKwong*> mixtureComponents,
scalar a0,
scalar b,
scalar Tcrit,
scalar n,
scalar rhostd,
scalar rhoMin,
scalar rhoMax
)
:
soaveRedlichKwong(srk),
singleComponent(1),
weigths(weigths),
mixtureComponents(mixtureComponents),
numOfComp(numOfComp)
{
a0_=a0;
b_=b;
Tcrit_=Tcrit;
n_=n;
rhostd_=rhostd;
rhoMin_=rhoMin;
rhoMax_=rhoMax;
}
// Construct from components
// CL: needed for operator+
inline mixtureSoaveRedlichKwong::mixtureSoaveRedlichKwong
(
const soaveRedlichKwong& srk,
label numOfComp,
DynamicList<scalar> weigths,
DynamicList<mixtureSoaveRedlichKwong*> mixtureComponents,
scalar rhoMin,
scalar rhoMax
)
:
soaveRedlichKwong(srk),
singleComponent(0),
weigths(weigths),
mixtureComponents(mixtureComponents),
numOfComp(numOfComp)
{
TSave=0.0;
rhoMin_=rhoMin;
rhoMax_=rhoMax;
rhostd_=this->rho(this->Pstd(), this->Tstd(), this->Pstd()*this->W()/(this->Tstd()*this->R()));
}
// Construct as named copy
inline mixtureSoaveRedlichKwong::mixtureSoaveRedlichKwong(const word& name, const mixtureSoaveRedlichKwong& srk)
:
soaveRedlichKwong(name, srk)
{}
// Construct and return a clone
inline autoPtr<mixtureSoaveRedlichKwong> mixtureSoaveRedlichKwong::clone() const
{
return autoPtr<mixtureSoaveRedlichKwong>(new mixtureSoaveRedlichKwong(*this));
}
// Selector from Istream
inline autoPtr<mixtureSoaveRedlichKwong> mixtureSoaveRedlichKwong::New(Istream& is)
{
return autoPtr<mixtureSoaveRedlichKwong>(new mixtureSoaveRedlichKwong(is));
}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
//CL: updates the coefficents of the model after the final construction of the mixture
//CL: uses the van der waals mixing rule
inline void mixtureSoaveRedlichKwong::updateModelCoefficients(const scalar T) const
{
// CL: to make sure that the coefficents are only updated if the mixture has more than 1 component
if (singleComponent==0)
{
// Checking if the mixture coefficient were already calculated for this temperature
if(TSave!=T)
{
label i,j,k;
aSave=0;
daSave=0;
d2aSave=0;
b_=0;
for (i=0;i<numOfComp;i++)
{
//CL: getting a(T), dadT(T) and d2adT2(T) for the mixture
//CL: using van der waals mixing rule
for (j=0;j<numOfComp;j++)
{
//CL: use mixture correction cofficients k_ij (see paper reference)
if(mixtureComponents[0]->realMixtureCorr_==true)
{
// first and second order temperature derivative of the van der waals mixing rule for a(T)
if(i==j)
{
aSave=aSave+amix(T,i,j);
daSave=daSave+dadTmix(T,i,j);
d2aSave=d2aSave+d2adT2mix(T,i,j);
}
else
{
//CL: gives the position of the correction factor in the vector realMixtureCorrCoef_
k=i+j-1;
aSave=aSave+amix(T,i,j)*(1-mixtureComponents[0]->realMixtureCorrCoef_[k]);
daSave=daSave+dadTmix(T,i,j)*(1-mixtureComponents[0]->realMixtureCorrCoef_[k]);
d2aSave=d2aSave+d2adT2mix(T,i,j)*(1-mixtureComponents[0]->realMixtureCorrCoef_[k]);
}
}
else
{
aSave=aSave+amix(T,i,j);
daSave=daSave+dadTmix(T,i,j);
d2aSave=d2aSave+d2adT2mix(T,i,j);
}
}
//CL: getting b for the mixture
//CL: using van der waals mixing rule
b_=b_+weigths[i]*mixtureComponents[i]->b()*this->W();
//CL: saving the temperature at which the mixture coefficients are valid
TSave=T;
}
b2_=b_*b_;
b3_=pow(b_,3);
b5_=pow(b_,5);
}
}
}
//CL: returns part of "a" for the mixture
//CL: Important: a is a function of T
//CL: TODO: RECHECK
inline scalar mixtureSoaveRedlichKwong::amix(const scalar T, const label i, const label j) const
{
return weigths[i]*weigths[j]*this->W()*this->W()*
pow(mixtureComponents[i]->a(T)*mixtureComponents[j]->a(T),0.5);
}
//CL: returns part of the first derivative of "a" for the mixture
//CL: Important: a is a function of T
//CL: TODO: RECHECK
inline scalar mixtureSoaveRedlichKwong::dadTmix(const scalar T, const label i, const label j) const
{
return 0.5*weigths[i]*weigths[j]*this->W()*this->W()
*(pow(mixtureComponents[i]->a(T)/mixtureComponents[j]->a(T),0.5)*mixtureComponents[j]->dadT(T)
+pow(mixtureComponents[j]->a(T)/mixtureComponents[i]->a(T),0.5)*mixtureComponents[i]->dadT(T));
}
//CL: returns part of the second derivative of "a" for the mixture
//CL: Important: a is a function of T
//CL: TODO: RECHECK
inline scalar mixtureSoaveRedlichKwong::d2adT2mix(const scalar T, const label i, const label j) const
{
return 0.5*weigths[i]*weigths[j]*this->W()*this->W()*
(
pow(mixtureComponents[i]->a(T)/mixtureComponents[j]->a(T),0.5)*mixtureComponents[j]->d2adT2(T)
+pow(mixtureComponents[j]->a(T)/mixtureComponents[i]->a(T),0.5)*mixtureComponents[i]->d2adT2(T)
+pow(mixtureComponents[i]->a(T)*mixtureComponents[j]->a(T),-0.5)
*mixtureComponents[i]->dadT(T)*mixtureComponents[j]->dadT(T)
-0.5*pow(mixtureComponents[i]->a(T)/mixtureComponents[j]->a(T),0.5)*1/mixtureComponents[j]->a(T)
*pow(mixtureComponents[j]->dadT(T),2)
-0.5*pow(mixtureComponents[j]->a(T)/mixtureComponents[i]->a(T),0.5)
*1/mixtureComponents[i]->a(T)*pow(mixtureComponents[i]->dadT(T),2)
);
}
//returns the pressure for a given density and temperature
inline scalar mixtureSoaveRedlichKwong::p(const scalar rho, const scalar T) const
{
updateModelCoefficients(T);
return soaveRedlichKwong::p(rho,T);
}
//Real deviative dp/dv at constant temperature
//(molar values)
inline scalar mixtureSoaveRedlichKwong::dpdv(const scalar rho, const scalar T) const
{
updateModelCoefficients(T);
return soaveRedlichKwong::dpdv(rho,T);
}
//Real deviative dp/dT at constant molar volume
//(molar values)
inline scalar mixtureSoaveRedlichKwong::dpdT(const scalar rho, const scalar T) const
{
updateModelCoefficients(T);
return soaveRedlichKwong::dpdT(rho,T);
}
//Real deviative dv/dT at constant pressure
//(molar values)
inline scalar mixtureSoaveRedlichKwong::dvdT(const scalar rho, const scalar T) const
{
updateModelCoefficients(T);
return soaveRedlichKwong::dvdT(rho,T);
}
//Real deviative dv/dp at constant temperature
//(molar values)
inline scalar mixtureSoaveRedlichKwong::dvdp(const scalar rho, const scalar T) const
{
updateModelCoefficients(T);
return soaveRedlichKwong::dvdp(rho,T);
}
//needed to calculate the internal energy
//(molar values)
inline scalar mixtureSoaveRedlichKwong::integral_p_dv
(
const scalar rho,
const scalar T
) const
{
updateModelCoefficients(T);
return soaveRedlichKwong::integral_p_dv(rho,T);
}
//needed to calculate the entropy
//(molar values)
inline scalar mixtureSoaveRedlichKwong::integral_dpdT_dv
(
const scalar rho,
const scalar T
) const
{
updateModelCoefficients(T);
return soaveRedlichKwong::integral_dpdT_dv(rho,T);
}
//(molar values)
inline scalar mixtureSoaveRedlichKwong::d2pdT2(const scalar rho, const scalar T) const
{
updateModelCoefficients(T);
return soaveRedlichKwong::d2pdT2(rho,T);
}
//(molar values)
inline scalar mixtureSoaveRedlichKwong::d2pdv2(const scalar rho, const scalar T) const
{
updateModelCoefficients(T);
return soaveRedlichKwong::d2pdv2(rho,T);
}
//(molar values)
inline scalar mixtureSoaveRedlichKwong::d2vdT2(const scalar rho, const scalar T) const
{
updateModelCoefficients(T);
return soaveRedlichKwong::d2vdT2(rho,T);
}
//(molar values)
inline scalar mixtureSoaveRedlichKwong::d2pdvdT(const scalar rho, const scalar T) const
{
updateModelCoefficients(T);
return soaveRedlichKwong::d2pdvdT(rho,T);
}
// the result of this intergal is needed for the nasa based cp polynomial
//(molar values)
inline scalar mixtureSoaveRedlichKwong::integral_d2pdT2_dv
(
const scalar rho,
const scalar T
) const
{
updateModelCoefficients(T);
return soaveRedlichKwong::integral_d2pdT2_dv(rho,T);
}
//Isobar expansion Coefficent beta = 1/v (dv/dt) at constant p
//(molar values)
inline scalar mixtureSoaveRedlichKwong::isobarExpCoef
(
const scalar rho,
const scalar T
) const
{
updateModelCoefficients(T);
return soaveRedlichKwong::isobarExpCoef(rho,T);
}
//isothemal compressiblity kappa (not Thermal conductivity)
//(molar values)
inline scalar mixtureSoaveRedlichKwong::isothermalCompressiblity
(
const scalar rho,
const scalar T
) const
{
updateModelCoefficients(T);
return soaveRedlichKwong::isothermalCompressiblity(rho,T);
}
//- Return density [kg/m^3]on
inline scalar mixtureSoaveRedlichKwong::rho
(
const scalar p,
const scalar T,
const scalar rho0
) const
{
updateModelCoefficients(T);
return soaveRedlichKwong::rho(p,T,rho0);
}
//- Return density [kg/m^3]on
inline scalar mixtureSoaveRedlichKwong::rho(const scalar p, const scalar T) const
{
updateModelCoefficients(T);
return soaveRedlichKwong::rho(p,T);
}
//- Return compressibility drho/dp at T=constant [s^2/m^2]
inline scalar mixtureSoaveRedlichKwong::psi(const scalar rho, const scalar T) const
{
this->updateModelCoefficients(T);
return soaveRedlichKwong::psi(rho,T);
}
//- Return compression factor []
inline scalar mixtureSoaveRedlichKwong::Z
(
const scalar p,
const scalar T,
const scalar rho0
) const
{
updateModelCoefficients(T);
return soaveRedlichKwong::Z(p,T,rho0);
}
// * * * * * * * * * * * * * * * Member Operators * * * * * * * * * * * * * //
inline void mixtureSoaveRedlichKwong::operator+=(const mixtureSoaveRedlichKwong& srk)
{
//CL: Resizes the DynamicList if number of compontens is larger than the DynamicList
if((weigths.size()<numOfComp+1)||(mixtureComponents.size()<numOfComp+1))
{
weigths.setSize(2*numOfComp);
mixtureComponents.setSize(2*numOfComp);
}
//Cl: oldClass+=srk
//CL: Saving the object pointer and weigths of pr (which is @ mixtureComponents[0] and weigths[0]) at the numOfComp's value of the oldClass
weigths[numOfComp]=srk.weigths[0];
mixtureComponents[numOfComp]=srk.mixtureComponents[0];
soaveRedlichKwong::operator+=(srk);
//CL: increase number of Components by 1
numOfComp=numOfComp+1;
//CL: setting TSave=0--> results that all coefficients (a, da, d2a) are calculated in updateModelCoefficients
TSave=0.0;
singleComponent=0;
//CL:setting rho boundaries
rhoMin_=min(rhoMin_,srk.rhoMin_);
rhoMax_=max(rhoMax_,srk.rhoMax_);
//CL: calculating rho @ std
rhostd_=this->rho(this->Pstd(), this->Tstd(), this->Pstd()*this->W()/(this->Tstd()*this->R()));
}
// * * * * * * * * * * * * * * * Friend Operators * * * * * * * * * * * * * //
inline mixtureSoaveRedlichKwong operator+
(
const mixtureSoaveRedlichKwong& srk1,
const mixtureSoaveRedlichKwong& srk2
)
{
//CL:save both a pointer to the objector srk1 and srk2 and the weights in the new object
DynamicList<scalar> weigths=srk1.weigths;
DynamicList<mixtureSoaveRedlichKwong*> mixtureComponents=srk1.mixtureComponents;
//CL: Resizes the DynamicList if number of compontens is larger than the DynamicList
if((weigths.size()<srk1.numOfComp+1)||(mixtureComponents.size()<srk1.numOfComp+1))
{
weigths.setSize(2*srk1.numOfComp);
mixtureComponents.setSize(2*srk1.numOfComp);
}
//CL: Getting the new weigths and mixtureComponents lists,
//CL: Saving the object pointer and weigths of srk2 (which is @ mixtureComponents[0] and weigths[0]) at the numOfComp's value of the new object
weigths[srk1.numOfComp]=srk2.weigths[0];
mixtureComponents[srk1.numOfComp]=srk2.mixtureComponents[0];
return mixtureSoaveRedlichKwong(static_cast<const soaveRedlichKwong&>(srk1)+static_cast<const soaveRedlichKwong&>(srk2),
srk1.numOfComp+1, weigths, mixtureComponents, min(srk1.rhoMin_,srk2.rhoMin_), max(srk1.rhoMax_,srk2.rhoMax_));
}
inline mixtureSoaveRedlichKwong operator*
(
const scalar s,
const mixtureSoaveRedlichKwong& srk
)
{
//CL: saving the "concentraction" of the component of the mixture in the vector weights
//CL: saved at the Position "numOfComp-1"
DynamicList<scalar> weigths=srk.weigths;
weigths[srk.numOfComp-1]=s*srk.nMoles();
return mixtureSoaveRedlichKwong(s*static_cast<const soaveRedlichKwong&>(srk), srk.numOfComp, weigths,
srk.mixtureComponents,srk.a0_,srk.b_,srk.Tcrit_,srk.n_,srk.rhostd_, srk.rhoMin_, srk.rhoMax_);
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// ************************************************************************* //

76
src/thermophysicalModels/specie/equationOfState/cubicEquationOfState/pengRobinson/pengRobinson.C
View file

@ -1,26 +1,25 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / F ield | foam-extend: Open Source CFD
\\ / O peration |
\\ / A nd | Copyright held by original author
\\ / A nd | For copyright notice see file Copyright
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
This file is part of foam-extend.
OpenFOAM is free software; you can redistribute it and/or modify it
foam-extend 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 2 of the License, or (at your
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.
foam-extend 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, write to the Free Software Foundation,
Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
along with foam-extend. If not, see <http://www.gnu.org/licenses/>.
Description
Peng Robinson equation of state.
@ -36,32 +35,31 @@ Germany
#include "pengRobinson.H"
#include "IOstreams.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
pengRobinson::pengRobinson(Istream& is)
Foam::pengRobinson::pengRobinson(Istream& is)
:
specie(is),
pcrit_(readScalar(is)),
Tcrit_(readScalar(is)),
azentricFactor_(readScalar(is)),
n_(0.37464+1.54226*azentricFactor_-0.26992*pow(azentricFactor_,2)),
a0_(0.457235*pow(this->RR(), 2)*pow(Tcrit_, 2)/pcrit_),
b_(0.077796*this->RR()*Tcrit_/pcrit_),
n_(0.37464 + 1.54226*azentricFactor_ - 0.26992*pow(azentricFactor_, 2)),
b2_(b_*b_),
b3_(b2_*b_),
b4_(b3_*b_),
b5_(b4_*b_),
b6_(b5_*b_),
//CL: Only uses the default values
b2_(pow(b_,2)),
b3_(pow(b_,3)),
b4_(pow(b_,4)),
b5_(pow(b_,5)),
rhoMax_(1500),
rhoMin_(1e-3),
TSave(0.0),
rhoMax_(1500),
// Starting GUESS for the density by ideal gas law
rhostd_(this->rho(this->Pstd(),this->Tstd(),this->Pstd()/(this->Tstd()*this->R())))
rhostd_(this->rho(this->Pstd(), this->Tstd(), this->Pstd()/(this->Tstd()*this->R()))),
aSave(0.0),
daSave(0.0),
d2aSave(0.0),
TSave(0.0)
{
is.check("pengRobinson::pengRobinson(Istream& is)");
}
@ -78,19 +76,25 @@ pengRobinson::pengRobinson(const dictionary& dict)
//CL: rhoMin and rhoMax are only used as boundaries for the bisection methode (see rho function)
//CL: important: rhoMin and rhoMax are not used as boundary for the newton solver
//CL: therefore, rho can be larger than rhoMax and smaller than rhoMin
rhoMin_(dict.subDict("equationOfState").lookupOrDefault("rhoMin",1e-3)),
rhoMax_(dict.subDict("equationOfState").lookupOrDefault("rhoMax",1500)),
a0_(0.457235*pow(this->RR(), 2)*pow(Tcrit_, 2)/pcrit_),
b_(0.077796*this->RR()*Tcrit_/pcrit_),
n_(0.37464 + 1.54226*azentricFactor_ - 0.26992*pow(azentricFactor_, 2)),
b2_(b_*b_),
b3_(b2_*b_),
b4_(b3_*b_),
b5_(b4_*b_),
b6_(b5_*b_),
rhoMin_(dict.subDict("equationOfState").lookupOrDefault("rhoMin",1e-3)),
rhoMax_(dict.subDict("equationOfState").lookupOrDefault("rhoMax",1500)),
aSave(0.0),
daSave(0.0),
d2aSave(0.0),
TSave(0.0),
b2_(pow(b_,2)),
b3_(pow(b_,3)),
b4_(pow(b_,4)),
b5_(pow(b_,5)),
rhostd_(this->rho(Pstd,Tstd,Pstd/(Tstd*this->R())))
// Starting GUESS for the density by ideal gas law
rhostd_(this->rho(this->Pstd(), this->Tstd(), this->Pstd()/(this->Tstd()*this->R())))
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
void Foam::pengRobinson::write(Ostream& os) const
@ -107,9 +111,11 @@ void Foam::pengRobinson::write(Ostream& os) const
os << indent << dict.dictName() << dict;
}
*/
// * * * * * * * * * * * * * * * Ostream Operator * * * * * * * * * * * * * //
Ostream& operator<<(Ostream& os, const pengRobinson& pr)
Foam::Ostream& Foam::operator<<(Ostream& os, const pengRobinson& pr)
{
os << static_cast<const specie&>(pr)<< token::SPACE
<< pr.pcrit_ << tab<< pr.Tcrit_<< tab << pr.azentricFactor_;
@ -119,8 +125,4 @@ Ostream& operator<<(Ostream& os, const pengRobinson& pr)
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// ************************************************************************* //

98
src/thermophysicalModels/specie/equationOfState/cubicEquationOfState/pengRobinson/pengRobinson.H
View file

@ -1,26 +1,25 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / F ield | foam-extend: Open Source CFD
\\ / O peration |
\\ / A nd | Copyright held by original author
\\ / A nd | For copyright notice see file Copyright
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
This file is part of foam-extend.
OpenFOAM is free software; you can redistribute it and/or modify it
foam-extend 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 2 of the License, or (at your
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.
foam-extend 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, write to the Free Software Foundation,
Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
along with foam-extend. If not, see <http://www.gnu.org/licenses/>.
Class
Foam::pengRobinson
@ -53,14 +52,13 @@ Germany
#include "specie.H"
#include "autoPtr.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
/*---------------------------------------------------------------------------*\
Class perfectGas Declaration
Class Peng Robinson Declaration
\*---------------------------------------------------------------------------*/
class pengRobinson
@ -68,27 +66,30 @@ class pengRobinson
public specie
{
protected:
// Protected data
// private data
//CL: data at critical point
scalar pcrit_;
scalar Tcrit_;
scalar azentricFactor_;
//-Peng Robinson factors
scalar n_;
scalar a0_;
mutable scalar b_;
scalar b_;
scalar n_;
//CL: pow of constants b_ used in the code e.g. b2_=b*b;
mutable scalar b2_;
mutable scalar b3_;
mutable scalar b4_;
mutable scalar b5_;
scalar b2_;
scalar b3_;
scalar b4_;
scalar b5_;
scalar b6_;
//CL: rhoMin and rhoMax are only used as boundaries for the bisection methode (see rho function)
scalar rhoMax_;
scalar rhoMin_;
scalar rhoMax_;
//- Density @STD, initialise after a0, b!
scalar rhostd_;
//CL: Variables to save the values of a, dadT and d2adT2 of the mixture
//CL: Variables must corrected for changing temperatures
@ -99,16 +100,12 @@ protected:
//CL: save the temperature for which the save coefficients (amix,dadTmix,d2adT2mix) are correct
mutable scalar TSave;
//- Density @STD, initialise after a0, b!
mutable scalar rhostd_;
//Protected functions
//CL: function updates the coefficients (aSave, daSave, d2aSave)
inline void updateModelCoefficients(const scalar T) const;
public:
// Constructors
//- Construct from components
@ -133,8 +130,22 @@ public:
inline static autoPtr<pengRobinson> New(Istream& is);
// Member functions
//Return Peng Robinson factors
inline scalar a0() const;
inline scalar b() const;
inline scalar n() const;
inline scalar rhostd() const;
inline scalar rhoMin() const;
inline scalar rhoMax() const;
inline scalar Tcrit() const;
//CL: Model coefficient a(T)
inline scalar a(const scalar T) const;
@ -144,22 +155,6 @@ public:
//CL: second order temperature deriviative of model coefficient a(T)
inline scalar d2adT2(const scalar T) const;
//Return Peng Robinson factors
inline scalar a0()const;
inline scalar b()const;
inline scalar n()const;
//CL: return power of constants b_
inline scalar b2()const;
inline scalar b3()const;
inline scalar b4()const;
inline scalar b5()const;
//CL: Equation of state
inline scalar p(const scalar rho, const scalar T) const;
@ -201,10 +196,18 @@ public:
inline scalar d2vdT2(const scalar rho, const scalar T) const;
//CL: Used for internal Energy
inline scalar integral_p_dv(const scalar rho,const scalar T) const;
inline scalar integral_p_dv
(
const scalar rho,
const scalar T
) const;
//CL: Used for Entropy
inline scalar integral_dpdT_dv(const scalar rho,const scalar T) const;
inline scalar integral_dpdT_dv
(
const scalar rho,
const scalar T
) const;
//- Return density [kg/m^3]
// rho0 is the starting point of the newton solver used to calculate rho
@ -239,11 +242,7 @@ public:
inline void operator+=(const pengRobinson&);
/*
inline void operator-=(const pengRobinson&);
inline void operator*=(const scalar);
*/
// Friend operators
@ -260,6 +259,7 @@ public:
const pengRobinson&
);
// Ostream Operator
friend Ostream& operator<<(Ostream&, const pengRobinson&);

233
src/thermophysicalModels/specie/equationOfState/cubicEquationOfState/pengRobinson/pengRobinsonI.H
View file

@ -57,15 +57,21 @@ inline pengRobinson::pengRobinson(const word& name, const pengRobinson& pr)
pcrit_(pr.pcrit_),
Tcrit_(pr.Tcrit_),
azentricFactor_(pr.azentricFactor_),
n_(pr.n_),
a0_(pr.a0_),
b_(pr.b_),
n_(pr.n_),
b2_(pr.b2_),
b3_(pr.b3_),
b4_(pr.b4_),
b5_(pr.b5_),
TSave(0),
rhostd_(pr.rhostd_)
b6_(pr.b6_),
rhoMin_(pr.rhoMin_),
rhoMax_(pr.rhoMax_),
rhostd_(pr.rhostd_),
aSave(0.0),
daSave(0.0),
d2aSave(0.0),
TSave(0.0)
{}
@ -83,25 +89,43 @@ inline autoPtr<pengRobinson> pengRobinson::New(Istream& is)
}
// * * * * * * * * * * * * * Member Functions * * * * * * * * * * * //
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
inline scalar pengRobinson::rhostd() const
{
return rhostd_;
}
inline scalar pengRobinson::rhoMin() const
{
return rhoMin_;
}
inline scalar pengRobinson::rhoMax() const
{
return rhoMax_;
}
inline scalar pengRobinson::Tcrit() const
{
return Tcrit_;
}
inline void pengRobinson::updateModelCoefficients(const scalar T) const
{
aSave=a0_*pow(1 + n_*(1 - pow(T/Tcrit_, 0.5)), 2);
daSave=a0_*n_*(n_*sqrt(T/Tcrit_) - n_ - 1)*sqrt(T/Tcrit_)/T;
d2aSave=a0_*n_*(n_+1)*sqrt(T/Tcrit_)/(2*pow(T,2));
d2aSave=a0_*n_*(n_ + 1)*sqrt(T/Tcrit_)/(2*T*T);
//CL: saving the temperature at which the coefficients are valid
TSave=T;
}
inline scalar pengRobinson::rhostd()const
{
return rhostd_;
}
//CL: Model coefficient a(T)
inline scalar pengRobinson::a(const scalar T) const
{
@ -170,33 +194,13 @@ inline scalar pengRobinson::n()const
return n_;
}
//CL: pow of constant b() used in the code e.g. b2_=b*b;
inline scalar pengRobinson::b2()const
{
return b2_;
}
inline scalar pengRobinson::b3()const
{
return b3_;
}
inline scalar pengRobinson::b4()const
{
return b4_;
}
inline scalar pengRobinson::b5()const
{
return b5_;
}
//returns the pressure for a given density and temperature
inline scalar pengRobinson::p(const scalar rho, const scalar T) const
{
scalar Vm = this->W()/rho;
return this->RR()*T/(Vm-b())-a(T)/(pow(Vm,2)+2*b()*Vm-b2());
scalar Vm2 = Vm*Vm;
return this->RR()*T/(Vm - b_) - a(T)/(Vm2 + 2*b_*Vm - b2_);
}
@ -206,20 +210,21 @@ inline scalar pengRobinson::dpdv(const scalar rho, const scalar T) const
{
scalar Vm = this->W()/rho;
scalar Vm2 = Vm*Vm;
scalar Vm3 = Vm*Vm*Vm;
scalar Vm3 = Vm2*Vm;
scalar Vm4 = Vm3*Vm;
return(
return
(
2*a(T)*
(
b3()-b2()*Vm-b()*Vm2+Vm3
b3_ - b2_*Vm - b_*Vm2 + Vm3
)
- this->RR()*T*
(
b4()-4*b3()*Vm+2*b2()*Vm2
+4*b()*Vm3+pow(Vm,4)
b4_ - 4*b3_*Vm + 2*b2_*Vm2 + 4*b_*Vm3 + Vm4
)
)
/(pow(b()-Vm,2)*pow(b2()-2*b()*Vm-Vm2,2));
/(pow(b_ - Vm, 2)*pow(b2_ - 2*b_*Vm - Vm2, 2));
}
@ -228,8 +233,8 @@ inline scalar pengRobinson::dpdv(const scalar rho, const scalar T) const
inline scalar pengRobinson::dpdT(const scalar rho, const scalar T) const
{
scalar Vm = this->W()/rho;
return this->RR()/(Vm-b())-dadT(T)/(pow(Vm,2)+2*b()*Vm-b2());
scalar Vm2 = Vm*Vm;
return this->RR()/(Vm - b_) - dadT(T)/(Vm2 + 2*b_*Vm - b2_);
}
@ -238,10 +243,11 @@ inline scalar pengRobinson::dpdT(const scalar rho, const scalar T) const
//(molar values)
inline scalar pengRobinson::dvdT(const scalar rho, const scalar T) const
{
return (-1)*this->dpdT(rho,T)/this->dpdv(rho,T);
return -this->dpdT(rho, T)/this->dpdv(rho, T);
}
//Real deviative dv/dp at constant temperature
//(molar values)
inline scalar pengRobinson::dvdp(const scalar rho, const scalar T) const
{
@ -249,25 +255,38 @@ inline scalar pengRobinson::dvdp(const scalar rho,const scalar T) const
}
//(molar values)
//needed to calculate the internal energy
inline scalar pengRobinson::integral_p_dv(const scalar rho,const scalar T) const
//(molar values)
inline scalar pengRobinson::integral_p_dv
(
const scalar rho,
const scalar T
) const
{
scalar Vm = this->W()/rho;
scalar root2=pow(2, 0.5);
return -pow(2,0.5)*a(T)*log(b()*(1-pow(2,0.5))+Vm)/(4*b())+this->RR()*T*log(Vm-b())
+pow(2,0.5)*a(T)*log(b()*(pow(2,0.5)+1)+Vm)/(4*b());
return
- root2*a(T)*log(b_*(1-root2) + Vm)/(4*b_)
+ this->RR()*T*log(Vm - b_)
+ root2*a(T)*log(b_*(root2 + 1) + Vm)/(4*b_);
}
//needed to calculate the entropy
//(molar values)
//needed to calculate the entropy
inline scalar pengRobinson::integral_dpdT_dv(const scalar rho,const scalar T) const
inline scalar pengRobinson::integral_dpdT_dv
(
const scalar rho,
const scalar T
) const
{
scalar Vm = this->W()/rho;
return -pow(2,0.5)*dadT(T)*log(b()*(1-pow(2,0.5))+Vm)/(4*b())
+this->RR()*log(Vm-b())+pow(2,0.5)*dadT(T)*log(b()*(pow(2,0.5)+1)+Vm)/(4*b());
scalar root2=pow(2, 0.5);
return
- root2*dadT(T)*log(b_*(1 - root2) + Vm)/(4*b_)
+ this->RR()*log(Vm - b_) + root2*dadT(T)*log(b_*(root2 + 1) + Vm)/(4*b_);
}
@ -275,8 +294,8 @@ inline scalar pengRobinson::integral_dpdT_dv(const scalar rho,const scalar T) co
inline scalar pengRobinson::d2pdT2(const scalar rho, const scalar T) const
{
scalar Vm = this->W()/rho;
return -d2adT2(T)/(pow(Vm,2)+2*b()*Vm-b2());
scalar Vm2 = Vm*Vm;
return -d2adT2(T)/(Vm2 + 2*b_*Vm-b2_);
}
@ -285,75 +304,95 @@ inline scalar pengRobinson::d2pdv2(const scalar rho,const scalar T) const
{
scalar Vm = this->W()/rho;
scalar Vm2 = Vm*Vm;
scalar Vm3 = Vm*Vm*Vm;
scalar Vm4 = Vm*Vm*Vm*Vm;
scalar Vm5 = Vm*Vm*Vm*Vm*Vm;
scalar Vm3 = Vm2*Vm;
scalar Vm4 = Vm3*Vm;
scalar Vm5 = Vm4*Vm;
scalar Vm6 = Vm5*Vm;
return 2*
(
a(T)*
(
5*b5()-9*b4()*Vm+4*b2()*Vm3+3*b()*Vm4-3*Vm5
5*b5_ - 9*b4_*Vm + 4*b2_*Vm3 + 3*b_*Vm4 - 3*Vm5
)
- this->RR()*T*
(
pow(b(),6)-6*b5()*Vm+9*b4()*Vm2+4*b3()*Vm3
-9*b2()*Vm4-6*b()*Vm5-pow(Vm,6)
b6_ - 6*b5_*Vm + 9*b4_*Vm2 + 4*b3_*Vm3 - 9*b2_*Vm4 - 6*b_*Vm5 - Vm6
)
)
/(pow(b()-Vm,3)*pow(b2()-2*b()*Vm-Vm2,3));
/(pow(b_ - Vm, 3)*pow(b2_ - 2*b_*Vm - Vm2, 3));
}
//(molar values)
//using second order implicit differentiation
inline scalar pengRobinson::d2vdT2(const scalar rho, const scalar T) const
inline scalar pengRobinson::d2vdT2
(
const scalar rho,
const scalar T
) const
{
scalar dpdT2=this->dpdT(rho, T)*this->dpdT(rho, T);
scalar dpdv2=this->dpdv(rho, T)*this->dpdv(rho, T);
scalar dpdv3=dpdv2*this->dpdv(rho, T);
return
-(
pow(this->dpdT(rho,T),2)*this->d2pdv2(rho,T)
+ pow(this->dpdv(rho,T),2)*this->d2pdT2(rho,T)
dpdT2*this->d2pdv2(rho, T)
+ dpdv2*this->d2pdT2(rho, T)
- 2*this->dpdv(rho, T)*this->dpdT(rho, T)*this->d2pdvdT(rho, T)
)
/(pow(this->dpdv(rho,T),3));
/dpdv3;
}
//(molar values)
inline scalar pengRobinson::d2pdvdT(const scalar rho, const scalar T) const
inline scalar pengRobinson::d2pdvdT
(
const scalar rho,
const scalar T
) const
{
scalar Vm = this->W()/rho;
scalar Vm2 = Vm*Vm;
scalar Vm3 = Vm*Vm*Vm;
scalar Vm3 = Vm2*Vm;
scalar Vm4 = Vm3*Vm;
return(
2*dadT(T)*
return
(
b3()-b2()*Vm-b()*Vm2+Vm3
2*dadT(T)*(b3_ - b2_*Vm - b_*Vm2 + Vm3)
- this->RR()*(b4_ - 4*b3_*Vm + 2*b2_*Vm2 + 4*b_*Vm3 + Vm4)
)
-this->RR()*
(
b4()-4*b3()*Vm+2*b2()*Vm2
+4*b()*Vm3+pow(Vm,4)
)
)
/(pow(b()-Vm,2)*pow(b2()-2*b()*Vm-Vm2,2));
/(pow(b_ - Vm, 2)*pow(b2_ - 2*b_*Vm - Vm2, 2));
}
// the result of this intergal is needed for the nasa based cp polynomial
//(molar values)
inline scalar pengRobinson::integral_d2pdT2_dv(const scalar rho,const scalar T) const
inline scalar pengRobinson::integral_d2pdT2_dv
(
const scalar rho,
const scalar T
) const
{
scalar Vm = this->W()/rho;
scalar root2=pow(2, 0.5);
return pow(2,0.5)*d2adT2(T)*log(b()*(pow(2,0.5)+1)+Vm)/(4*b())-pow(2,0.5)*d2adT2(T)*log(b()*(1-pow(2,0.5))+Vm)/(4*b());
return root2*d2adT2(T)/(4*b_)
*(log(b_*(root2 + 1) + Vm) - log(b_*(1 - root2) + Vm));
// root2*d2adT2(T)*log(b_*(root2 + 1) + Vm)/(4*b_)
//- root2*d2adT2(T)*log(b_*(1 - root2)+Vm)/(4*b_);
}
//Isobar expansion Coefficent beta = 1/v (dv/dt) at constant p
//(molar values)
inline scalar pengRobinson::isobarExpCoef(const scalar rho,const scalar T) const
inline scalar pengRobinson::isobarExpCoef
(
const scalar rho,
const scalar T
) const
{
return this->dvdT(rho, T)*rho/this->W();
}
@ -361,7 +400,11 @@ inline scalar pengRobinson::isobarExpCoef(const scalar rho,const scalar T) cons
//isothemal compressiblity kappa (not Thermal conductivity)
//(molar values)
inline scalar pengRobinson::isothermalCompressiblity(const scalar rho,const scalar T) const
inline scalar pengRobinson::isothermalCompressiblity
(
const scalar rho,
const scalar T
) const
{
return this->isobarExpCoef(rho, T)/this->dpdT(rho, T);
//CL: also possible
@ -370,13 +413,13 @@ inline scalar pengRobinson::isothermalCompressiblity(const scalar rho,const scal
//- Return density [kg/m^3]
inline scalar pengRobinson::rho(
inline scalar pengRobinson::rho
(
const scalar p,
const scalar T,
const scalar rho0
) const
{
scalar molarVolumePrevIteration;
scalar molarVolume;
label iter = 0;
@ -394,6 +437,7 @@ inline scalar pengRobinson::rho(
label i=0;
do
{
//CL: modified Newton solver
molarVolume=molarVolumePrevIteration
-(
(this->p((this->W()/molarVolumePrevIteration),T) - p)
@ -404,8 +448,7 @@ inline scalar pengRobinson::rho(
if (i>8)
{
//CL: using bisection methode as backup,
//CL: solution must be between rho=0.001 to rho=1500;
//CL: if not, change rhoMax_ and rhoMin_
//CL: solution must be between rhoMin_ to rhoMax
for(i=0; i<200; i++)
{
scalar f1 = this->p(rho1, T) - p;
@ -449,7 +492,8 @@ inline scalar pengRobinson::rho(
{
FatalErrorIn
(
"inline scalar pengRobinson::rho(const scalar p, const scalar T, const scalar rho0) const "
"inline scalar pengRobinson::rho"
"(const scalar p, const scalar T, const scalar rho0) const "
) << "Maximum number of iterations exceeded"
<< abort(FatalError);
}
@ -463,7 +507,7 @@ inline scalar pengRobinson::rho(
//- Return density [kg/m^3]on
inline scalar pengRobinson::rho(const scalar p, const scalar T) const
{
// using perfect gas equation as starting point
//CL: using perfect gas equation as starting point
return rho(p, T, p/(this->R()*T));
}
@ -471,12 +515,17 @@ inline scalar pengRobinson::rho(const scalar p,const scalar T) const
//- Return compressibility drho/dp at T=constant [s^2/m^2]
inline scalar pengRobinson::psi(const scalar rho, const scalar T) const
{
return -this->dvdp(rho,T)*pow(rho,2)/this->W();
return -this->dvdp(rho, T)*rho*rho/this->W();
}
//- Return compression factor []
inline scalar pengRobinson::Z( const scalar p, const scalar T,const scalar rho0) const
inline scalar pengRobinson::Z
(
const scalar p,
const scalar T,
const scalar rho0
) const
{
return p/(this->R()*T*this->rho(p, T, rho0));
}
@ -484,24 +533,18 @@ inline scalar pengRobinson::Z( const scalar p, const scalar T,const scalar rho0)
// * * * * * * * * * * * * * * * Member Operators * * * * * * * * * * * * * //
inline void pengRobinson::operator+=(const pengRobinson& pr)
{
specie::operator+=(pr);
}
/*
inline void pengRobinson::operator-=(const pengRobinson& pr)
{
specie::operator-=(pr);
}
inline void pengRobinson::operator*=(const scalar s)
{
specie::operator*=(s);
}
*/
// * * * * * * * * * * * * * * * Friend Operators * * * * * * * * * * * * * //
inline pengRobinson operator+
@ -517,6 +560,7 @@ inline pengRobinson operator+
);
}
inline pengRobinson operator*
(
const scalar s,
@ -526,6 +570,7 @@ inline pengRobinson operator*
return pengRobinson(s*static_cast<const specie&>(pr));
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam

37
src/thermophysicalModels/specie/equationOfState/cubicEquationOfState/redlichKwong/redlichKwong.C
View file

@ -36,53 +36,52 @@ Germany
#include "redlichKwong.H"
#include "IOstreams.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
redlichKwong::redlichKwong(Istream& is)
Foam::redlichKwong::redlichKwong(Istream& is)
:
specie(is),
pcrit_(readScalar(is)),
Tcrit_(readScalar(is)),
a_(0.42748*pow(this->RR(), 2)*pow(Tcrit_, 2.5)/pcrit_),
b_(0.08664*this->RR()*Tcrit_/pcrit_),
//CL: Only uses the default values
b2_(pow(b_,2)),
b3_(pow(b_,3)),
b5_(pow(b_,5)),
rhoMax_(1500),
//CL: Only uses the default values
rhoMin_(1e-3),
rhoMax_(1500),
// Starting GUESS for the density by ideal gas law
rhostd_(this->rho(this->Pstd(), this->Tstd(), this->Pstd()/(this->Tstd()*this->R())))
{
is.check("redlichKwong::redlichKwong(Istream& is)");
}
//CL: Constructed needed in OpenFOAM 2.x.x
//CL: Code works fine, but compiling problem in OpenFOAM 1.6.ext
//CL: because specie has no constructor using dict
/*
redlichKwong::redlichKwong(const dictionary& dict)
Foam::redlichKwong::redlichKwong(const dictionary& dict)
:
specie(dict),
pcrit_(readScalar(dict.subDict("equationOfState").lookup("pCritical"))),
Tcrit_(readScalar(dict.subDict("equationOfState").lookup("TCritical"))),
//CL: rhoMin and rhoMax are only used as boundaries for the bisection methode (see rho function)
//CL: important: rhoMin and rhoMax are not used as boundary for the newton solver
//CL: therefore, rho can be larger than rhoMax and smaller than rhoMin
rhoMin_(dict.subDict("equationOfState").lookupOrDefault("rhoMin",1e-3)),
rhoMax_(dict.subDict("equationOfState").lookupOrDefault("rhoMax",1500)),
a_(0.42748*pow(this->RR(),2)*pow(Tcrit_,2.5)/pcrit_),
b_(0.08664*this->RR()*Tcrit_/pcrit_),
b2_(pow(b_,2)),
b3_(pow(b_,3)),
b5_(pow(b_,5)),
//CL: rhoMin and rhoMax are only used as boundaries for the bisection methode (see rho function)
//CL: important: rhoMin and rhoMax are not used as boundary for the newton solver
//CL: therefore, rho can be larger than rhoMax and smaller than rhoMin
rhoMin_(dict.subDict("equationOfState").lookupOrDefault("rhoMin",1e-3)),
rhoMax_(dict.subDict("equationOfState").lookupOrDefault("rhoMax",1500)),
// Starting GUESS for the density by ideal gas law
rhostd_(this->rho(Pstd, Tstd, Pstd/(Tstd*this->R())))
{}
rhostd_(this->rho(this->Pstd(), this->Tstd(), this->Pstd()/(this->Tstd()*this->R())))
{
is.check("redlichKwong::redlichKwong(Istream& is)");
}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
@ -100,9 +99,10 @@ void Foam::redlichKwong::write(Ostream& os) const
}
*/
// * * * * * * * * * * * * * * * Ostream Operator * * * * * * * * * * * * * //
Ostream& operator<<(Ostream& os, const redlichKwong& rk)
Foam::Ostream& Foam::operator<<(Ostream& os, const redlichKwong& rk)
{
os << static_cast<const specie&>(rk)<< token::SPACE
<< rk.pcrit_ << tab<< rk.Tcrit_;
@ -111,8 +111,5 @@ Ostream& operator<<(Ostream& os, const redlichKwong& rk)
return os;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// ************************************************************************* //

54
src/thermophysicalModels/specie/equationOfState/cubicEquationOfState/redlichKwong/redlichKwong.H
View file

@ -59,29 +59,26 @@ class redlichKwong
:
public specie
{
protected:
// Protected data
// private data
//CL: data at critical point
scalar pcrit_;
scalar Tcrit_;
//-Redlich Kwong factors
mutable scalar a_;
mutable scalar b_;
//CL: Redlich Kwong factors
scalar a_;
scalar b_;
//CL: pow of constants b_ used in the code e.g. b2_=b*b;
mutable scalar b2_;
mutable scalar b3_;
mutable scalar b5_;
scalar b2_;
scalar b3_;
scalar b5_;
//CL: rhoMin and rhoMax are only used as boundaries for the bisection methode (see rho function)
scalar rhoMax_;
scalar rhoMin_;
scalar rhoMax_;
//- Density @STD, initialise after a, b!
mutable scalar rhostd_;
scalar rhostd_;
public:
@ -110,24 +107,23 @@ public:
// Member functions
inline scalar rhostd() const;
//Return Redlich Kwong factors
inline scalar a() const;
inline scalar b() const;
inline scalar rhostd() const;
inline scalar rhoMin() const;
inline scalar rhoMax() const;
inline scalar Tcrit() const;
//CL: Equation of state
inline scalar p(const scalar rho, const scalar T) const;
//CL: return power of constants b_
inline scalar b2()const;
inline scalar b3()const;
inline scalar b5()const;
//first order derivatives
//CL: first order derivatives
inline scalar dpdv(const scalar rho, const scalar T) const;
inline scalar dpdT(const scalar rho, const scalar T) const;
@ -144,7 +140,7 @@ public:
const scalar T
) const;
// Used for cv
//CL: Used for cv
inline scalar integral_d2pdT2_dv
(
const scalar rho,
@ -157,7 +153,7 @@ public:
// Used for Entropy
inline scalar integral_dpdT_dv(const scalar rho, const scalar T) const;
// second order derivatives, not Used At The Moment
//CL: second order derivatives, not Used At The Moment
inline scalar d2pdv2(const scalar rho, const scalar T) const;
inline scalar d2pdT2(const scalar rho, const scalar T) const;
@ -198,10 +194,9 @@ public:
// Member operators
inline void operator+=(const redlichKwong&);
/*
inline void operator-=(const redlichKwong&);
inline void operator*=(const scalar);
*/
// Friend operators
@ -217,6 +212,7 @@ public:
const redlichKwong&
);
// Ostream Operator
friend Ostream& operator<<(Ostream&, const redlichKwong&);

121
src/thermophysicalModels/specie/equationOfState/cubicEquationOfState/redlichKwong/redlichKwongI.H
View file

@ -60,6 +60,8 @@ inline redlichKwong::redlichKwong(const word& name, const redlichKwong& rk)
b2_(rk.b2_),
b3_(rk.b3_),
b5_(rk.b5_),
rhoMin_(rk.rhoMin_),
rhoMax_(rk.rhoMax_),
rhostd_(rk.rhostd_)
{}
@ -86,6 +88,24 @@ inline scalar redlichKwong::rhostd()const
}
inline scalar redlichKwong::rhoMin() const
{
return rhoMin_;
}
inline scalar redlichKwong::rhoMax() const
{
return rhoMax_;
}
inline scalar redlichKwong::Tcrit() const
{
return Tcrit_;
}
inline scalar redlichKwong::a() const
{
return a_;
@ -97,22 +117,6 @@ inline scalar redlichKwong::b()const
return b_;
}
inline scalar redlichKwong::b2()const
{
return b2_;
}
inline scalar redlichKwong::b3()const
{
return b3_;
}
inline scalar redlichKwong::b5()const
{
return b5_;
}
//returns the pressure for a given density and temperature
inline scalar redlichKwong::p(const scalar rho, const scalar T) const
@ -128,9 +132,12 @@ inline scalar redlichKwong::dpdv(const scalar rho, const scalar T) const
{
scalar Vm = this->W()/rho;
scalar Vm2 = Vm*Vm;
scalar Vm3 = Vm2*Vm;
return (a_*(b3() - 3*b_*Vm2 + 2*pow(Vm,3))
- this->RR()*pow(T,1.5)*Vm2*(b2() + 2*b_*Vm + Vm2))
return
(
a_*(b3_ - 3*b_*Vm2 + 2*Vm3) - this->RR()*pow(T, 1.5)*Vm2*(b2_ + 2*b_*Vm + Vm2)
)
/(sqrt(T)*Vm2*pow((b_ + Vm), 2)*pow( (b_ - Vm), 2));
}
@ -170,9 +177,8 @@ inline scalar redlichKwong::integral_p_dv
) const
{
scalar Vm = this->W()/rho;
return this->RR()*T*log(Vm - b_)
+ (a_*log(b_ + Vm))/(b_*sqrt(T))
- (a_*log(Vm))/(b_*sqrt(T));
return this->RR()*T*log(Vm - b_) + a_/(b_*sqrt(T))*(log(b_ + Vm) - log(Vm));
//return this->RR()*T*log(Vm - b_) + (a_*log(b_ + Vm))/(b_*sqrt(T)) - (a_*log(Vm))/(b_*sqrt(T));
}
@ -185,9 +191,8 @@ inline scalar redlichKwong::integral_dpdT_dv
) const
{
scalar Vm = this->W()/rho;
return this->RR()*log(Vm - b_)
-(a_*log(b_ + Vm))/(2*b_*pow(T,1.5))
+(a_*log(Vm))/(2*b_*pow(T,1.5));
return this->RR()*log(Vm - b_) - a_/(2*b_*pow(T, 1.5))*(log(b_ + Vm) - log(Vm));
//return this->RR()*log(Vm - b_) - (a_*log(b_ + Vm))/(2*b_*T15_) + (a_*log(Vm))/(2*b_*T15_);
}
@ -204,19 +209,16 @@ inline scalar redlichKwong::d2pdv2(const scalar rho, const scalar T) const
{
scalar Vm = this->W()/rho;
scalar Vm2 = Vm*Vm;
scalar Vm3 = Vm*Vm*Vm;
scalar Vm3 = Vm2*Vm;
scalar Vm4 = Vm3*Vm;
scalar Vm5 = Vm4*Vm;
return
(
2*(
a_*(
b5()-3*b3()*Vm2
- b2()*Vm3
+ 6*b_*pow(Vm,4)-3*pow(Vm,5)
)
+ this->RR()*pow(T,1.5)*Vm3*(b3()
+ 3*b2()*Vm
+ 3*b_*Vm2+Vm3)
2*
(
a_*(b5_ - 3*b3_*Vm2 - b2_*Vm3 + 6*b_*Vm4 - 3*Vm5)
+ this->RR()*pow(T, 1.5)*Vm3*(b3_ + 3*b2_*Vm + 3*b_*Vm2 + Vm3)
)
/(sqrt(T)*Vm3*pow((b_ + Vm), 3)*pow(Vm - b_, 3))
);
@ -224,16 +226,20 @@ inline scalar redlichKwong::d2pdv2(const scalar rho, const scalar T) const
//(molar values)
//using second Order implicit differentiation
//using second order implicit differentiation
inline scalar redlichKwong::d2vdT2(const scalar rho, const scalar T) const
{
scalar dpdT2=this->dpdT(rho, T)*this->dpdT(rho, T);
scalar dpdv2=this->dpdv(rho, T)*this->dpdv(rho, T);
scalar dpdv3=dpdv2*this->dpdv(rho, T);
return
-(
pow(this->dpdT(rho,T),2)*this->d2pdv2(rho,T)
+ pow(this->dpdv(rho,T),2)*this->d2pdT2(rho,T)
dpdT2*this->d2pdv2(rho, T)
+ dpdv2*this->d2pdT2(rho, T)
- 2*this->dpdv(rho, T)*this->dpdT(rho, T)*this->d2pdvdT(rho, T)
)
/(pow(this->dpdv(rho,T),3));
/dpdv3;
}
@ -242,13 +248,18 @@ inline scalar redlichKwong::d2pdvdT(const scalar rho, const scalar T) const
{
scalar Vm = this->W()/rho;
scalar Vm2 = Vm*Vm;
scalar Vm3 = Vm2*Vm;
scalar T15_ = pow(T, 1.5);
return
-(0.5*(
a_*(b3() - 3*b_*Vm2 + 2*pow(Vm,3))
+ 2*this->RR()*pow(T,1.5)*Vm2*(b2() + 2*b_*Vm + Vm2)
))
/(pow(T,1.5)*Vm2*pow(b_ + Vm,2)*pow(b_ - Vm,2));
-(
0.5*
(
a_*(b3_ - 3*b_*Vm2 + 2*Vm3)
+ 2*this->RR()*T15_*Vm2*(b2_ + 2*b_*Vm + Vm2)
)
)
/(T15_*Vm2*pow(b_ + Vm, 2)*pow(b_ - Vm, 2));
}
@ -260,9 +271,9 @@ inline scalar redlichKwong::integral_d2pdT2_dv
const scalar T
) const
{
scalar T25_=pow(T,2.5);
scalar Vm = this->W()/rho;
return 0.75*a_*log(b_ + Vm)/(pow(T,2.5)*b_)
- 0.75*a_*log(Vm)/(pow(T,2.5)*b_);
return 0.75*a_*log(b_ + Vm)/(T25_*b_) - 0.75*a_*log(Vm)/(T25_*b_);
}
@ -299,7 +310,6 @@ inline scalar redlichKwong::rho
const scalar rho0
) const
{
scalar molarVolumePrevIteration;
scalar molarVolume;
label iter = 0;
@ -372,7 +382,8 @@ inline scalar redlichKwong::rho
{
FatalErrorIn
(
"inline scalar redlichKwong::rho(const scalar p, const scalar T, const scalar rho0) const "
"inline scalar redlichKwong::rho"
"(const scalar p, const scalar T, const scalar rho0) const "
) << "Maximum number of iterations exceeded"
<< abort(FatalError);
}
@ -383,7 +394,7 @@ inline scalar redlichKwong::rho
}
//- Return density [kg/m^3]on
//- Return density [kg/m^3]
inline scalar redlichKwong::rho(const scalar p, const scalar T) const
{
// using perfect gas equation as starting point
@ -394,7 +405,7 @@ inline scalar redlichKwong::rho(const scalar p, const scalar T) const
//- Return compressibility drho/dp at T=constant [s^2/m^2]
inline scalar redlichKwong::psi(const scalar rho, const scalar T) const
{
return -this->dvdp(rho,T)*pow(rho,2)/this->W();
return -this->dvdp(rho, T)*rho*rho/this->W();
}
@ -412,27 +423,19 @@ inline scalar redlichKwong::Z
// * * * * * * * * * * * * * * * Member Operators * * * * * * * * * * * * * //
inline void redlichKwong::operator+=(const redlichKwong& rk)
{
specie::operator+=(rk);
}
/*
inline void redlichKwong::operator-=(const redlichKwong& rk)
{
specie::operator-=(rk);
}
inline void redlichKwong::operator*=(const scalar s)
{
specie::operator*=(s);
}
*/
// * * * * * * * * * * * * * * * Friend Operators * * * * * * * * * * * * * //
// * * * * * * * * * * * * * * * Friend Operators * * * * * * * * * * * * * //
inline redlichKwong operator+
(
@ -447,6 +450,7 @@ inline redlichKwong operator+
);
}
inline redlichKwong operator*
(
const scalar s,
@ -456,6 +460,7 @@ inline redlichKwong operator*
return redlichKwong(s*static_cast<const specie&>(rk));
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam

66
src/thermophysicalModels/specie/equationOfState/cubicEquationOfState/soaveRedlichKwong/soaveRedlichKwong.C
View file

@ -36,34 +36,33 @@ Germany
#include "soaveRedlichKwong.H"
#include "IOstreams.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
soaveRedlichKwong::soaveRedlichKwong(Istream& is)
Foam::soaveRedlichKwong::soaveRedlichKwong(Istream& is)
:
specie(is),
pcrit_(readScalar(is)),
Tcrit_(readScalar(is)),
azentricFactor_(readScalar(is)),
n_(0.48508+1.55171*azentricFactor_-0.15613*pow(azentricFactor_,2)),
a0_(0.42747*pow(this->RR(),2)*pow(Tcrit_,2)/(pcrit_)),
a0_(0.42747*pow(this->RR(), 2)*pow(Tcrit_, 2)/pcrit_),
b_(0.08664*this->RR()*Tcrit_/pcrit_),
n_(0.48508 + 1.55171*azentricFactor_ - 0.15613*pow(azentricFactor_, 2)),
b2_(b_*b_),
b3_(b2_*b_),
b5_(b2_*b3_),
//CL: Only uses the default values
b2_(pow(b_,2)),
b3_(pow(b_,3)),
b5_(pow(b_,5)),
rhoMax_(1500),
rhoMin_(1e-3),
TSave(0.0),
rhoMax_(1500),
// Starting GUESS for the density by ideal gas law
rhostd_(this->rho(this->Pstd(),this->Tstd(),this->Pstd()/(this->Tstd()*this->R())))
rhostd_(this->rho(this->Pstd(), this->Tstd(), this->Pstd()/(this->Tstd()*this->R()))),
aSave(0.0),
daSave(0.0),
d2aSave(0.0),
TSave(0.0)
{
is.check("soaveRedlichKwong::soaveRedlichKwong(Istream& is)");
}
//CL: Constructed needed in OpenFOAM 2.x.x
//CL: Code works fine, but compiling problem in OpenFOAM 1.6.ext
//CL: because specie has no constructor using dict
@ -74,21 +73,32 @@ soaveRedlichKwong::soaveRedlichKwong(const dictionary& dict)
pcrit_(readScalar(dict.subDict("equationOfState").lookup("pCritical"))),
Tcrit_(readScalar(dict.subDict("equationOfState").lookup("TCritical"))),
azentricFactor_(readScalar(dict.subDict("equationOfState").lookup("azentricFactor"))),
//CL: rhoMin and rhoMax are only used as boundaries for the bisection methode (see rho function)
a0_(0.42747*pow(this->RR(), 2)*pow(Tcrit_, 2)/pcrit_),
b_(0.08664*this->RR()*Tcrit_/pcrit_),
n_(0.48508 + 1.55171*azentricFactor_ - 0.15613*pow(azentricFactor_, 2)),
b2_(b_*b_),
b3_(b_*b2_),
b5_(b3_*b2_),
//CL: rhoMin and rhoMax are only used as boundaries for the bisection method (see rho function)
//CL: important: rhoMin and rhoMax are not used as boundary for the newton solver
//CL: therefore, rho can be larger than rhoMax and smaller than rhoMin
rhoMin_(dict.subDict("equationOfState").lookupOrDefault("rhoMin",1e-3)),
rhoMax_(dict.subDict("equationOfState").lookupOrDefault("rhoMax",1500)),
//CL: rhoMin and rhoMax are only used as boundaries for the bisection method (see rho function)
//CL: important: rhoMin and rhoMax are not used as boundary for the newton solver
//CL: therefore, rho can be larger than rhoMax and smaller than rhoMin
rhoMin_(dict.subDict("equationOfState").lookupOrDefault("rhoMin",1e-3)),
rhoMax_(dict.subDict("equationOfState").lookupOrDefault("rhoMax",1500)),
a0_(0.42747*pow(this->RR(),2)*pow(Tcrit_,2)/(pcrit_)),
b_(0.08664*this->RR()*Tcrit_/pcrit_),
n_(0.48508+1.55171*azentricFactor_-0.15613*pow(azentricFactor_,2)),
TSave(0.0),
b2_(pow(b_,2)),
b3_(pow(b_,3)),
b5_(pow(b_,5)),
// Starting GUESS for the density by ideal gas law
rhostd_(this->rho(this->Pstd(),this->Tstd(),this->Pstd()/(this->Tstd()*this->R())))
{}
rhostd_(this->rho(this->Pstd(), this->Tstd(), this->Pstd()/(this->Tstd()*this->R()))),
aSave(0.0),
daSave(0.0),
d2aSave(0.0),
TSave(0.0)
{
is.check("soaveRedlichKwong::soaveRedlichKwong(Istream& is)");
}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
@ -106,9 +116,11 @@ void Foam::soaveRedlichKwong::write(Ostream& os) const
os << indent << dict.dictName() << dict;
}
*/
// * * * * * * * * * * * * * * * Ostream Operator * * * * * * * * * * * * * //
Ostream& operator<<(Ostream& os, const soaveRedlichKwong& srk)
Foam::Ostream& Foam::operator<<(Ostream& os, const soaveRedlichKwong& srk)
{
os << static_cast<const specie&>(srk)<< token::SPACE
<< srk.pcrit_ << tab<< srk.Tcrit_<<tab<<srk.azentricFactor_;
@ -118,8 +130,4 @@ Ostream& operator<<(Ostream& os, const soaveRedlichKwong& srk)
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// ************************************************************************* //

89
src/thermophysicalModels/specie/equationOfState/cubicEquationOfState/soaveRedlichKwong/soaveRedlichKwong.H
View file

@ -51,7 +51,6 @@ Germany
#include "specie.H"
#include "autoPtr.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
@ -64,28 +63,29 @@ namespace Foam
class soaveRedlichKwong
:
public specie
{
protected:
// Protected data
// private data
//CL: data at critical point
scalar pcrit_;
scalar Tcrit_;
scalar azentricFactor_;
//-Soave Redlich Kwong
scalar n_;
scalar a0_;
mutable scalar b_;
scalar b_;
scalar n_;
//CL: pow of constants b_ used in the code e.g. b2_=b*b;
mutable scalar b2_;
mutable scalar b3_;
mutable scalar b5_;
scalar b2_;
scalar b3_;
scalar b5_;
//CL: rhoMin and rhoMax are only used as boundaries for the bisection methode (see rho function)
scalar rhoMax_;
scalar rhoMin_;
scalar rhoMax_;
//- Density @STD, initialise after a0, b!
scalar rhostd_;
//CL: Variables to save the values of a, dadT and d2adT2 of the mixture
//CL: Variables must corrected for changing temperatures
@ -96,23 +96,18 @@ protected:
//CL: save the temperature for which the save coefficients (amix,dadTmix,d2adT2mix) are correct
mutable scalar TSave;
//- Density @STD, initialise after a0, b!
mutable scalar rhostd_;
//Protected functions
//CL: function updates the coefficients (aSave, daSave, d2aSave)
inline void updateModelCoefficients(const scalar T) const;
public:
// Constructors
//- Construct from components
inline soaveRedlichKwong
(
const specie& s
const specie& sp
);
//- Construct from Istream
@ -131,8 +126,22 @@ public:
inline static autoPtr<soaveRedlichKwong> New(Istream& is);
// Member functions
//Return Soave Redlich Kwong factors
inline scalar a0() const;
inline scalar b() const;
inline scalar n() const;
inline scalar rhostd() const;
inline scalar rhoMin() const;
inline scalar rhoMax() const;
inline scalar Tcrit() const;
//CL: Model coefficient a(T)
inline scalar a(const scalar T) const;
@ -142,20 +151,6 @@ public:
//CL: second order temperature deriviative of model coefficient a(T)
inline scalar d2adT2(const scalar T) const;
//Return Soave Redlich Kwong factors
inline scalar a0()const;
inline scalar b()const;
inline scalar n()const;
//CL: return power of constants b_
inline scalar b2()const;
inline scalar b3()const;
inline scalar b5()const;
//CL: Equation of state
inline scalar p(const scalar rho, const scalar T) const;
@ -168,11 +163,7 @@ public:
inline scalar dvdp(const scalar rho, const scalar T) const;
inline scalar isobarExpCoef
(
const scalar rho,
const scalar T
) const;
inline scalar isobarExpCoef(const scalar rho, const scalar T) const;
inline scalar isothermalCompressiblity
(
@ -181,7 +172,11 @@ public:
) const;
//CL: Used for cv
inline scalar integral_d2pdT2_dv(const scalar rho,const scalar T) const ;
inline scalar integral_d2pdT2_dv
(
const scalar rho,
const scalar T
) const;
//CL: second order derivatives, not Used At The Moment
inline scalar d2pdv2(const scalar rho, const scalar T) const;
@ -193,10 +188,18 @@ public:
inline scalar d2vdT2(const scalar rho, const scalar T) const;
//CL: Used for internal Energy
inline scalar integral_p_dv(const scalar rho,const scalar T) const;
inline scalar integral_p_dv
(
const scalar rho,
const scalar T
) const;
//Used for Entropy
inline scalar integral_dpdT_dv(const scalar rho,const scalar T) const;
//CL: Used for Entropy
inline scalar integral_dpdT_dv
(
const scalar rho,
const scalar T
) const;
//- Return density [kg/m^3]
// rho0 is the starting point of the newton solver used to calculate rho
@ -232,11 +235,8 @@ public:
inline void operator+=(const soaveRedlichKwong&);
/*
inline void operator-=(const soaveRedlichKwong&);
inline void operator*=(const scalar);
*/
// Friend operators
@ -252,6 +252,7 @@ public:
const soaveRedlichKwong&
);
// Ostream Operator
friend Ostream& operator<<(Ostream&, const soaveRedlichKwong&);

201
src/thermophysicalModels/specie/equationOfState/cubicEquationOfState/soaveRedlichKwong/soaveRedlichKwongI.H
View file

@ -28,7 +28,6 @@ Institut für Thermodynamik
Technische Universität Braunschweig
Germany
\*---------------------------------------------------------------------------*/
#include "soaveRedlichKwong.H"
@ -58,14 +57,19 @@ inline soaveRedlichKwong::soaveRedlichKwong(const word& name,const soaveRedlich
pcrit_(srk.pcrit_),
Tcrit_(srk.Tcrit_),
azentricFactor_(srk.azentricFactor_),
n_(srk.n_),
a0_(srk.a0_),
b_(srk.b_),
n_(srk.n_),
b2_(srk.b2_),
b3_(srk.b3_),
b5_(srk.b5_),
TSave(0),
rhostd_(srk.rhostd_)
rhoMin_(srk.rhoMin_),
rhoMax_(srk.rhoMax_),
rhostd_(srk.rhostd_),
aSave(0.0),
daSave(0.0),
d2aSave(0.0),
TSave(0.0)
{}
@ -85,23 +89,41 @@ inline autoPtr<soaveRedlichKwong> soaveRedlichKwong::New(Istream& is)
// * * * * * * * * * * * * * Member Functions * * * * * * * * * * * //
inline scalar soaveRedlichKwong::rhostd() const
{
return rhostd_;
}
inline scalar soaveRedlichKwong::rhoMin() const
{
return rhoMin_;
}
inline scalar soaveRedlichKwong::rhoMax() const
{
return rhoMax_;
}
inline scalar soaveRedlichKwong::Tcrit() const
{
return Tcrit_;
}
inline void soaveRedlichKwong::updateModelCoefficients(const scalar T) const
{
aSave=a0_*pow(1 + n_*(1 - pow(T/Tcrit_, 0.5)), 2);
daSave=a0_*n_*(n_*sqrt(T/Tcrit_) - n_ - 1)*sqrt(T/Tcrit_)/T;
d2aSave=a0_*n_*(n_+1)*sqrt(T/Tcrit_)/(2*pow(T,2));
d2aSave=a0_*n_*(n_ + 1)*sqrt(T/Tcrit_)/(2*T*T);
//CL: saving the temperature at which the coefficients are valid
TSave=T;
}
inline scalar soaveRedlichKwong::rhostd()const
{
return rhostd_;
}
//CL: Model coefficient a(T)
inline scalar soaveRedlichKwong::a(const scalar T) const
{
@ -170,33 +192,12 @@ inline scalar soaveRedlichKwong::n()const
return n_;
}
//CL: pow of constant b() used in the code e.g. b2_=b*b;
inline scalar soaveRedlichKwong::b2()const
{
return b2_;
}
inline scalar soaveRedlichKwong::b3()const
{
return b3_;
}
inline scalar soaveRedlichKwong::b5()const
{
return b5_;
}
//returns the pressure for a given density and temperature
inline scalar soaveRedlichKwong::p(const scalar rho, const scalar T) const
{
scalar Vm = this->W()/rho;
return
(
this->RR()*T/(Vm-b_)
-a(T)/(Vm*(Vm+b_))
);
return this->RR()*T/(Vm - b_) - a(T)/(Vm*(Vm + b_));
}
@ -206,11 +207,11 @@ inline scalar soaveRedlichKwong::dpdv(const scalar rho, const scalar T) const
{
scalar Vm = this->W()/rho;
scalar Vm2 = Vm*Vm;
scalar Vm3 = Vm2*Vm;
return
(
a(T)*(b3()-3*b_*Vm2+2*pow(Vm,3))
-this->RR()*T*Vm2*(b2()+2*b_*Vm+Vm2)
a(T)*(b3_ - 3*b_*Vm2 + 2*Vm3) - this->RR()*T*Vm2*(b2_ + 2*b_*Vm + Vm2)
)
/(Vm2*pow(b_ + Vm, 2)*pow(b_ - Vm, 2));
}
@ -221,12 +222,7 @@ inline scalar soaveRedlichKwong::dpdv(const scalar rho, const scalar T) const
inline scalar soaveRedlichKwong::dpdT(const scalar rho, const scalar T) const
{
scalar Vm = this->W()/rho;
return
(
this->RR()/(Vm-b_)
-dadT(T)/(Vm*(Vm+b_))
);
return this->RR()/(Vm - b_) - dadT(T)/(Vm*(Vm + b_));
}
@ -235,7 +231,7 @@ inline scalar soaveRedlichKwong::dpdT(const scalar rho, const scalar T) const
// (molar values)
inline scalar soaveRedlichKwong::dvdT(const scalar rho, const scalar T) const
{
return (-1)*this->dpdT(rho,T)/this->dpdv(rho,T);
return -this->dpdT(rho, T)/this->dpdv(rho, T);
}
@ -249,21 +245,30 @@ inline scalar soaveRedlichKwong::dvdp(const scalar rho,const scalar T) const
//needed to calculate the internal energy
//(molar values)
inline scalar soaveRedlichKwong::integral_p_dv(const scalar rho,const scalar T) const
inline scalar soaveRedlichKwong::integral_p_dv
(
const scalar rho,
const scalar T
) const
{
scalar Vm = this->W()/rho;
return this->RR()*T*log(Vm-b_)+a(T)*log(b_+Vm)/b_-a(T)*log(Vm)/b_;
return this->RR()*T*log(Vm - b_) + a(T)/b_*(log(b_ + Vm) - log(Vm));
//return this->RR()*T*log(Vm - b_) + a(T)*log(b_ + Vm)/b_ - a(T)*log(Vm)/b_;
}
//needed to calculate the entropy
//(molar values)
inline scalar soaveRedlichKwong::integral_dpdT_dv(const scalar rho,const scalar T) const
//needed to calculate the entropy
inline scalar soaveRedlichKwong::integral_dpdT_dv
(
const scalar rho,
const scalar T
) const
{
scalar Vm = this->W()/rho;
return this->RR()*log(Vm-b_)+dadT(T)*log(b_+Vm)/b_-dadT(T)*log(Vm)/b_;
return this->RR()*log(Vm - b_) + dadT(T)/b_*(log(b_ + Vm) - log(Vm));
//return this->RR()*log(Vm - b_) + dadT(T)*log(b_ + Vm)/b_ - dadT(T)*log(Vm)/b_;
}
@ -271,7 +276,6 @@ inline scalar soaveRedlichKwong::integral_dpdT_dv(const scalar rho,const scalar
inline scalar soaveRedlichKwong::d2pdT2(const scalar rho, const scalar T) const
{
scalar Vm = this->W()/rho;
return -d2adT2(T)/(Vm*(Vm + b_));
}
@ -281,48 +285,55 @@ inline scalar soaveRedlichKwong::d2pdv2(const scalar rho,const scalar T) const
{
scalar Vm = this->W()/rho;
scalar Vm2 = Vm*Vm;
scalar Vm3 = Vm*Vm*Vm;
scalar Vm3 = Vm2*Vm;
scalar Vm4 = Vm3*Vm;
scalar Vm5 = Vm4*Vm;
return
2*
return 2*
(
a(T)*
(
b5()-3*b3()*Vm2-b2()*Vm3+6*b_*pow(Vm,4)-3*pow(Vm,5)
)
+this->RR()*T*Vm3*
(
b3()+3*b2()*Vm+3*b_*Vm2+Vm3
)
a(T)*(b5_ - 3*b3_*Vm2 - b2_*Vm3 + 6*b_*Vm4 - 3*Vm5)
+ this->RR()*T*Vm3*(b3_ + 3*b2_*Vm + 3*b_*Vm2 + Vm3)
)
/(Vm3*pow(b_ + Vm, 3)*pow(Vm - b_, 3));
}
//(molar values)
// using second Order implicit differentiation
inline scalar soaveRedlichKwong::d2vdT2(const scalar rho, const scalar T) const
//using second order implicit differentiation
inline scalar soaveRedlichKwong::d2vdT2
(
const scalar rho,
const scalar T
) const
{
scalar dpdT2=this->dpdT(rho, T)*this->dpdT(rho, T);
scalar dpdv2=this->dpdv(rho, T)*this->dpdv(rho, T);
scalar dpdv3=dpdv2*this->dpdv(rho, T);
return
-(
pow(this->dpdT(rho,T),2)*this->d2pdv2(rho,T)
+ pow(this->dpdv(rho,T),2)*this->d2pdT2(rho,T)
dpdT2*this->d2pdv2(rho, T)
+ dpdv2*this->d2pdT2(rho, T)
- 2*this->dpdv(rho, T)*this->dpdT(rho, T)*this->d2pdvdT(rho, T)
)
/(pow(this->dpdv(rho,T),3));
/dpdv3;
}
//(molar values)
inline scalar soaveRedlichKwong::d2pdvdT(const scalar rho, const scalar T) const
inline scalar soaveRedlichKwong::d2pdvdT
(
const scalar rho,
const scalar T
) const
{
scalar Vm = this->W()/rho;
scalar Vm2 = Vm*Vm;
scalar Vm3 = Vm2*Vm;
return
(
dadT(T)*(b3()-3*b_*Vm2+2*pow(Vm,3))
-this->RR()*Vm2*(b2()+2*b_*Vm+Vm2)
dadT(T)*(b3_ - 3*b_*Vm2 + 2*Vm3) - this->RR()*Vm2*(b2_ + 2*b_*Vm + Vm2)
)
/(Vm2*pow(b_ + Vm, 2)*pow(b_ - Vm, 2));
}
@ -330,17 +341,24 @@ inline scalar soaveRedlichKwong::d2pdvdT(const scalar rho, const scalar T) const
// the result of this intergal is needed for the nasa based cp polynomial
//(molar values)
inline scalar soaveRedlichKwong::integral_d2pdT2_dv(const scalar rho,const scalar T) const
inline scalar soaveRedlichKwong::integral_d2pdT2_dv
(
const scalar rho,
const scalar T
) const
{
scalar Vm = this->W()/rho;
return d2adT2(T)*log(b_ + Vm)/b_ - d2adT2(T)*log(Vm)/b_;
}
//Isobar expansion Coefficent beta = 1/v (dv/dt) at constant p
//(molar values)
inline scalar soaveRedlichKwong::isobarExpCoef(const scalar rho,const scalar T) const
inline scalar soaveRedlichKwong::isobarExpCoef
(
const scalar rho,
const scalar T
) const
{
return this->dvdT(rho, T)*rho/this->W();
}
@ -348,20 +366,26 @@ inline scalar soaveRedlichKwong::isobarExpCoef(const scalar rho,const scalar T)
//isothemal compressiblity kappa (not Thermal conductivity)
//(molar values)
inline scalar soaveRedlichKwong::isothermalCompressiblity(const scalar rho,const scalar T) const
inline scalar soaveRedlichKwong::isothermalCompressiblity
(
const scalar rho,
const scalar T
) const
{
return this->isobarExpCoef(rho, T)/this->dpdT(rho, T);
//CL: also possible
//CL: return -this->dvdp(rho, T)*rho/this->W();
}
//- Return density [kg/m^3]
inline scalar soaveRedlichKwong::rho(
inline scalar soaveRedlichKwong::rho
(
const scalar p,
const scalar T,
const scalar rho0
) const
{
scalar molarVolumePrevIteration;
scalar molarVolume;
label iter = 0;
@ -379,6 +403,7 @@ inline scalar soaveRedlichKwong::rho(
label i=0;
do
{
//CL: modified Newton solver
molarVolume=molarVolumePrevIteration
-(
(this->p((this->W()/molarVolumePrevIteration),T) - p)
@ -389,8 +414,7 @@ inline scalar soaveRedlichKwong::rho(
if (i>8)
{
//CL: using bisection methode as backup,
//CL: solution must be between rho=0.001 to rho=1500;
//CL: if not, change rhoMax_ and rhoMin_
//CL: solution must be between rhoMin_ to rhoMax
for(i=0; i<200; i++)
{
scalar f1 = this->p(rho1, T) - p;
@ -434,7 +458,8 @@ inline scalar soaveRedlichKwong::rho(
{
FatalErrorIn
(
"inline scalar soaveRedlichKwong::rho(const scalar p, const scalar T, const scalar rho0) const "
"inline scalar soaveRedlichKwong::rho"
"(const scalar p, const scalar T, const scalar rho0) const "
) << "Maximum number of iterations exceeded"
<< abort(FatalError);
}
@ -444,7 +469,8 @@ inline scalar soaveRedlichKwong::rho(
return this->W()/molarVolume;
}
//- Return density [kg/m^3]on
//- Return density [kg/m^3]
inline scalar soaveRedlichKwong::rho(const scalar p, const scalar T) const
{
//CL: using perfect gas equation as starting point
@ -455,16 +481,22 @@ inline scalar soaveRedlichKwong::rho(const scalar p,const scalar T) const
//- Return compressibility drho/dp at T=constant [s^2/m^2]
inline scalar soaveRedlichKwong::psi(const scalar rho, const scalar T) const
{
return -this->dvdp(rho,T)*pow(rho,2)/this->W();
return -this->dvdp(rho, T)*rho*rho/this->W();
}
//- Return compression factor []
inline scalar soaveRedlichKwong::Z( const scalar p, const scalar T,const scalar rho0) const
inline scalar soaveRedlichKwong::Z
(
const scalar p,
const scalar T,
const scalar rho0
) const
{
return p/(this->R()*T*this->rho(p, T, rho0));
}
// * * * * * * * * * * * * * * * Member Operators * * * * * * * * * * * * * //
inline void soaveRedlichKwong::operator+=(const soaveRedlichKwong& srk)
@ -472,18 +504,13 @@ inline void soaveRedlichKwong::operator+=(const soaveRedlichKwong& srk)
specie::operator+=(srk);
}
/*
inline void soaveRedlichKwong::operator-=(const soaveRedlichKwong& srk)
{
specie::operator-=(srk);
}
inline void soaveRedlichKwong::operator*=(const scalar s)
{
specie::operator*=(s);
}
*/
// * * * * * * * * * * * * * * * Friend Operators * * * * * * * * * * * * * //
inline soaveRedlichKwong operator+
@ -499,6 +526,7 @@ inline soaveRedlichKwong operator+
);
}
inline soaveRedlichKwong operator*
(
const scalar s,
@ -508,6 +536,7 @@ inline soaveRedlichKwong operator*
return soaveRedlichKwong(s*static_cast<const specie&>(srk));
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam

71
src/thermophysicalModels/specie/thermo/realGasThermo/constantHeatCapacity/constantHeatCapacity.C Executable file
View file

@ -0,0 +1,71 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | foam-extend: Open Source CFD
\\ / O peration |
\\ / A nd | For copyright notice see file Copyright
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of foam-extend.
foam-extend 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.
foam-extend 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 foam-extend. If not, see <http://www.gnu.org/licenses/>.
Author
Christian Lucas
Institut für Thermodynamik
Technische Universität Braunschweig
Germany
\*---------------------------------------------------------------------------*/
#include "constantHeatCapacity.H"
#include "IOstreams.H"
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
template<class equationOfState>
Foam::constantHeatCapacity<equationOfState>::constantHeatCapacity(Istream& is)
:
equationOfState(is),
Cp0_(readScalar(is)),
cp0_(Cp0_*this->W()),
//values for some need terms at std
e0_std(e0(this->Tstd())),
s0_std(s0(this->Tstd())),
integral_p_dv_std(this->integral_p_dv(this->rhostd(),this->Tstd())),
integral_dpdT_dv_std(this->integral_dpdT_dv(this->rhostd(),this->Tstd())),
// cp @ STD (needed to limit cp for stability
cp_std(this->cp_nonLimited(this->rhostd(),this->Tstd()))
{
is.check("constantHeatCapacity::constantHeatCapacity(Istream& is)");
}
// * * * * * * * * * * * * * * * Ostream Operator * * * * * * * * * * * * * //
template<class equationOfState>
Foam::Ostream& Foam::operator<<
(
Ostream& os,
const constantHeatCapacity<equationOfState>& np
)
{
os << static_cast<const equationOfState&>(np) << tab
<< np.Cp0_;
os.check("Ostream& operator<<(Ostream& os, const constantHeatCapacity& np)");
return os;
}
// ************************************************************************* //

256
src/thermophysicalModels/specie/thermo/realGasThermo/constantHeatCapacity/constantHeatCapacity.H Executable file
View file

@ -0,0 +1,256 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | foam-extend: Open Source CFD
\\ / O peration |
\\ / A nd | For copyright notice see file Copyright
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of foam-extend.
foam-extend 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.
foam-extend 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 foam-extend. If not, see <http://www.gnu.org/licenses/>.
Class
Foam::constantHeatCapacity
Description
real gas thermodynamic class --> constant perfect gas heat capacity
Important: the perfect gas heat capacity is constant, the real heat capacity is not constant due to real gas corrections
templated into the equationOfState
-> uses the equation of state to calculate all real Gas properties like Enthalpy, Entropy ...
-> can not be used with the perfectGas equation of state
Equations for the real gas correction: Have a look at thermodnamics books e.g. Thermodynamics:
An Engineering Approch, 5 Edition, Chapter 12
SourceFiles
constantHeatCapacityI.H
constantHeatCapacity.C
Author
Christian Lucas
Institut für Thermodynamik
Technische Universität Braunschweig
Germany
\*---------------------------------------------------------------------------*/
#ifndef constantHeatCapacity_H
#define constantHeatCapacity_H
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
// Forward declaration of friend functions and operators
template<class equationOfState> class constantHeatCapacity;
template<class equationOfState>
inline constantHeatCapacity<equationOfState> operator+
(
const constantHeatCapacity<equationOfState>&,
const constantHeatCapacity<equationOfState>&
);
template<class equationOfState>
inline constantHeatCapacity<equationOfState> operator-
(
const constantHeatCapacity<equationOfState>&,
const constantHeatCapacity<equationOfState>&
);
template<class equationOfState>
inline constantHeatCapacity<equationOfState> operator*
(
const scalar,
const constantHeatCapacity<equationOfState>&
);
template<class equationOfState>
inline constantHeatCapacity<equationOfState> operator==
(
const constantHeatCapacity<equationOfState>&,
const constantHeatCapacity<equationOfState>&
);
template<class equationOfState>
Ostream& operator<<
(
Ostream&,
const constantHeatCapacity<equationOfState>&
);
/*---------------------------------------------------------------------------*\
Class constantHeatCapacity Thermo Declaration
\*---------------------------------------------------------------------------*/
template<class equationOfState>
class constantHeatCapacity
:
public equationOfState
{
// Private data
//CL: spec. cp
scalar Cp0_;
//CL: molar values
scalar cp0_;
scalar e0_std;
scalar s0_std;
scalar integral_p_dv_std;
scalar integral_dpdT_dv_std;
scalar cp_std;
// Private member functions
//- Construct from components
//CL: used for the operator+
inline constantHeatCapacity
(
const equationOfState& st,
const scalar cp0_
);
//- Construct from components
//CL: used for the operator*
inline constantHeatCapacity
(
const equationOfState& st,
const scalar cp0_,
const scalar e0_std_,
const scalar s0_std_,
const scalar integral_p_dv_std_,
const scalar integral_dpdT_dv_std_,
const scalar cp_std_
);
public:
//Variable
// Constructors
//- Construct from Istream
constantHeatCapacity(Istream&);
//- Construct as named copy
inline constantHeatCapacity(const word&, const constantHeatCapacity&);
//- Construct and return a clone
inline autoPtr<constantHeatCapacity> clone() const;
//- Selector from Istream
inline static autoPtr<constantHeatCapacity> New(Istream& is);
// Member Functions
//- perfect Gas Enthalpy [J/kmol]
inline scalar h0(const scalar T) const;
//- perfect Gas Entropy [J/(kmol K)]
inline scalar s0(const scalar T) const;
//- perfect Gas internal Energy [J/kmol]
inline scalar e0(const scalar T) const;
//- perfect gas Heat capacity at constant pressure [J/(kmol K)]
inline scalar cv0(const scalar T) const;
//- perfect gas Heat capacity at constant pressure [J/(kmol K)]
inline scalar cp0(const scalar T) const;
//- Limited Heat capacity at constant pressure [J/(kmol K)]
inline scalar cp(const scalar rho, const scalar T) const;
//- non Limited Heat capacity at constant pressure [J/(kmol K)]
inline scalar cp_nonLimited(const scalar rho, const scalar T) const;
//- Heat capacity at constant pressure [J/(kmol K)]
inline scalar cv(const scalar rho, const scalar T) const;
//- Enthalpy [J/kmol]
inline scalar h(const scalar rho, const scalar T) const;
//- Entropy [J/(kmol K)]
inline scalar s(const scalar rho,const scalar T) const;
//- Internal Energy [J/kmol]
inline scalar e(const scalar rho, const scalar T) const;
// Member operators
inline void operator+=(const constantHeatCapacity&);
inline void operator-=(const constantHeatCapacity&);
// Friend operators
friend constantHeatCapacity operator+ <equationOfState>
(
const constantHeatCapacity&,
const constantHeatCapacity&
);
friend constantHeatCapacity operator- <equationOfState>
(
const constantHeatCapacity&,
const constantHeatCapacity&
);
friend constantHeatCapacity operator* <equationOfState>
(
const scalar,
const constantHeatCapacity&
);
friend constantHeatCapacity operator== <equationOfState>
(
const constantHeatCapacity&,
const constantHeatCapacity&
);
// IOstream Operators
friend Ostream& operator<< <equationOfState>
(
Ostream&,
const constantHeatCapacity&
);
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#include "constantHeatCapacityI.H"
#ifdef NoRepository
# include "constantHeatCapacity.C"
#endif
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

407
src/thermophysicalModels/specie/thermo/realGasThermo/constantHeatCapacity/constantHeatCapacityI.H Executable file
View file

@ -0,0 +1,407 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | foam-extend: Open Source CFD
\\ / O peration |
\\ / A nd | For copyright notice see file Copyright
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of foam-extend.
foam-extend 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.
foam-extend 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 foam-extend. If not, see <http://www.gnu.org/licenses/>.
Author
Christian Lucas
Institut für Thermodynamik
Technische Universität Braunschweig
Germany
\*---------------------------------------------------------------------------*/
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
//- Construct from components
//CL: used for the operator+
template<class equationOfState>
inline Foam::constantHeatCapacity<equationOfState>::constantHeatCapacity
(
const equationOfState& st,
const scalar cp0_
)
:
equationOfState(st),
cp0_(cp0_),
e0_std(e0(this->Tstd)),
s0_std(s0(this->Tstd)),
integral_p_dv_std(this->integral_p_dv(this->rhostd(),this->Tstd)),
integral_dpdT_dv_std(this->integral_dpdT_dv(this->rhostd(),this->Tstd)),
cp_std(this->cp_nonLimited(this->rhostd(),this->Tstd))
{}
//- Construct from components
//CL: used for the operator*
template<class equationOfState>
inline Foam::constantHeatCapacity<equationOfState>::constantHeatCapacity
(
const equationOfState& st,
const scalar cp0_,
const scalar e0_std_,
const scalar s0_std_,
const scalar integral_p_dv_std_,
const scalar integral_dpdT_dv_std_,
const scalar cp_std_
)
:
equationOfState(st),
cp0_(cp0_),
e0_std(e0_std_),
s0_std(s0_std_),
integral_p_dv_std(integral_p_dv_std_),
integral_dpdT_dv_std(integral_dpdT_dv_std_),
cp_std(cp_std_)
{}
template<class equationOfState>
inline Foam::constantHeatCapacity<equationOfState>::constantHeatCapacity
(
const word& name,
const constantHeatCapacity& np
)
:
equationOfState(name, np),
cp0_(np.cp0_),
e0_std(np.e0_std),
s0_std(np.s0_std),
integral_p_dv_std(np.integral_p_dv_std),
integral_dpdT_dv_std(np.integral_dpdT_dv_std),
cp_std(np.cp_std)
{}
template<class equationOfState>
inline Foam::autoPtr<Foam::constantHeatCapacity<equationOfState> >
Foam::constantHeatCapacity<equationOfState>::clone() const
{
return autoPtr<constantHeatCapacity<equationOfState> >
(
new constantHeatCapacity<equationOfState>(*this)
);
}
template<class equationOfState>
inline Foam::autoPtr<Foam::constantHeatCapacity<equationOfState> >
Foam::constantHeatCapacity<equationOfState>::New(Istream& is)
{
return autoPtr<constantHeatCapacity<equationOfState> >
(
new constantHeatCapacity<equationOfState>(is)
);
}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
//used to calculate the internal energy
//perfect gas enthalpy
template<class equationOfState>
inline Foam::scalar Foam::constantHeatCapacity<equationOfState>::h0
(
const scalar T
) const
{
return cp0_*T;
}
//used to calculate the internal energy
//perfect gas internal energy
template<class equationOfState>
inline Foam::scalar Foam::constantHeatCapacity<equationOfState>::e0
(
const scalar T
) const
{
return this->h0(T) - this->RR()*T;
}
// used to calculate the entropy
// perfect gas entropy
template<class equationOfState>
inline Foam::scalar Foam::constantHeatCapacity<equationOfState>::s0
(
const scalar T
) const
{
return cp0_*log(T);
}
//perfect gas cp
template<class equationOfState>
inline Foam::scalar Foam::constantHeatCapacity<equationOfState>::cp0
(
const scalar T
) const
{
return cp0_;
}
//perfect gas cv
template<class equationOfState>
inline Foam::scalar Foam::constantHeatCapacity<equationOfState>::cv0
(
const scalar T
) const
{
return this->cp0(T)-this->RR();
}
//function to calculate real gas cp
//using cp=cv+(dp/dT)^2/(dp/dv)
template<class equationOfState>
inline Foam::scalar Foam::constantHeatCapacity<equationOfState>::cp
(
const scalar rho,
const scalar T
) const
{
// Problem --> dpdv(rho,T) is =0 at some points within the vapour dome. To increase stability, (dp/dv) has to be limited
// cp can be negative within the vapor dome. To avoid this nonphysical result, the absolute value is used.
// within the vapourdome and at the critical point, cp increases to very high values --> infinity,
// this would decrease the stability, so cp will be limited to 20 times the cp @ STD
return
min
(
cp_std*20,
fabs
(
this->cv(rho,T)
- T*pow((this->dpdT(rho, T)), 2)/min(this->dpdv(rho, T), -1)
)
);
}
// this function is needed to get cp @ STD (without the limit imposed in the function above),
// which in turn is needed to limit the cp in the function above
template<class equationOfState>
inline Foam::scalar Foam::constantHeatCapacity<equationOfState>::cp_nonLimited
(
const scalar rho,
const scalar T
) const
{
return
fabs(this->cv(rho, T)
- T*pow((this->dpdT(rho, T)), 2)/min(this->dpdv(rho, T), -1));
}
//function to calculate real gas c
//cv=cv0+T*integral d2p/dT2 dv
template<class equationOfState>
inline Foam::scalar Foam::constantHeatCapacity<equationOfState>::cv
(
const scalar rho,
const scalar T
) const
{
return this->cv0(T) + T*this->integral_d2pdT2_dv(rho, T);
}
//function to calculate real gas enthalpy
template<class equationOfState>
inline Foam::scalar Foam::constantHeatCapacity<equationOfState>::h
(
const scalar rho,
const scalar T
) const
{
return
this->e(rho, T)
+ this->p(rho, T)/rho*this->W()
- this->Pstd()/this->rhostd()*this->W();
}
// function to calculate real gas internal energy
// important assumption used: internal Energie is 0 at STD conditions.
// equation: du= cv0 dT +[T*dp/dT -p]dv
template<class equationOfState>
inline Foam::scalar Foam::constantHeatCapacity<equationOfState>::e
(
const scalar rho,
const scalar T
) const
{
return
- this->Tstd()*integral_dpdT_dv_std
+ integral_p_dv_std
+ this->e0(T)
- e0_std
+ T*this->integral_dpdT_dv(rho, T)
- this->integral_p_dv(rho, T);
}
//function to calculate real gas entropy
// important assumption used: the Entropy is 0 at STD conditions.
// equation: ds= cv0/T * dT + dp/dT *dv
// --> integral cv0/T dT = s0(T) -s0(Tstd) - R*ln(T/Tstd) --> due to s0(T)-s0(Tstd)=integral cp0/T dT
template<class equationOfState>
inline Foam::scalar Foam::constantHeatCapacity<equationOfState>::s
(
const scalar rho,
const scalar T
) const
{
return
- integral_dpdT_dv_std
+ this->s0(T)
- s0_std
- this->RR*log(T/this->Tstd())
+ this->integral_dpdT_dv(rho, T);
}
// * * * * * * * * * * * * * * * Member Operators * * * * * * * * * * * * * //
template<class equationOfState>
inline void Foam::constantHeatCapacity<equationOfState>::operator+=
(
const constantHeatCapacity<equationOfState>& np
)
{
scalar molr1 = this->nMoles();
equationOfState::operator+=(np);
molr1 /= this->nMoles();
scalar molr2 = np.nMoles()/this->nMoles();
cp0_ = molr1*cp0_ + molr2*np.cp0_;
}
template<class equationOfState>
inline void Foam::constantHeatCapacity<equationOfState>::operator-=
(
const constantHeatCapacity<equationOfState>& np
)
{
scalar molr1 = this->nMoles();
constantHeatCapacity::operator-=(np);
molr1 /= this->nMoles();
scalar molr2 = np.nMoles()/this->nMoles();
cp0_ = molr1*cp0_ - molr2*np.cp0_;
}
// * * * * * * * * * * * * * * * Friend Operators * * * * * * * * * * * * * //
template<class equationOfState>
inline Foam::constantHeatCapacity<equationOfState> Foam::operator+
(
const constantHeatCapacity<equationOfState>& np1,
const constantHeatCapacity<equationOfState>& np2
)
{
equationOfState eofs
(
static_cast<const equationOfState&>(np1)
+ static_cast<const equationOfState&>(np2)
);
//CL: Important, calls a different constructor as operator*
//CL: the coefficients as well as the EOS (coefficients) changed
//CL: therefore, the values at STD needs to be recalculated
return constantHeatCapacity<equationOfState>
(
eofs,
np1.nMoles()/eofs.nMoles()*np1.cp0_
+ np2.nMoles()/eofs.nMoles()*np2.cp0_
);
}
template<class equationOfState>
inline Foam::constantHeatCapacity<equationOfState> Foam::operator-
(
const constantHeatCapacity<equationOfState>& np1,
const constantHeatCapacity<equationOfState>& np2
)
{
equationOfState eofs
(
static_cast<const equationOfState&>(np1)
- static_cast<const equationOfState&>(np2)
);
return constantHeatCapacity<equationOfState>
(
eofs,
np1.nMoles()/eofs.nMoles()*np1.cp0_
- np2.nMoles()/eofs.nMoles()*np2.cp0_
);
}
template<class equationOfState>
inline Foam::constantHeatCapacity<equationOfState> Foam::operator*
(
const scalar s,
const constantHeatCapacity<equationOfState>& np
)
{
//CL: values at STD don't need to be recalculated,
//CL: therefore, providing the values in the constructor
return constantHeatCapacity<equationOfState>
(
s*static_cast<const equationOfState&>(np),
np.cp0_,
np.e0_std,
np.s0_std,
np.integral_p_dv_std,
np.integral_dpdT_dv_std,
np.cp_std
);
}
template<class equationOfState>
inline Foam::constantHeatCapacity<equationOfState> Foam::operator==
(
const constantHeatCapacity<equationOfState>& np1,
const constantHeatCapacity<equationOfState>& np2
)
{
return np2 - np1;
}
// ************************************************************************* //

3
src/thermophysicalModels/specie/thermo/realGasThermo/nasaHeatCapacityPolynomial/nasaHeatCapacityPolynomial.H
View file

@ -164,6 +164,9 @@ public:
//- Construct from Istream
nasaHeatCapacityPolynomial(Istream&);
//- Construct from dictionary
nasaHeatCapacityPolynomial(const dictionary& dict);
//- Construct as named copy
inline nasaHeatCapacityPolynomial(const word&, const nasaHeatCapacityPolynomial&);

10
src/thermophysicalModels/specie/thermo/realGasThermo/nasaHeatCapacityPolynomial/nasaHeatCapacityPolynomialI.H
View file

@ -146,6 +146,7 @@ Foam::nasaHeatCapacityPolynomial<equationOfState>::New(Istream& is)
);
}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
//used to calculate the internal energy
@ -200,8 +201,6 @@ inline Foam::scalar Foam::nasaHeatCapacityPolynomial<equationOfState>::s0
}
//perfect gas cp
template<class equationOfState>
inline Foam::scalar Foam::nasaHeatCapacityPolynomial<equationOfState>::cp0
@ -234,7 +233,6 @@ inline Foam::scalar Foam::nasaHeatCapacityPolynomial<equationOfState>::cv0
}
//function to calculate real gas cp
//using cp=cv+(dp/dT)^2/(dp/dv)
template<class equationOfState>
@ -281,15 +279,12 @@ inline Foam::scalar Foam::nasaHeatCapacityPolynomial<equationOfState>::cv
(
const scalar rho,
const scalar T
) const
{
return this->cv0(T)+T*this->integral_d2pdT2_dv(rho, T);
}
//function to calculate real gas enthalpy
template<class equationOfState>
inline Foam::scalar Foam::nasaHeatCapacityPolynomial<equationOfState>::h
@ -332,7 +327,6 @@ inline Foam::scalar Foam::nasaHeatCapacityPolynomial<equationOfState>::s
(
const scalar rho,
const scalar T
) const
{
return -integral_dpdT_dv_std
@ -341,8 +335,8 @@ inline Foam::scalar Foam::nasaHeatCapacityPolynomial<equationOfState>::s
+ this->integral_dpdT_dv(rho,T);
}
// * * * * * * * * * * * * * * * Member Operators * * * * * * * * * * * * * //
// * * * * * * * * * * * * * * * Member Operators * * * * * * * * * * * * * //
template<class equationOfState>
inline void Foam::nasaHeatCapacityPolynomial<equationOfState>::operator+=

4
src/thermophysicalModels/specie/thermo/realGasThermo/realGasSpecieThermo/realGasSpecieThermo.C
View file

@ -47,8 +47,8 @@ template<class thermo>
const Foam::debug::optimisationSwitch
Foam::realGasSpecieThermo<thermo>::maxIter_
(
"speciesThermoMaxIter",
100
"realGasSpecieThermoMaxIter",
500
);
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //

4
src/thermophysicalModels/specie/thermo/realGasThermo/realGasSpecieThermo/realGasSpecieThermo.H
View file

@ -34,8 +34,6 @@ Description
cp, h, s obtained from the template argument type thermo. All other
properties are derived from these primitive functions.
Some function copied from the "orginal" specieThermo!
SourceFiles
realGasSpecieThermoI.H
realGasSpecieThermo.C
@ -46,7 +44,6 @@ Institut für Thermodynamik
Technische Universität Braunschweig
Germany
\*---------------------------------------------------------------------------*/
#ifndef realGasSpecieThermo_H
@ -226,7 +223,6 @@ public:
inline void operator+=(const realGasSpecieThermo&);
inline void operator-=(const realGasSpecieThermo&);
inline void operator*=(const scalar);

7
src/thermophysicalModels/specie/thermo/realGasThermo/realGasSpecieThermo/realGasSpecieThermoI.H
View file

@ -28,7 +28,6 @@ Institut für Thermodynamik
Technische Universität Braunschweig
Germany
\*---------------------------------------------------------------------------*/
#include "realGasSpecieThermo.H"
@ -84,11 +83,11 @@ inline void Foam::realGasSpecieThermo<thermo>::T
(
(i<20)
&&
((
(
mag((this->*F)(rho, Tnew) - f)
>
mag((this->*F)(rho, Test) - f)
))
)
);
if (iter++ > maxIter_)
@ -114,6 +113,7 @@ inline void Foam::realGasSpecieThermo<thermo>::T
T0=Tnew;
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
template<class thermo>
@ -129,7 +129,6 @@ inline Foam::realGasSpecieThermo<thermo>::realGasSpecieThermo
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
template<class thermo>
inline Foam::scalar Foam::realGasSpecieThermo<thermo>::gamma(const scalar rho, const scalar T ) const
{

4
src/thermophysicalModels/specie/transport/const/constTransport.H
View file

@ -137,10 +137,10 @@ public:
//- Thermal diffusivity for enthalpy [kg/ms]
inline scalar alpha(const scalar T) const;
//- Thermal conductivity [W/mK]
//- Thermal conductivity [W/mK] for real gas
inline scalar kappa(const scalar rho, const scalar T) const;
//- Thermal diffusivity for enthalpy [kg/ms]
//- Thermal diffusivity for enthalpy [kg/ms] for real gas
inline scalar alpha(const scalar rho, const scalar T) const;
// Species diffusivity

11
src/thermophysicalModels/specie/transport/const/constTransportI.H
View file

@ -127,13 +127,20 @@ inline scalar constTransport<thermo>::kappa(const scalar rho,const scalar T) con
// CL: for real gas thermo
// Thermal diffusivity for enthalpy [kg/ms]
template<class thermo>
inline scalar constTransport<thermo>::alpha(const scalar rho,const scalar T) const
inline scalar constTransport<thermo>::alpha
(
const scalar rho,
const scalar T
) const
{
scalar Cp_ = this->Cp(rho, T);
scalar deltaT = T - specie::Tstd();
scalar CpBar =
(deltaT*(this->H(rho,T) - this->H(this->rhostd(),specie::Tstd())) + Cp_)/(sqr(deltaT) + 1);
(
deltaT*(this->H(rho, T) - this->H(this->rhostd(), specie::Tstd()))
+ Cp_
)/(sqr(deltaT) + 1);
return Cp_*mu(T)*rPr/CpBar;
}

56
src/thermophysicalModels/specie/transport/constRealGas/constRealGasTransport.C Normal file
View file

@ -0,0 +1,56 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | foam-extend: Open Source CFD
\\ / O peration |
\\ / A nd | For copyright notice see file Copyright
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of foam-extend.
foam-extend 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.
foam-extend 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 foam-extend. If not, see <http://www.gnu.org/licenses/>.
\*---------------------------------------------------------------------------*/
#include "constRealGasTransport.H"
#include "IOstreams.H"
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
template<class Thermo>
Foam::constRealGasTransport<Thermo>::constRealGasTransport(Istream& is)
:
Thermo(is),
mu_(readScalar(is)),
rPr_(1.0/readScalar(is))
{
is.check("constRealGasTransport::constRealGasTransport(Istream& is)");
}
// * * * * * * * * * * * * * * * IOstream Operators * * * * * * * * * * * * //
template<class Thermo>
Foam::Ostream& Foam::operator<<(Ostream& os, const constRealGasTransport<Thermo>& ct)
{
operator<<(os, static_cast<const Thermo&>(ct));
os << tab << ct.mu_ << tab << 1.0/ct.rPr_;
os.check("Ostream& operator<<(Ostream&, const constRealGasTransport&)");
return os;
}
// ************************************************************************* //

202
src/thermophysicalModels/specie/transport/constRealGas/constRealGasTransport.H Normal file
View file

@ -0,0 +1,202 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | foam-extend: Open Source CFD
\\ / O peration |
\\ / A nd | For copyright notice see file Copyright
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of foam-extend.
foam-extend 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.
foam-extend 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 foam-extend. If not, see <http://www.gnu.org/licenses/>.
Class
Foam::constRealGasTransport
Description
Constant properties Transport package.
Templated into a given thermodynamics package (needed for thermal
conductivity).
SourceFiles
constRealGasTransportI.H
constRealGasTransport.C
\*---------------------------------------------------------------------------*/
#ifndef constRealGasTransport_H
#define constRealGasTransport_H
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
// Forward declaration of friend functions and operators
template<class Thermo> class constRealGasTransport;
template<class Thermo>
inline constRealGasTransport<Thermo> operator+
(
const constRealGasTransport<Thermo>&,
const constRealGasTransport<Thermo>&
);
template<class Thermo>
inline constRealGasTransport<Thermo> operator-
(
const constRealGasTransport<Thermo>&,
const constRealGasTransport<Thermo>&
);
template<class Thermo>
inline constRealGasTransport<Thermo> operator*
(
const scalar,
const constRealGasTransport<Thermo>&
);
template<class Thermo>
inline constRealGasTransport<Thermo> operator==
(
const constRealGasTransport<Thermo>&,
const constRealGasTransport<Thermo>&
);
template<class Thermo>
Ostream& operator<<
(
Ostream&,
const constRealGasTransport<Thermo>&
);
/*---------------------------------------------------------------------------*\
Class constRealGasTransport Declaration
\*---------------------------------------------------------------------------*/
template<class Thermo>
class constRealGasTransport
:
public Thermo
{
// Private data
//- Constant dynamic viscosity [Pa.s]
scalar mu_;
//- Reciprocal Prandtl Number []
scalar rPr_;
// Private Member Functions
//- Construct from components
inline constRealGasTransport
(
const Thermo& t,
const scalar mu,
const scalar Pr
);
public:
// Constructors
//- Construct as named copy
inline constRealGasTransport(const word&, const constRealGasTransport&);
//- Construct from Istream
constRealGasTransport(Istream&);
// Member functions
//- Dynamic viscosity [kg/ms]
inline scalar mu(const scalar T) const;
//- Thermal conductivity [W/mK]
inline scalar kappa(const scalar rho, const scalar T) const;
//- Thermal diffusivity for enthalpy [kg/ms]
inline scalar alpha(const scalar rho, const scalar T) const;
// Species diffusivity
//inline scalar D(const scalar T) const;
// Member operators
inline constRealGasTransport& operator=
(
const constRealGasTransport&
);
// Friend operators
friend constRealGasTransport operator+ <Thermo>
(
const constRealGasTransport&,
const constRealGasTransport&
);
friend constRealGasTransport operator- <Thermo>
(
const constRealGasTransport&,
const constRealGasTransport&
);
friend constRealGasTransport operator* <Thermo>
(
const scalar,
const constRealGasTransport&
);
friend constRealGasTransport operator== <Thermo>
(
const constRealGasTransport&,
const constRealGasTransport&
);
// Ostream Operator
friend Ostream& operator<< <Thermo>
(
Ostream&,
const constRealGasTransport&
);
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#include "constRealGasTransportI.H"
#ifdef NoRepository
# include "constRealGasTransport.C"
#endif
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

217
src/thermophysicalModels/specie/transport/constRealGas/constRealGasTransportI.H Normal file
View file

@ -0,0 +1,217 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | foam-extend: Open Source CFD
\\ / O peration |
\\ / A nd | For copyright notice see file Copyright
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of foam-extend.
foam-extend 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.
foam-extend 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 foam-extend. If not, see <http://www.gnu.org/licenses/>.
\*---------------------------------------------------------------------------*/
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
template<class Thermo>
inline Foam::constRealGasTransport<Thermo>::constRealGasTransport
(
const Thermo& t,
const scalar mu,
const scalar Pr
)
:
Thermo(t),
mu_(mu),
rPr_(1.0/Pr)
{}
template<class Thermo>
inline Foam::constRealGasTransport<Thermo>::constRealGasTransport
(
const word& name,
const constRealGasTransport& ct
)
:
Thermo(name, ct),
mu_(ct.mu_),
rPr_(ct.rPr_)
{}
// Construct and return a clone
template<class thermo>
inline autoPtr<constRealGasTransport<Thermo> > constTransport<thermo>::clone
() const
{
return autoPtr<constRealGasTransport<Thermo> >
(
new constTransport<thermo>(*this)
);
}
// Selector from Istream
template<class thermo>
inline autoPtr<constRealGasTransport<Thermo> > constTransport<thermo>::New
(
Istream& is
)
{
return autoPtr<constRealGasTransport<Thermo> >
(
new constRealGasTransport<thermo>(is)
);
}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
template<class Thermo>
inline Foam::scalar Foam::constRealGasTransport<Thermo>::mu(const scalar) const
{
return mu_;
}
// CL: for real gas thermo
// Thermal conductivity [W/mK]
template<class thermo>
inline Foam::scalar Foam::constRealGasTransport<thermo>::kappa(const scalar rho, const scalar T) const
{
return this->Cp(rho, T)*mu(T)*rPr_;
}
// CL: for real gas thermo
// Thermal diffusivity for enthalpy [kg/ms]
template<class thermo>
inline Foam::scalar Foam::constRealGasTransport<thermo>::alpha(const scalar rho, const scalar T) const
{
scalar Cp_ = this->Cp(rho, T);
scalar deltaT = T - specie::Tstd();
scalar CpBar =
(deltaT*(this->H(rho, T) - this->H(this->rhostd(), specie::Tstd())) + Cp_)/(sqr(deltaT) + 1);
return Cp_*mu(T)*rPr_/CpBar;
}
// * * * * * * * * * * * * * * * Member Operators * * * * * * * * * * * * * //
template<class Thermo>
inline Foam::constRealGasTransport<Thermo>& Foam::constRealGasTransport<Thermo>::operator=
(
const constRealGasTransport<Thermo>& ct
)
{
Thermo::operator=(ct);
mu_ = ct.mu_;
rPr_ = ct.rPr_;
return *this;
}
// * * * * * * * * * * * * * * * Friend Operators * * * * * * * * * * * * * //
template<class Thermo>
inline Foam::constRealGasTransport<Thermo> Foam::operator+
(
const constRealGasTransport<Thermo>& ct1,
const constRealGasTransport<Thermo>& ct2
)
{
Thermo t
(
static_cast<const Thermo&>(ct1) + static_cast<const Thermo&>(ct2)
);
scalar molr1 = ct1.nMoles()/t.nMoles();
scalar molr2 = ct2.nMoles()/t.nMoles();
return constRealGasTransport<Thermo>
(
t,
molr1*ct1.mu_ + molr2*ct2.mu_,
molr1*ct1.rPr_ + molr2*ct2.rPr_
);
}
template<class Thermo>
inline Foam::constRealGasTransport<Thermo> Foam::operator-
(
const constRealGasTransport<Thermo>& ct1,
const constRealGasTransport<Thermo>& ct2
)
{
Thermo t
(
static_cast<const Thermo&>(ct1) - static_cast<const Thermo&>(ct2)
);
scalar molr1 = ct1.nMoles()/t.nMoles();
scalar molr2 = ct2.nMoles()/t.nMoles();
return constRealGasTransport<Thermo>
(
t,
molr1*ct1.mu_ - molr2*ct2.mu_,
molr1*ct1.rPr_ - molr2*ct2.rPr_
);
}
template<class Thermo>
inline Foam::constRealGasTransport<Thermo> Foam::operator*
(
const scalar s,
const constRealGasTransport<Thermo>& ct
)
{
return constRealGasTransport<Thermo>
(
s*static_cast<const Thermo&>(ct),
ct.mu_,
ct.rPr_
);
}
template<class Thermo>
inline Foam::constRealGasTransport<Thermo> Foam::operator==
(
const constRealGasTransport<Thermo>& ct1,
const constRealGasTransport<Thermo>& ct2
)
{
return ct2 - ct1;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// ************************************************************************* //

4
src/thermophysicalModels/specie/transport/sutherland/sutherlandTransport.H
View file

@ -159,10 +159,10 @@ public:
//- Thermal diffusivity for enthalpy [kg/ms]
inline scalar alpha(const scalar T) const;
//- Thermal conductivity [W/mK]
//- Thermal conductivity [W/mK] for real gas
inline scalar kappa(const scalar rho, const scalar T) const;
//- Thermal diffusivity for enthalpy [kg/ms]
//- Thermal diffusivity for enthalpy [kg/ms] for real gas
inline scalar alpha(const scalar rho, const scalar T) const;
// Species diffusivity

24
src/thermophysicalModels/specie/transport/sutherland/sutherlandTransportI.H
View file

@ -150,31 +150,45 @@ inline scalar sutherlandTransport<thermo>::alpha(const scalar T) const
scalar deltaT = T - specie::Tstd();
scalar CpBar =
(deltaT*(this->H(T) - this->H(specie::Tstd())) + Cp_)/(sqr(deltaT) + 1);
(deltaT*(this->H(T) - this->H(specie::Tstd())) + Cp_)
/(sqr(deltaT) + 1);
return mu(T)*Cv_*(1.32 + 1.77*this->R()/Cv_)/CpBar;
return mu(T)*Cv_*(1.32 + 1.77*R_/Cv_)/CpBar;
}
// CL: for real gas thermo
// Thermal conductivity [W/mK]
template<class thermo>
inline scalar sutherlandTransport<thermo>::kappa(const scalar rho, const scalar T) const
inline scalar sutherlandTransport<thermo>::kappa
(
const scalar rho,
const scalar T
) const
{
scalar Cv_ = this->Cv(rho,T);
return mu(T)*Cv_*(1.32 + 1.77*this->R()/Cv_);
}
// CL: for real gas thermo
// Thermal diffusivity for enthalpy [kg/ms]
template<class thermo>
inline scalar sutherlandTransport<thermo>::alpha(const scalar rho, const scalar T) const
inline scalar sutherlandTransport<thermo>::alpha
(
const scalar rho,
const scalar T
) const
{
scalar Cv_ = this->Cv(rho, T);
scalar Cp_ = this->Cp(rho, T);
scalar deltaT = T - specie::Tstd();
scalar CpBar =
(deltaT*(this->H(rho, T) - this->H(this->rhostd(),specie::Tstd())) + Cp_)/(sqr(deltaT) + 1);
(
deltaT*(this->H(rho, T) - this->H(this->rhostd(), specie::Tstd()))
+ Cp_
)/(sqr(deltaT) + 1);
return mu(T)*Cv_*(1.32 + 1.77*this->R()/Cv_)/CpBar;
}

2
tutorials/compressible/realFluidPisoFoam/ras/backStep/0/U
View file

@ -24,7 +24,7 @@ boundaryField
inlet
{
type pressureInletVelocity;
value uniform (3.5 0 0);
value uniform (1 0 0);
}
outlet
{

42
tutorials/compressible/realFluidPisoFoam/ras/backStep/0/epsilon
View file

@ -1,9 +1,9 @@
/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: 2.1.1 |
| \\ / A nd | Web: www.OpenFOAM.org |
| \\/ M anipulation | |
| \\ / F ield | foam-extend: Open Source CFD |
| \\ / O peration | Version: 4.0 |
| \\ / A nd | Web: http://www.foam-extend.org |
| \\/ M anipulation | For copyright notice see file Copyright |
\*---------------------------------------------------------------------------*/
FoamFile
{
@ -21,22 +21,6 @@ internalField uniform 1120;
boundaryField
{
upperWall
{
type compressible::epsilonWallFunction;
Cmu 0.09;
kappa 0.41;
E 9.8;
value uniform 1120;
}
lowerWall
{
type compressible::epsilonWallFunction;
Cmu 0.09;
kappa 0.41;
E 9.8;
value uniform 1120;
}
inlet
{
type fixedValue;
@ -48,6 +32,24 @@ boundaryField
inletValue uniform 1120;
value uniform 1120;
}
upperWall
{
type compressible::epsilonWallFunction;
refValue uniform 0;
value uniform 1120;
Cmu 0.09;
kappa 0.41;
E 9.8;
}
lowerWall
{
type compressible::epsilonWallFunction;
refValue uniform 0;
value uniform 1120;
Cmu 0.09;
kappa 0.41;
E 9.8;
}
frontAndBack
{
type empty;

28
tutorials/compressible/realFluidPisoFoam/ras/backStep/0/k
View file

@ -1,9 +1,9 @@
/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: 2.1.1 |
| \\ / A nd | Web: www.OpenFOAM.org |
| \\/ M anipulation | |
| \\ / F ield | foam-extend: Open Source CFD |
| \\ / O peration | Version: 4.0 |
| \\ / A nd | Web: http://www.foam-extend.org |
| \\/ M anipulation | For copyright notice see file Copyright |
\*---------------------------------------------------------------------------*/
FoamFile
{
@ -21,16 +21,6 @@ internalField uniform 5;
boundaryField
{
upperWall
{
type compressible::kqRWallFunction;
value uniform 5;
}
lowerWall
{
type compressible::kqRWallFunction;
value uniform 5;
}
inlet
{
type turbulentIntensityKineticEnergyInlet;
@ -43,6 +33,16 @@ boundaryField
{
type zeroGradient;
}
upperWall
{
type compressible::kqRWallFunction;
value uniform 5;
}
lowerWall
{
type compressible::kqRWallFunction;
value uniform 5;
}
frontAndBack
{
type empty;

12
tutorials/compressible/realFluidPisoFoam/ras/backStep/0/p
View file

@ -1,9 +1,9 @@
/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: 1.7.1 |
| \\ / A nd | Web: www.OpenFOAM.com |
| \\/ M anipulation | |
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | foam-extend: Open Source CFD
\\ / O peration |
\\ / A nd | For copyright notice see file Copyright
\\/ M anipulation |
\*---------------------------------------------------------------------------*/
FoamFile
{

12
tutorials/compressible/realFluidPisoFoam/ras/backStep/constant/RASProperties
View file

@ -1,9 +1,9 @@
/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: 1.6 |
| \\ / A nd | Web: www.OpenFOAM.org |
| \\/ M anipulation | |
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | foam-extend: Open Source CFD
\\ / O peration |
\\ / A nd | For copyright notice see file Copyright
\\/ M anipulation |
\*---------------------------------------------------------------------------*/
FoamFile
{

12
tutorials/compressible/realFluidPisoFoam/ras/backStep/constant/polyMesh/blockMeshDict
View file

@ -1,9 +1,9 @@
/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: 2.1.1 |
| \\ / A nd | Web: www.OpenFOAM.org |
| \\/ M anipulation | |
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | foam-extend: Open Source CFD
\\ / O peration |
\\ / A nd | For copyright notice see file Copyright
\\/ M anipulation |
\*---------------------------------------------------------------------------*/
FoamFile
{

8
tutorials/compressible/realFluidPisoFoam/ras/backStep/constant/polyMesh/boundary
View file

@ -1,9 +1,9 @@
/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: 2.1.1 |
| \\ / A nd | Web: www.OpenFOAM.org |
| \\/ M anipulation | |
| \\ / F ield | foam-extend: Open Source CFD |
| \\ / O peration | Version: 4.0 |
| \\ / A nd | Web: http://www.foam-extend.org |
| \\/ M anipulation | For copyright notice see file Copyright |
\*---------------------------------------------------------------------------*/
FoamFile
{

55
tutorials/compressible/realFluidPisoFoam/ras/backStep/constant/thermophysicalProperties
View file

@ -1,9 +1,9 @@
/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: 1.6 |
| \\ / A nd | Web: www.OpenFOAM.org |
| \\/ M anipulation | |
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | foam-extend: Open Source CFD
\\ / O peration |
\\ / A nd | For copyright notice see file Copyright
\\/ M anipulation |
\*---------------------------------------------------------------------------*/
FoamFile
{
@ -16,13 +16,13 @@ FoamFile
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
//CL: List of possible real gas models
thermoType realGasHThermo<pureMixture<sutherlandTransport<realGasSpecieThermo<nasaHeatCapacityPolynomial<redlichKwong>>>>>;
mixture CO2 1 44.01 73.773e5 304.13 49436.5054 -626.411601 5.30172524 0.002503813816 -0.0000002127308728 -0.000000000768998878 2.849677801e-13 1.4792e-06 116;
thermoType realGasHThermo<pureMixture<sutherlandTransport<realGasSpecieThermo<nasaHeatCapacityPolynomial<aungierRedlichKwong>>>>>;
mixture CO2 1 44.01 73.773e5 304.13 0.22394 467.6 49436.5054 -626.411601 5.30172524 0.002503813816 -0.0000002127308728 -0.000000000768998878 2.849677801E-13 1.4792e-06 116;
//thermoType realGasHThermo<pureMixture<sutherlandTransport<realGasSpecieThermo<nasaHeatCapacityPolynomial<redlichKwong>>>>>;
//mixture CO2 1 44.01 73.773e5 304.13 0.22394 49436.5054 -626.411601 5.30172524 0.002503813816 -0.0000002127308728 -0.000000000768998878 2.849677801e-13 1.4792e-06 116;
//thermoType realGasHThermo<pureMixture<sutherlandTransport<realGasSpecieThermo<nasaHeatCapacityPolynomial<aungierRedlichKwong>>>>>;
//mixture CO2 1 44.01 73.773e5 304.13 0.22394 467.6 49436.5054 -626.411601 5.30172524 0.002503813816 -0.0000002127308728 -0.000000000768998878 2.849677801E-13 1.4792e-06 116;
//thermoType realGasHThermo<pureMixture<sutherlandTransport<realGasSpecieThermo<nasaHeatCapacityPolynomial<pengRobinson>>>>>;
//mixture CO2 1 44.01 73.773e5 304.13 0.22394 49436.5054 -626.411601 5.30172524 0.002503813816 -0.0000002127308728 -0.000000000768998878 2.849677801e-13 1.4792e-06 116;
@ -30,30 +30,39 @@ mixture CO2 1 44.01 73.773e5 304.13 0.22394 467.6 49436.5054 -626.411601
//thermoType realGasHThermo<pureMixture<sutherlandTransport<realGasSpecieThermo<nasaHeatCapacityPolynomial<soaveRedlichKwong>>>>>;
//mixture CO2 1 44.01 73.773e5 304.13 0.22394 49436.5054 -626.411601 5.30172524 0.002503813816 -0.0000002127308728 -0.000000000768998878 2.849677801e-13 1.4792e-06 116;
//thermoType realGasHThermo<pureMixture<constTransport<realGasSpecieThermo<nasaHeatCapacityPolynomial<aungierRedlichKwong>>>>>;
//mixture CO2 1 44.01 73.773e5 304.13 0.22394 467.6 49436.5054 -626.411601 5.30172524 0.002503813816 -0.0000002127308728 -0.000000000768998878 2.849677801E-13 1e-06 0.7;
//thermoType realGasHThermo<pureMixture<constTransport<realGasSpecieThermo<nasaHeatCapacityPolynomial<redlichKwong>>>>>;
//mixture CO2 1 44.01 73.773e5 304.13 0.22394 49436.5054 -626.411601 5.30172524 0.002503813816 -0.0000002127308728 -0.000000000768998878 2.849677801e-13 1e-06 0.7;
//thermoType realGasHThermo<pureMixture<sutherlandTransport<realGasSpecieThermo<constantHeatCapacity<redlichKwong>>>>>;
//mixture CO2 1 44.01 73.773e5 304.13 839 1.4792e-06 116;
//thermoType realGasHThermo<pureMixture<constTransport<realGasSpecieThermo<nasaHeatCapacityPolynomial<pengRobinson>>>>>;
//mixture CO2 1 44.01 73.773e5 304.13 0.22394 49436.5054 -626.411601 5.30172524 0.002503813816 -0.0000002127308728 -0.000000000768998878 2.849677801e-13 1e-06 0.7;
//thermoType realGasHThermo<pureMixture<sutherlandTransport<realGasSpecieThermo<constantHeatCapacity<aungierRedlichKwong>>>>>;
//thermoType realGasHThermo<pureMixture<constTransport<realGasSpecieThermo<nasaHeatCapacityPolynomial<soaveRedlichKwong>>>>>;
//mixture CO2 1 44.01 73.773e5 304.13 0.22394 49436.5054 -626.411601 5.30172524 0.002503813816 -0.0000002127308728 -0.000000000768998878 2.849677801e-13 1e-06 0.7;
//thermoType realGasHThermo<pureMixture<sutherlandTransport<realGasSpecieThermo<constantHeatCapacity<pengRobinson>>>>>;
//thermoType realGasHThermo<pureMixture<sutherlandTransport<realGasSpecieThermo<constantHeatCapacity<soaveRedlichKwong>>>>>;
//thermoType realGasHThermo<pureMixture<constRealGasTransport<realGasSpecieThermo<nasaHeatCapacityPolynomial<redlichKwong>>>>>;
//thermoType realGasHThermo<pureMixture<constRealGasTransport<realGasSpecieThermo<nasaHeatCapacityPolynomial<pengRobinson>>>>>;
//thermoType realGasHThermo<pureMixture<constRealGasTransport<realGasSpecieThermo<nasaHeatCapacityPolynomial<aungierRedlichKwong>>>>>;
//thermoType realGasHThermo<pureMixture<constRealGasTransport<realGasSpecieThermo<nasaHeatCapacityPolynomial<soaveRedlichKwong>>>>>;
//thermoType realGasHThermo<pureMixture<constRealGasTransport<realGasSpecieThermo<constantHeatCapacity<redlichKwong>>>>>;
//thermoType realGasHThermo<pureMixture<constRealGasTransport<realGasSpecieThermo<constantHeatCapacity<pengRobinson>>>>>;
//thermoType realGasHThermo<pureMixture<constRealGasTransport<realGasSpecieThermo<constantHeatCapacity<aungierRedlichKwong>>>>>;
//thermoType realGasHThermo<pureMixture<constRealGasTransport<realGasSpecieThermo<constantHeatCapacity<soaveRedlichKwong>>>>>;
//CL: description of coefficients
// *********************************************************************************************************************** //
// Coefficient
// Coefficient:
// CO2 --> Name
// 1
// 44.01 --> Molar Volume
// 77.773e5 --> critical pressure
// 304.13 --> critical temperatur
// 0.22394 --> acentric factor
// 0.22394 --> acentric factor (not needed for redlich kwong)
// 467.6 --> critical density (only for aungier redlich kwong)
// 49436.5054 --> 2.849677801e-13 --> 7 heat capacity polynomial coefficent's
// .... --> two coefficent's for sutherland model or for the constTransport model
// .... --> two coefficent's for sutherlandTransport or for the constRealGasTransport model
// 839 --> perfect gas heat capacity Cp0 (J/kgK), needed for constantHeatCapacity
// *********************************************************************************************************************** //
// ************************************************************************* //

12
tutorials/compressible/realFluidPisoFoam/ras/backStep/constant/turbulenceProperties
View file

@ -1,9 +1,9 @@
/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: 1.6 |
| \\ / A nd | Web: www.OpenFOAM.org |
| \\/ M anipulation | |
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | foam-extend: Open Source CFD
\\ / O peration |
\\ / A nd | For copyright notice see file Copyright
\\/ M anipulation |
\*---------------------------------------------------------------------------*/
FoamFile
{

16
tutorials/compressible/realFluidPisoFoam/ras/backStep/system/controlDict
View file

@ -1,9 +1,9 @@
/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: 1.6 |
| \\ / A nd | Web: www.OpenFOAM.org |
| \\/ M anipulation | |
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | foam-extend: Open Source CFD
\\ / O peration |
\\ / A nd | For copyright notice see file Copyright
\\/ M anipulation |
\*---------------------------------------------------------------------------*/
FoamFile
{
@ -19,6 +19,8 @@ application realFluidPisoFoam;
startFrom latestTime;
//startTime 0;
stopAt endTime;
endTime 0.5;
@ -27,7 +29,7 @@ deltaT 1e-5;
writeControl runTime;
writeInterval 1e-2;
writeInterval 0.1;
purgeWrite 0;

12
tutorials/compressible/realFluidPisoFoam/ras/backStep/system/fvSchemes
View file

@ -1,9 +1,9 @@
/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM Extend Project: Open Source CFD |
| \\ / O peration | Version: 1.6-ext |
| \\ / A nd | Web: www.extend-project.de |
| \\/ M anipulation | |
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | foam-extend: Open Source CFD
\\ / O peration |
\\ / A nd | For copyright notice see file Copyright
\\/ M anipulation |
\*---------------------------------------------------------------------------*/
FoamFile
{

14
tutorials/compressible/realFluidPisoFoam/ras/backStep/system/fvSolution
View file

@ -1,9 +1,9 @@
/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | |
| \\ / A nd |
| \\/ M anipulation | |
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | foam-extend: Open Source CFD
\\ / O peration |
\\ / A nd | For copyright notice see file Copyright
\\/ M anipulation |
\*---------------------------------------------------------------------------*/
FoamFile
{
@ -17,7 +17,6 @@ FoamFile
solvers
{
p
{
solver GAMG;
@ -88,7 +87,6 @@ p
PISO
{
realFluid true;
nNonOrthogonalCorrectors 0;
nCorrectors 2;
momentumPredictor yes;

49
tutorials/compressible/realFluidPisoFoam/ras/backStep_IAPWS97/0.org/T Executable file
View file

@ -0,0 +1,49 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | foam-extend: Open Source CFD
\\ / O peration |
\\ / A nd | For copyright notice see file Copyright
\\/ M anipulation |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class volScalarField;
location "0";
object T;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
dimensions [0 0 0 1 0 0 0];
internalField uniform 500;
boundaryField
{
inlet
{
type fixedValue;
value uniform 500;
}
outlet
{
type zeroGradient;
}
upperWall
{
type zeroGradient;
}
lowerWall
{
type zeroGradient;
}
frontAndBack
{
type empty;
}
}
// ************************************************************************* //

51
tutorials/compressible/realFluidPisoFoam/ras/backStep_IAPWS97/0.org/U Executable file
View file

@ -0,0 +1,51 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | foam-extend: Open Source CFD
\\ / O peration |
\\ / A nd | For copyright notice see file Copyright
\\/ M anipulation |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class volVectorField;
location "0";
object U;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
dimensions [0 1 -1 0 0 0 0];
internalField uniform (0 0 0);
boundaryField
{
inlet
{
type fixedValue;
value uniform (1 0 0);
}
outlet
{
type zeroGradient;
}
upperWall
{
type fixedValue;
value uniform (0 0 0);
}
lowerWall
{
type fixedValue;
value uniform (0 0 0);
}
frontAndBack
{
type empty;
}
}
// ************************************************************************* //

58
tutorials/compressible/realFluidPisoFoam/ras/backStep_IAPWS97/0.org/epsilon Normal file
View file

@ -0,0 +1,58 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | foam-extend: Open Source CFD
\\ / O peration |
\\ / A nd | For copyright notice see file Copyright
\\/ M anipulation |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class volScalarField;
location "0";
object epsilon;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
dimensions [0 2 -3 0 0 0 0];
internalField uniform 1120;
boundaryField
{
upperWall
{
type compressible::epsilonWallFunction;
Cmu 0.09;
kappa 0.41;
E 9.8;
value uniform 1120;
}
lowerWall
{
type compressible::epsilonWallFunction;
Cmu 0.09;
kappa 0.41;
E 9.8;
value uniform 1120;
}
inlet
{
type fixedValue;
value uniform 1120;
}
outlet
{
type inletOutlet;
inletValue uniform 1120;
value uniform 1120;
}
frontAndBack
{
type empty;
}
}
// ************************************************************************* //

53
tutorials/compressible/realFluidPisoFoam/ras/backStep_IAPWS97/0.org/k Normal file
View file

@ -0,0 +1,53 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | foam-extend: Open Source CFD
\\ / O peration |
\\ / A nd | For copyright notice see file Copyright
\\/ M anipulation |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class volScalarField;
location "0";
object k;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
dimensions [0 2 -2 0 0 0 0];
internalField uniform 5;
boundaryField
{
upperWall
{
type compressible::kqRWallFunction;
value uniform 5;
}
lowerWall
{
type compressible::kqRWallFunction;
value uniform 5;
}
inlet
{
type turbulentIntensityKineticEnergyInlet;
intensity 0.01;
U U;
phi phi;
value uniform 5;
}
outlet
{
type zeroGradient;
}
frontAndBack
{
type empty;
}
}
// ************************************************************************* //

49
tutorials/compressible/realFluidPisoFoam/ras/backStep_IAPWS97/0.org/p Executable file
View file

@ -0,0 +1,49 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | foam-extend: Open Source CFD
\\ / O peration |
\\ / A nd | For copyright notice see file Copyright
\\/ M anipulation |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class volScalarField;
location "0";
object p;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
dimensions [1 -1 -2 0 0 0 0];
internalField uniform 10e5;
boundaryField
{
inlet
{
type zeroGradient;
}
outlet
{
type fixedValue;
value uniform 1e+06;
}
upperWall
{
type zeroGradient;
}
lowerWall
{
type zeroGradient;
}
frontAndBack
{
type empty;
}
}
// ************************************************************************* //

49
tutorials/compressible/realFluidPisoFoam/ras/backStep_IAPWS97/0/T Executable file
View file

@ -0,0 +1,49 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | foam-extend: Open Source CFD
\\ / O peration |
\\ / A nd | For copyright notice see file Copyright
\\/ M anipulation |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class volScalarField;
location "0";
object T;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
dimensions [0 0 0 1 0 0 0];
internalField uniform 500;
boundaryField
{
inlet
{
type fixedValue;
value uniform 500;
}
outlet
{
type zeroGradient;
}
upperWall
{
type zeroGradient;
}
lowerWall
{
type zeroGradient;
}
frontAndBack
{
type empty;
}
}
// ************************************************************************* //

51
tutorials/compressible/realFluidPisoFoam/ras/backStep_IAPWS97/0/U Executable file
View file

@ -0,0 +1,51 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | foam-extend: Open Source CFD
\\ / O peration |
\\ / A nd | For copyright notice see file Copyright
\\/ M anipulation |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class volVectorField;
location "0";
object U;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
dimensions [0 1 -1 0 0 0 0];
internalField uniform (0 0 0);
boundaryField
{
inlet
{
type fixedValue;
value uniform (1 0 0);
}
outlet
{
type zeroGradient;
}
upperWall
{
type fixedValue;
value uniform (0 0 0);
}
lowerWall
{
type fixedValue;
value uniform (0 0 0);
}
frontAndBack
{
type empty;
}
}
// ************************************************************************* //

60
tutorials/compressible/realFluidPisoFoam/ras/backStep_IAPWS97/0/epsilon Normal file
View file

@ -0,0 +1,60 @@
/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | foam-extend: Open Source CFD |
| \\ / O peration | Version: 4.0 |
| \\ / A nd | Web: http://www.foam-extend.org |
| \\/ M anipulation | For copyright notice see file Copyright |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class volScalarField;
location "0";
object epsilon;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
dimensions [0 2 -3 0 0 0 0];
internalField uniform 1120;
boundaryField
{
inlet
{
type fixedValue;
value uniform 1120;
}
outlet
{
type inletOutlet;
inletValue uniform 1120;
value uniform 1120;
}
upperWall
{
type compressible::epsilonWallFunction;
refValue uniform 0;
value uniform 1120;
Cmu 0.09;
kappa 0.41;
E 9.8;
}
lowerWall
{
type compressible::epsilonWallFunction;
refValue uniform 0;
value uniform 1120;
Cmu 0.09;
kappa 0.41;
E 9.8;
}
frontAndBack
{
type empty;
}
}
// ************************************************************************* //

53
tutorials/compressible/realFluidPisoFoam/ras/backStep_IAPWS97/0/k Normal file
View file

@ -0,0 +1,53 @@
/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | foam-extend: Open Source CFD |
| \\ / O peration | Version: 4.0 |
| \\ / A nd | Web: http://www.foam-extend.org |
| \\/ M anipulation | For copyright notice see file Copyright |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class volScalarField;
location "0";
object k;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
dimensions [0 2 -2 0 0 0 0];
internalField uniform 5;
boundaryField
{
inlet
{
type turbulentIntensityKineticEnergyInlet;
intensity 0.01;
U U;
phi phi;
value uniform 5;
}
outlet
{
type zeroGradient;
}
upperWall
{
type compressible::kqRWallFunction;
value uniform 5;
}
lowerWall
{
type compressible::kqRWallFunction;
value uniform 5;
}
frontAndBack
{
type empty;
}
}
// ************************************************************************* //

49
tutorials/compressible/realFluidPisoFoam/ras/backStep_IAPWS97/0/p Executable file
View file

@ -0,0 +1,49 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | foam-extend: Open Source CFD
\\ / O peration |
\\ / A nd | For copyright notice see file Copyright
\\/ M anipulation |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class volScalarField;
location "0";
object p;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
dimensions [1 -1 -2 0 0 0 0];
internalField uniform 10e5;
boundaryField
{
inlet
{
type zeroGradient;
}
outlet
{
type fixedValue;
value uniform 1e+06;
}
upperWall
{
type zeroGradient;
}
lowerWall
{
type zeroGradient;
}
frontAndBack
{
type empty;
}
}
// ************************************************************************* //

9
tutorials/compressible/realFluidPisoFoam/ras/backStep_IAPWS97/Allrun Executable file
View file

@ -0,0 +1,9 @@
#!/bin/sh
cd ${0%/*} || exit 1 # run from this directory
# Source tutorial run functions
. $WM_PROJECT_DIR/bin/tools/RunFunctions
runApplication blockMesh
runApplication `getApplication`

25
tutorials/compressible/realFluidPisoFoam/ras/backStep_IAPWS97/constant/RASProperties Executable file
View file

@ -0,0 +1,25 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | foam-extend: Open Source CFD
\\ / O peration |
\\ / A nd | For copyright notice see file Copyright
\\/ M anipulation |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
location "constant";
object RASProperties;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
RASModel kEpsilon;
turbulence on;
printCoeffs on;
// ************************************************************************* //

173
tutorials/compressible/realFluidPisoFoam/ras/backStep_IAPWS97/constant/polyMesh/blockMeshDict Normal file
View file

@ -0,0 +1,173 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | foam-extend: Open Source CFD
\\ / O peration |
\\ / A nd | For copyright notice see file Copyright
\\/ M anipulation |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
object blockMeshDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
convertToMeters 0.001;
vertices
(
(-20.6 0 -0.5)
(-20.6 3 -0.5)
(-20.6 12.7 -0.5)
(-20.6 25.4 -0.5)
(0 -25.4 -0.5)
(0 -5 -0.5)
(0 0 -0.5)
(0 3 -0.5)
(0 12.7 -0.5)
(0 25.4 -0.5)
(206 -25.4 -0.5)
(206 -8.5 -0.5)
(206 0 -0.5)
(206 6.5 -0.5)
(206 17 -0.5)
(206 25.4 -0.5)
(290 -16.6 -0.5)
(290 -6.3 -0.5)
(290 0 -0.5)
(290 4.5 -0.5)
(290 11 -0.5)
(290 16.6 -0.5)
(-20.6 0 0.5)
(-20.6 3 0.5)
(-20.6 12.7 0.5)
(-20.6 25.4 0.5)
(0 -25.4 0.5)
(0 -5 0.5)
(0 0 0.5)
(0 3 0.5)
(0 12.7 0.5)
(0 25.4 0.5)
(206 -25.4 0.5)
(206 -8.5 0.5)
(206 0 0.5)
(206 6.5 0.5)
(206 17 0.5)
(206 25.4 0.5)
(290 -16.6 0.5)
(290 -6.3 0.5)
(290 0 0.5)
(290 4.5 0.5)
(290 11 0.5)
(290 16.6 0.5)
);
blocks
(
hex (0 6 7 1 22 28 29 23) (18 7 1) simpleGrading (0.5 1.8 1)
hex (1 7 8 2 23 29 30 24) (18 10 1) simpleGrading (0.5 4 1)
hex (2 8 9 3 24 30 31 25) (18 13 1) simpleGrading (0.5 0.25 1)
hex (4 10 11 5 26 32 33 27) (180 18 1) simpleGrading (4 1 1)
hex (5 11 12 6 27 33 34 28) (180 9 1) edgeGrading (4 4 4 4 0.5 1 1 0.5 1 1 1 1)
hex (6 12 13 7 28 34 35 29) (180 7 1) edgeGrading (4 4 4 4 1.8 1 1 1.8 1 1 1 1)
hex (7 13 14 8 29 35 36 30) (180 10 1) edgeGrading (4 4 4 4 4 1 1 4 1 1 1 1)
hex (8 14 15 9 30 36 37 31) (180 13 1) simpleGrading (4 0.25 1)
hex (10 16 17 11 32 38 39 33) (25 18 1) simpleGrading (2.5 1 1)
hex (11 17 18 12 33 39 40 34) (25 9 1) simpleGrading (2.5 1 1)
hex (12 18 19 13 34 40 41 35) (25 7 1) simpleGrading (2.5 1 1)
hex (13 19 20 14 35 41 42 36) (25 10 1) simpleGrading (2.5 1 1)
hex (14 20 21 15 36 42 43 37) (25 13 1) simpleGrading (2.5 0.25 1)
);
edges
(
);
boundary
(
inlet
{
type patch;
faces
(
(0 22 23 1)
(1 23 24 2)
(2 24 25 3)
);
}
outlet
{
type patch;
faces
(
(16 17 39 38)
(17 18 40 39)
(18 19 41 40)
(19 20 42 41)
(20 21 43 42)
);
}
upperWall
{
type wall;
faces
(
(3 25 31 9)
(9 31 37 15)
(15 37 43 21)
);
}
lowerWall
{
type wall;
faces
(
(0 6 28 22)
(6 5 27 28)
(5 4 26 27)
(4 10 32 26)
(10 16 38 32)
);
}
frontAndBack
{
type empty;
faces
(
(22 28 29 23)
(23 29 30 24)
(24 30 31 25)
(26 32 33 27)
(27 33 34 28)
(28 34 35 29)
(29 35 36 30)
(30 36 37 31)
(32 38 39 33)
(33 39 40 34)
(34 40 41 35)
(35 41 42 36)
(36 42 43 37)
(0 1 7 6)
(1 2 8 7)
(2 3 9 8)
(4 5 11 10)
(5 6 12 11)
(6 7 13 12)
(7 8 14 13)
(8 9 15 14)
(10 11 17 16)
(11 12 18 17)
(12 13 19 18)
(13 14 20 19)
(14 15 21 20)
);
}
);
mergePatchPairs
(
);
// ************************************************************************* //

52
tutorials/compressible/realFluidPisoFoam/ras/backStep_IAPWS97/constant/polyMesh/boundary Normal file
View file

@ -0,0 +1,52 @@
/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | foam-extend: Open Source CFD |
| \\ / O peration | Version: 4.0 |
| \\ / A nd | Web: http://www.foam-extend.org |
| \\/ M anipulation | For copyright notice see file Copyright |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class polyBoundaryMesh;
location "constant/polyMesh";
object boundary;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
5
(
inlet
{
type patch;
nFaces 30;
startFace 24170;
}
outlet
{
type patch;
nFaces 57;
startFace 24200;
}
upperWall
{
type wall;
nFaces 223;
startFace 24257;
}
lowerWall
{
type wall;
nFaces 250;
startFace 24480;
}
frontAndBack
{
type empty;
nFaces 24450;
startFace 24730;
}
)
// ************************************************************************* //

23
tutorials/compressible/realFluidPisoFoam/ras/backStep_IAPWS97/constant/thermophysicalProperties Executable file
View file

@ -0,0 +1,23 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | foam-extend: Open Source CFD
\\ / O peration |
\\ / A nd | For copyright notice see file Copyright
\\/ M anipulation |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
location "constant";
object thermophysicalProperties;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
//CL: this is all
thermoType IAPWSThermo;
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

21
tutorials/compressible/realFluidPisoFoam/ras/backStep_IAPWS97/constant/turbulenceProperties Executable file
View file

@ -0,0 +1,21 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | foam-extend: Open Source CFD
\\ / O peration |
\\ / A nd | For copyright notice see file Copyright
\\/ M anipulation |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
location "constant";
object turbulenceProperties;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
simulationType RASModel;
// ************************************************************************* //

61
tutorials/compressible/realFluidPisoFoam/ras/backStep_IAPWS97/system/controlDict Executable file
View file

@ -0,0 +1,61 @@
/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: 1.6 |
| \\ / A nd | Web: www.OpenFOAM.org |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
location "system";
object controlDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
application realFluidPisoFoam;
startFrom latestTime;
//startTime 0;
stopAt endTime;
endTime 0.5;
deltaT 1e-5;
writeControl runTime;
writeInterval 1e-1;
purgeWrite 0;
writeFormat ascii;
writePrecision 6;
writeCompression uncompressed;
timeFormat general;
timePrecision 10;
adjustTimeStep yes;
maxCo 0.75;
maxDeltaT 0.01;
runTimeModifiable yes;
libs
(
"libIAPWSThermo.so"
"libfreesteam.so"
);
// ************************************************************************* //

70
tutorials/compressible/realFluidPisoFoam/ras/backStep_IAPWS97/system/fvSchemes Normal file
View file

@ -0,0 +1,70 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | foam-extend: Open Source CFD
\\ / O peration |
\\ / A nd | For copyright notice see file Copyright
\\/ M anipulation |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
location "system";
object fvSchemes;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
ddtSchemes
{
default Euler;
}
gradSchemes
{
default Gauss linear;
grad(p) Gauss linear;
}
divSchemes
{
default none;
div(phi,U) Gauss limitedLinearV 1;
div(phid,p) Gauss limitedLinear 1;
div(phiU,p) Gauss linear;
div(phi,h) Gauss limitedLinear 1;
div(phi,k) Gauss limitedLinear 1;
div(phi,epsilon) Gauss limitedLinear 1;
div(phi,R) Gauss limitedLinear 1;
div(phi,omega) Gauss limitedLinear 1;
div((rho*R)) Gauss linear;
div(R) Gauss linear;
div(U) Gauss linear;
div((muEff*dev2(T(grad(U))))) Gauss linear;
}
laplacianSchemes
{
default none;
laplacian(muEff,U) Gauss linear corrected;
laplacian(mut,U) Gauss linear corrected;
laplacian(DkEff,k) Gauss linear corrected;
laplacian(DepsilonEff,epsilon) Gauss linear corrected;
laplacian(DREff,R) Gauss linear corrected;
laplacian(DomegaEff,omega) Gauss linear corrected;
laplacian((rho*(1|A(U))),p) Gauss linear corrected;
laplacian(alphaEff,h) Gauss linear corrected;
}
interpolationSchemes
{
default linear;
}
snGradSchemes
{
default corrected;
}
// ************************************************************************* //

96
tutorials/compressible/realFluidPisoFoam/ras/backStep_IAPWS97/system/fvSolution Executable file
View file

@ -0,0 +1,96 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | foam-extend: Open Source CFD
\\ / O peration |
\\ / A nd | For copyright notice see file Copyright
\\/ M anipulation |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
location "system";
object fvSolution;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
solvers
{
p
{
solver GAMG;
tolerance 1e-14;
relTol 0.001;
smoother GaussSeidel;
minIter 4;
//maxIter 100;
cacheAgglomeration true;
nPreSweeps 1;
nPostSweeps 3;
nFinestSweeps 3;
nCellsInCoarsestLevel 20;
agglomerator faceAreaPair;
mergeLevels 1;
}
U
{
solver smoothSolver;
smoother GaussSeidel;
nSweeps 2;
tolerance 1e-14;
relTol 0.001;
}
rho
{
solver PCG;
preconditioner DIC;
tolerance 1e-10;
relTol 0;
}
htot
{
solver smoothSolver;
smoother GaussSeidel;
nSweeps 2;
tolerance 1e-14;
relTol 0.001;
}
h
{
solver smoothSolver;
smoother GaussSeidel;
nSweeps 2;
tolerance 1e-14;
relTol 0.001;
}
k
{
solver PBiCG;
preconditioner DILU;
tolerance 1e-10;
relTol 0;
}
epsilon
{
solver PBiCG;
preconditioner DILU;
tolerance 1e-10;
relTol 0;
}
}
PISO
{
nNonOrthogonalCorrectors 0;
nCorrectors 2;
momentumPredictor yes;
}
// ************************************************************************* //

12
tutorials/compressible/realFluidPisoFoam/ras/cavity_IAPWS97/0.org/T
View file

@ -1,9 +1,9 @@
/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: 2.1.x |
| \\ / A nd | Web: www.OpenFOAM.org |
| \\/ M anipulation | |
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | foam-extend: Open Source CFD
\\ / O peration |
\\ / A nd | For copyright notice see file Copyright
\\/ M anipulation |
\*---------------------------------------------------------------------------*/
FoamFile
{

12
tutorials/compressible/realFluidPisoFoam/ras/cavity_IAPWS97/0.org/U
View file

@ -1,9 +1,9 @@
/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: 2.1.x |
| \\ / A nd | Web: www.OpenFOAM.org |
| \\/ M anipulation | |
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | foam-extend: Open Source CFD
\\ / O peration |
\\ / A nd | For copyright notice see file Copyright
\\/ M anipulation |
\*---------------------------------------------------------------------------*/
FoamFile
{

20
tutorials/compressible/realFluidPisoFoam/ras/cavity_IAPWS97/0.org/epsilon
View file

@ -1,8 +1,8 @@
/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: 2.1.x |
| \\ / A nd | Web: www.OpenFOAM.org |
| \\ / F ield | foam-extend: Open Source CFD |
| \\ / O peration | Version: 3.1 |
| \\ / A nd | Web: http://www.extend-project.de |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/
FoamFile
@ -24,15 +24,21 @@ boundaryField
movingWall
{
type compressible::epsilonWallFunction;
value uniform 0;
refValue uniform 0;
value uniform 0.000765;
Cmu 0.09;
kappa 0.41;
E 9.8;
}
fixedWalls
{
type compressible::epsilonWallFunction;
value uniform 0;
refValue uniform 0;
value uniform 0.000765;
Cmu 0.09;
kappa 0.41;
E 9.8;
}
frontAndBack
{
type empty;

6
tutorials/compressible/realFluidPisoFoam/ras/cavity_IAPWS97/0.org/k
View file

@ -1,8 +1,8 @@
/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: 2.1.x |
| \\ / A nd | Web: www.OpenFOAM.org |
| \\ / F ield | foam-extend: Open Source CFD |
| \\ / O peration | Version: 3.1 |
| \\ / A nd | Web: http://www.extend-project.de |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/
FoamFile

12
tutorials/compressible/realFluidPisoFoam/ras/cavity_IAPWS97/0.org/p
View file

@ -1,9 +1,9 @@
/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: 2.1.x |
| \\ / A nd | Web: www.OpenFOAM.org |
| \\/ M anipulation | |
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | foam-extend: Open Source CFD
\\ / O peration |
\\ / A nd | For copyright notice see file Copyright
\\/ M anipulation |
\*---------------------------------------------------------------------------*/
FoamFile
{

12
tutorials/compressible/realFluidPisoFoam/ras/cavity_IAPWS97/0/T
View file

@ -1,9 +1,9 @@
/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: 2.1.x |
| \\ / A nd | Web: www.OpenFOAM.org |
| \\/ M anipulation | |
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | foam-extend: Open Source CFD
\\ / O peration |
\\ / A nd | For copyright notice see file Copyright
\\/ M anipulation |
\*---------------------------------------------------------------------------*/
FoamFile
{

12
tutorials/compressible/realFluidPisoFoam/ras/cavity_IAPWS97/0/U
View file

@ -1,9 +1,9 @@
/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: 2.1.x |
| \\ / A nd | Web: www.OpenFOAM.org |
| \\/ M anipulation | |
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | foam-extend: Open Source CFD
\\ / O peration |
\\ / A nd | For copyright notice see file Copyright
\\/ M anipulation |
\*---------------------------------------------------------------------------*/
FoamFile
{

22
tutorials/compressible/realFluidPisoFoam/ras/cavity_IAPWS97/0/epsilon
View file

@ -1,9 +1,9 @@
/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: 2.1.x |
| \\ / A nd | Web: www.OpenFOAM.org |
| \\/ M anipulation | |
| \\ / F ield | foam-extend: Open Source CFD |
| \\ / O peration | Version: 4.0 |
| \\ / A nd | Web: http://www.foam-extend.org |
| \\/ M anipulation | For copyright notice see file Copyright |
\*---------------------------------------------------------------------------*/
FoamFile
{
@ -24,15 +24,21 @@ boundaryField
movingWall
{
type compressible::epsilonWallFunction;
value uniform 0;
refValue uniform 0;
value uniform 0.000765;
Cmu 0.09;
kappa 0.41;
E 9.8;
}
fixedWalls
{
type compressible::epsilonWallFunction;
value uniform 0;
refValue uniform 0;
value uniform 0.000765;
Cmu 0.09;
kappa 0.41;
E 9.8;
}
frontAndBack
{
type empty;

8
tutorials/compressible/realFluidPisoFoam/ras/cavity_IAPWS97/0/k
View file

@ -1,9 +1,9 @@
/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: 2.1.x |
| \\ / A nd | Web: www.OpenFOAM.org |
| \\/ M anipulation | |
| \\ / F ield | foam-extend: Open Source CFD |
| \\ / O peration | Version: 4.0 |
| \\ / A nd | Web: http://www.foam-extend.org |
| \\/ M anipulation | For copyright notice see file Copyright |
\*---------------------------------------------------------------------------*/
FoamFile
{

12
tutorials/compressible/realFluidPisoFoam/ras/cavity_IAPWS97/0/p
View file

@ -1,9 +1,9 @@
/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: 2.1.x |
| \\ / A nd | Web: www.OpenFOAM.org |
| \\/ M anipulation | |
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | foam-extend: Open Source CFD
\\ / O peration |
\\ / A nd | For copyright notice see file Copyright
\\/ M anipulation |
\*---------------------------------------------------------------------------*/
FoamFile
{

61
tutorials/compressible/realFluidPisoFoam/ras/cavity_IAPWS97/README
View file

@ -1,61 +0,0 @@
This tutorial is to show how the IAPWSThermo work in OpenFOAM
//***********************INFORMATION TO TEST CASE*******************************//
Important files in this case:
1. thermophysicalProperties
--> see thermotype
1. fvSolution
PISO subdict --> see the realFluid --> this must be true (see modification in the pEqn of realFluidPisoFoam)
2. controlDict:
the dynamic libraries used for the IAPWS 97 water properties are not directly link into OpenFOAM.
Therefore, the dynamic libraries need to be defined in controlDict
//************************NEW CODE NEEDED FOR THIS CASE*************************//
1. Code in OpenFOAM: $WM_PROJECT_DIR/src/thermophysicalModels/externalMedia/
just run wmake libso in the folder
2. FreeSteam: download the code @ http://freesteam.sourceforge.net/example.php
get the sourceCode --> I'm not sure of the packages work with this code
see README.txt file to see what libraries you need to compile freeSteam
//***********************Problems*******************************//
If problems occur using the IAPWS 97 water properties, please check the following points
1. have you suggesfully compiled freeSteam.
a: If no, have you installed gsl, scons and so on (see README.txt in the freeSteam)
b: If yes, is the libfreesteam.so in a folder where OpenFOAM finds it e.g. $WM_PROJECT_DIR/lib/linux64GccDPOpt
(last folder name depends on your machine, compiler ...)
2. have you compiled the thermophysical models that connect freeSteam to OpenFOAM. They can be found in this folder: $WM_PROJECT_DIR/src/
thermophysicalModels/externalMedia/
3. IS THE SOLVER STARTIGN TO RUN THE CASE? NO
check b23.c file at line 39 (or at least somewhere in the code:
if you find the following code:
return B23_N[4] + sqrt((pi - B23_N[5])/B23_N[3]) /* * 1{K} */;
change it to this:
return B23_N[4] + sqrt(fabs(pi - B23_N[5])/B23_N[3]) /* * 1{K} */;

12
tutorials/compressible/realFluidPisoFoam/ras/cavity_IAPWS97/constant/RASProperties
View file

@ -1,9 +1,9 @@
/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: 1.6 |
| \\ / A nd | Web: www.OpenFOAM.org |
| \\/ M anipulation | |
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | foam-extend: Open Source CFD
\\ / O peration |
\\ / A nd | For copyright notice see file Copyright
\\/ M anipulation |
\*---------------------------------------------------------------------------*/
FoamFile
{

Some files were not shown because too many files have changed in this diff Show more