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

Removing trailing whitespace, tabs and DOS CR

Browse source
This commit is contained in:
Henrik Rusche 2016-05-30 12:23:04 +02:00
parent 40b0d8adf1
commit d16b06a15a
47 changed files with 1328 additions and 1328 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

2
src/thermophysicalModels/basic/basicThermo/basicThermo.H
View file

@ -247,7 +247,7 @@ public:
//- Thermal diffusivity for enthalpy of mixture [kg/m/s]
virtual const volScalarField& alpha() const;
//- Read thermophysicalProperties dictionary
virtual bool read();
};

4
src/thermophysicalModels/basic/psiThermo/basicPsiThermo/basicPsiThermo.H
View file

@ -109,10 +109,10 @@ public:
return p_*psi();
}
//CL: drhodh needed for pressure equation of the real gas solver
//CL: drhodh needed for pressure equation of the real gas solver
virtual const volScalarField& drhodh() const;
//CL: drhode needed for pressure equation of the real gas solver
//CL: drhode needed for pressure equation of the real gas solver
virtual const volScalarField& drhode() const;
};

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

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd |
\\ / A nd |
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -25,7 +25,7 @@ License
Author
Christian Lucas
Institut für Thermodynamik
Technische Universität Braunschweig
Technische Universität Braunschweig
Germany
\*---------------------------------------------------------------------------*/
@ -170,7 +170,7 @@ Foam::realGasEThermo<MixtureType>::realGasEThermo
forAll(rho_.boundaryField(), patchi)
{
rho_.boundaryField()[patchi] ==
rho(this->T_.boundaryField()[patchi], patchi);
rho(this->T_.boundaryField()[patchi], patchi);
}
forAll(eCells, celli)
@ -229,7 +229,7 @@ Foam::tmp<Foam::scalarField> Foam::realGasEThermo<MixtureType>::e
) const
{
//CL: need the pressure of the internal field to calculate the realGas internal energy
//CL: this is done this way to assure compatibility to old OF Thermo-Versions
//CL: this is done this way to assure compatibility to old OF Thermo-Versions
const scalarField& pCells = this->p_.internalField();
tmp<scalarField> te(new scalarField(T.size()));
@ -241,7 +241,7 @@ Foam::tmp<Foam::scalarField> Foam::realGasEThermo<MixtureType>::e
}
return te;
}
}
template<class MixtureType>
Foam::tmp<Foam::scalarField> Foam::realGasEThermo<MixtureType>::e
@ -251,7 +251,7 @@ Foam::tmp<Foam::scalarField> Foam::realGasEThermo<MixtureType>::e
) const
{
//CL: need the pressure at the patch to calculate the realGas internal energy
//CL: this is done this way to assure compatibility to old OF Thermo-Versions
//CL: this is done this way to assure compatibility to old OF Thermo-Versions
const fvPatchScalarField& pp = this->p_.boundaryField()[patchi];
tmp<scalarField> te(new scalarField(T.size()));
@ -273,7 +273,7 @@ Foam::tmp<Foam::scalarField> Foam::realGasEThermo<MixtureType>::rho
) const
{
//CL: need the pressure at the patch to calculate the realGas density
//CL: this is done this way to assure compatibility to old OF Thermo-Versions
//CL: this is done this way to assure compatibility to old OF Thermo-Versions
const fvPatchScalarField& pp = this->p_.boundaryField()[patchi];
tmp<scalarField> trho(new scalarField(T.size()));
@ -295,7 +295,7 @@ Foam::tmp<Foam::scalarField> Foam::realGasEThermo<MixtureType>::Cp
) const
{
//CL: need the pressure at the patch to calculate the realGas cp
//CL: this is done this way to assure compatibility to old OF Thermo-Versions
//CL: this is done this way to assure compatibility to old OF Thermo-Versions
const fvPatchScalarField& pp = this->p_.boundaryField()[patchi];
tmp<scalarField> tCp(new scalarField(T.size()));
@ -362,7 +362,7 @@ Foam::tmp<Foam::scalarField> Foam::realGasEThermo<MixtureType>::Cv
) const
{
//CL: need the pressure at the patch to calculate the realGas internal energy
//CL: this is done this way to assure compatibility to old OF Thermo-Versions
//CL: this is done this way to assure compatibility to old OF Thermo-Versions
const fvPatchScalarField& pp = this->p_.boundaryField()[patchi];
tmp<scalarField> tCv(new scalarField(T.size()));
@ -384,7 +384,7 @@ Foam::tmp<Foam::volScalarField> Foam::realGasEThermo<MixtureType>::rho() const
const fvMesh& mesh = this->T_.mesh();
//CL: create an rho Field, which will be return
//CL: the problem is that this function is "const",
//CL: the problem is that this function is "const",
//CL: so a new variabel is needed
tmp<volScalarField> trho
(
@ -402,11 +402,11 @@ Foam::tmp<Foam::volScalarField> Foam::realGasEThermo<MixtureType>::rho() const
dimDensity
)
);
//CL: copy "old" rho value onto the new rho field as start point
//CL: for the newton solver used in this->TE( ... )
//CL: copy "old" rho value onto the new rho field as start point
//CL: for the newton solver used in this->TE( ... )
trho()=rho_;
volScalarField& rho = trho();
const scalarField& eCells = e_.internalField();
@ -426,8 +426,8 @@ Foam::tmp<Foam::volScalarField> Foam::realGasEThermo<MixtureType>::rho() const
{
fvPatchScalarField pp = this->p_.boundaryField()[patchi];
fvPatchScalarField pe = e_.boundaryField()[patchi];
fvPatchScalarField pT = this->T_.boundaryField()[patchi];
fvPatchScalarField pT = this->T_.boundaryField()[patchi];
fvPatchScalarField& prho_ = rho.boundaryField()[patchi];
forAll(pp, facei)

4
src/thermophysicalModels/basic/psiThermo/realGasEThermo/realGasEThermo.H
View file

@ -35,7 +35,7 @@ SourceFiles
Author
Christian Lucas
Institut für Thermodynamik
Technische Universität Braunschweig
Technische Universität Braunschweig
Germany
@ -130,7 +130,7 @@ public:
return e_;
}
//CL: drhode needed for pressure equation of the real gas solver
//CL: drhode needed for pressure equation of the real gas solver
virtual const volScalarField& drhode() const
{
return drhode_;

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

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd |
\\ / A nd |
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -25,7 +25,7 @@ License
Author
Christian Lucas
Institut für Thermodynamik
Technische Universität Braunschweig
Technische Universität Braunschweig
Germany

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

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd |
\\ / A nd |
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -25,7 +25,7 @@ License
Author
Christian Lucas
Institut für Thermodynamik
Technische Universität Braunschweig
Technische Universität Braunschweig
Germany
\*---------------------------------------------------------------------------*/
@ -169,12 +169,12 @@ Foam::realGasHThermo<MixtureType>::realGasHThermo
forAll(rho_.boundaryField(), patchi)
{
rho_.boundaryField()[patchi] ==
rho(this->T_.boundaryField()[patchi], patchi);
rho(this->T_.boundaryField()[patchi], patchi);
}
forAll(hCells, celli)
{
hCells[celli] = this->cellMixture(celli).H(rhoCells[celli],TCells[celli]);
hCells[celli] = this->cellMixture(celli).H(rhoCells[celli],TCells[celli]);
}
forAll(h_.boundaryField(), patchi)
@ -228,7 +228,7 @@ Foam::tmp<Foam::scalarField> Foam::realGasHThermo<MixtureType>::h
) const
{
//CL: need the pressure of the internal field to calculate the realGas enthalpy
//CL: this is done this way to assure compatibility to old OF Thermo-Versions
//CL: this is done this way to assure compatibility to old OF Thermo-Versions
const scalarField& pCells = this->p_.internalField();
tmp<scalarField> th(new scalarField(T.size()));
@ -240,7 +240,7 @@ Foam::tmp<Foam::scalarField> Foam::realGasHThermo<MixtureType>::h
}
return th;
}
}
template<class MixtureType>
Foam::tmp<Foam::scalarField> Foam::realGasHThermo<MixtureType>::h
@ -250,7 +250,7 @@ Foam::tmp<Foam::scalarField> Foam::realGasHThermo<MixtureType>::h
) const
{
//CL: need the pressure at the patch to calculate the realGas enthalpy
//CL: this is done this way to assure compatibility to old OF Thermo-Versions
//CL: this is done this way to assure compatibility to old OF Thermo-Versions
const fvPatchScalarField& pp = this->p_.boundaryField()[patchi];
tmp<scalarField> th(new scalarField(T.size()));
@ -272,7 +272,7 @@ Foam::tmp<Foam::scalarField> Foam::realGasHThermo<MixtureType>::rho
) const
{
//CL: need the pressure at the patch to calculate the realGas enthalpy
//CL: this is done this way to assure compatibility to old OF Thermo-Versions
//CL: this is done this way to assure compatibility to old OF Thermo-Versions
const fvPatchScalarField& pp = this->p_.boundaryField()[patchi];
tmp<scalarField> trho(new scalarField(T.size()));
@ -294,7 +294,7 @@ Foam::tmp<Foam::scalarField> Foam::realGasHThermo<MixtureType>::Cp
) const
{
//CL: need the pressure at the patch to calculate the realGas enthalpy
//CL: this is done this way to assure compatibility to old OF Thermo-Versions
//CL: this is done this way to assure compatibility to old OF Thermo-Versions
const fvPatchScalarField& pp = this->p_.boundaryField()[patchi];
tmp<scalarField> tCp(new scalarField(T.size()));
@ -361,7 +361,7 @@ Foam::tmp<Foam::scalarField> Foam::realGasHThermo<MixtureType>::Cv
) const
{
//CL: need the pressure at the patch to calculate the realGas enthalpy
//CL: this is done this way to assure compatibility to old OF Thermo-Versions
//CL: this is done this way to assure compatibility to old OF Thermo-Versions
const fvPatchScalarField& pp = this->p_.boundaryField()[patchi];
tmp<scalarField> tCv(new scalarField(T.size()));
@ -383,7 +383,7 @@ Foam::tmp<Foam::volScalarField> Foam::realGasHThermo<MixtureType>::rho() const
const fvMesh& mesh = this->T_.mesh();
//CL: create a rho Field, which will be return
//CL: the problem is that this function is "const",
//CL: the problem is that this function is "const",
//CL: so a new variabel is needed
tmp<volScalarField> trho
(
@ -401,11 +401,11 @@ Foam::tmp<Foam::volScalarField> Foam::realGasHThermo<MixtureType>::rho() const
dimDensity
)
);
//CL: copy "old" rho value onto the new rho field as start point
//CL: for the newton solver used in this->TH( ... )
//CL: copy "old" rho value onto the new rho field as start point
//CL: for the newton solver used in this->TH( ... )
trho()=rho_;
volScalarField& rho = trho();
const scalarField& hCells = h_.internalField();
@ -425,8 +425,8 @@ Foam::tmp<Foam::volScalarField> Foam::realGasHThermo<MixtureType>::rho() const
{
fvPatchScalarField pp = this->p_.boundaryField()[patchi];
fvPatchScalarField ph = h_.boundaryField()[patchi];
fvPatchScalarField pT = this->T_.boundaryField()[patchi];
fvPatchScalarField pT = this->T_.boundaryField()[patchi];
fvPatchScalarField& prho_ = rho.boundaryField()[patchi];
forAll(pp, facei)

4
src/thermophysicalModels/basic/psiThermo/realGasHThermo/realGasHThermo.H
View file

@ -35,7 +35,7 @@ SourceFiles
Author
Christian Lucas
Institut für Thermodynamik
Technische Universität Braunschweig
Technische Universität Braunschweig
Germany
@ -130,7 +130,7 @@ public:
return h_;
}
//CL: drhodh needed for pressure equation of the real gas solver
//CL: drhodh needed for pressure equation of the real gas solver
virtual const volScalarField& drhodh() const
{
return drhodh_;

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

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd |
\\ / A nd |
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -25,7 +25,7 @@ License
Author
Christian Lucas
Institut für Thermodynamik
Technische Universität Braunschweig
Technische Universität Braunschweig
Germany

1178
src/thermophysicalModels/externalMedia/IAPWS_Waterproperties/IAPWSThermo/IAPWS-IF97.C
View file

File diff suppressed because it is too large Load diff

462
src/thermophysicalModels/externalMedia/IAPWS_Waterproperties/IAPWSThermo/IAPWS-IF97.H
View file

@ -1,231 +1,231 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / 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
IAPWS-IF97 (water) based thermodynamic class. Water properties calculated by freeSteam.
This code connects OpenFoam with freeSteam and provides the basic functions needed in OpenFOAM
For more information about freeSteam and its authors have a look @ http://freesteam.sourceforge.net/example.php
SourceFiles
IAPWS-IF97.C
Author
Christian Lucas
Institut für Thermodynamik
Technische Universität Braunschweig
Germany
\*---------------------------------------------------------------------------*/
#ifndef IAPWSIF97_H
#define IAPWSIF97_H
#include "basicPsiThermo.H"
#include "steam.H"
#ifdef __cplusplus
#define EXTERN extern "C"
#else
#define EXTERN extern
#endif
//CL:
EXTERN double freesteam_p(SteamState S);
EXTERN double freesteam_T(SteamState S);
EXTERN double freesteam_rho(SteamState S);
EXTERN double freesteam_v(SteamState S);
EXTERN double freesteam_u(SteamState S);
EXTERN double freesteam_h(SteamState S);
EXTERN double freesteam_s(SteamState S);
EXTERN double freesteam_cp(SteamState S);
EXTERN double freesteam_cv(SteamState S);
EXTERN double freesteam_w(SteamState S);
EXTERN double freesteam_x(SteamState S);
EXTERN double freesteam_mu(SteamState S);
EXTERN double freesteam_k(SteamState S);
//CL: getting SteamState for two given properties e.g. pressure and temperatur
EXTERN SteamState freesteam_set_pv(double,double);
EXTERN SteamState freesteam_set_pu(double,double);
EXTERN SteamState freesteam_set_pT(double,double);
EXTERN SteamState freesteam_set_ph(double,double);
//CL: getting region of the SteamState
EXTERN int freesteam_region(SteamState);
//CL: transport properties
EXTERN double freesteam_mu_rhoT(double,double);
EXTERN double freesteam_k_rhoT(double,double);
//CL: Region 1 --> see region1.h (freesteam)
EXTERN double freesteam_region1_v_pT(double,double);
EXTERN double freesteam_region1_h_pT(double,double);
EXTERN double freesteam_region1_kappaT_pT(double,double);
EXTERN double freesteam_region1_alphav_pT(double,double);
EXTERN double freesteam_region1_cp_pT(double,double);
EXTERN double freesteam_region1_u_pT(double,double);
EXTERN double freesteam_region1_s_pT(double,double);
EXTERN double freesteam_region1_cv_pT(double,double);
//CL: Region 2 --> see region2.h (freesteam)
EXTERN double freesteam_region2_v_pT(double,double);
EXTERN double freesteam_region2_u_pT(double,double);
EXTERN double freesteam_region2_s_pT(double,double);
EXTERN double freesteam_region2_h_pT(double,double);
EXTERN double freesteam_region2_cp_pT(double,double);
EXTERN double freesteam_region2_cv_pT(double,double);
EXTERN double freesteam_region2_alphav_pT(double,double);
EXTERN double freesteam_region2_kappaT_pT(double,double);
//CL: Region 3 --> see region3.h (freesteam)
EXTERN double freesteam_region3_p_rhoT(double,double);
EXTERN double freesteam_region3_u_rhoT(double,double);
EXTERN double freesteam_region3_s_rhoT(double,double);
EXTERN double freesteam_region3_h_rhoT(double,double);
EXTERN double freesteam_region3_cp_rhoT(double,double);
EXTERN double freesteam_region3_cv_rhoT(double,double);
EXTERN double freesteam_region3_alphap_rhoT(double,double);
EXTERN double freesteam_region3_betap_rhoT(double,double);
//CL: Region 4 --> see region4.h (freesteam)
EXTERN double freesteam_region4_psat_T(double);
EXTERN double freesteam_region4_Tsat_p(double);
EXTERN double freesteam_region4_rhof_T(double);
EXTERN double freesteam_region4_rhog_T(double);
EXTERN double freesteam_region4_v_Tx(double,double);
EXTERN double freesteam_region4_u_Tx(double,double);
EXTERN double freesteam_region4_h_Tx(double,double);
EXTERN double freesteam_region4_s_Tx(double,double);
EXTERN double freesteam_region4_cp_Tx(double,double);
EXTERN double freesteam_region4_cv_Tx(double,double);
EXTERN double freesteam_region4_dpsatdT_T(double);
namespace Foam
{
//CL: Functions to caluculate all fluid properties
void calculateProperties_h
(
SteamState S,
scalar &rho,
scalar &h,
scalar &T,
scalar &p,
scalar &psi,
scalar &drhodh,
scalar &mu,
scalar &alpha,
scalar &x
);
//CL: This functions returns all (minimal) needed propeties (p,T,h,rho,psi,drhodh,mu and alpha) for given p and T
void calculateProperties_pT
(
scalar &p,
scalar &T,
scalar &h,
scalar &rho,
scalar &psi,
scalar &drhodh,
scalar &mu,
scalar &alpha
);
//CL: This function returns the same values as the function above for given p and T
//CL: Additionally, the vapor mass fraction x is return
//CL: NOTE: This function is only included to have the possibility to update x at the fixedValue (Temperature) BC
//CL: can only return x=0 and x=1 because it is not possible to describe the vapour dome with p and T
void calculateProperties_pT
(
scalar &p,
scalar &T,
scalar &h,
scalar &rho,
scalar &psi,
scalar &drhodh,
scalar &mu,
scalar &alpha,
scalar &x
);
//CL: This functions returns all (minimal) needed properties (p,T,h,rho,psi,drhodh,mu and alpha) for given p and h
void calculateProperties_ph
(
scalar &p,
scalar &h,
scalar &T,
scalar &rho,
scalar &psi,
scalar &drhodh,
scalar &mu,
scalar &alpha
);
//CL: This function returns the same values as the function above for given p and h
//CL: Additionally, the vapor mass fraction x is return
void calculateProperties_ph
(
scalar &p,
scalar &h,
scalar &T,
scalar &rho,
scalar &psi,
scalar &drhodh,
scalar &mu,
scalar &alpha,
scalar &x
);
//CL: Return density for given pT or ph;
scalar rho_pT(scalar p,scalar T);
scalar rho_ph(scalar p,scalar h);
//CL: Return cp for given pT or ph;
scalar cp_pT(scalar p,scalar T);
scalar cp_ph(scalar p,scalar h);
//CL: Return cv for given pT or ph;
scalar cv_pT(scalar p,scalar T);
scalar cv_ph(scalar p,scalar h);
//CL: Return enthalpy for given pT;
scalar h_pT(scalar p,scalar T);
//CL: Return temperature for given ph;
scalar T_ph(scalar p,scalar T);
//CL: Return psiH=(drho/dp)_h=constant for given pT or ph;
scalar psiH_pT(scalar p,scalar T);
scalar psiH_ph(scalar p,scalar h);
scalar psiH(SteamState S);
//CL: Return drhodh=(drho/dh)_p=constant for given pT or ph;
scalar drhodh_pT(scalar p,scalar T);
scalar drhodh_ph(scalar p,scalar h);
scalar drhodh(SteamState S);
}
#endif //IAPWSIF97_C_
/*---------------------------------------------------------------------------*\
========= |
\\ / 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
IAPWS-IF97 (water) based thermodynamic class. Water properties calculated by freeSteam.
This code connects OpenFoam with freeSteam and provides the basic functions needed in OpenFOAM
For more information about freeSteam and its authors have a look @ http://freesteam.sourceforge.net/example.php
SourceFiles
IAPWS-IF97.C
Author
Christian Lucas
Institut für Thermodynamik
Technische Universität Braunschweig
Germany
\*---------------------------------------------------------------------------*/
#ifndef IAPWSIF97_H
#define IAPWSIF97_H
#include "basicPsiThermo.H"
#include "steam.H"
#ifdef __cplusplus
#define EXTERN extern "C"
#else
#define EXTERN extern
#endif
//CL:
EXTERN double freesteam_p(SteamState S);
EXTERN double freesteam_T(SteamState S);
EXTERN double freesteam_rho(SteamState S);
EXTERN double freesteam_v(SteamState S);
EXTERN double freesteam_u(SteamState S);
EXTERN double freesteam_h(SteamState S);
EXTERN double freesteam_s(SteamState S);
EXTERN double freesteam_cp(SteamState S);
EXTERN double freesteam_cv(SteamState S);
EXTERN double freesteam_w(SteamState S);
EXTERN double freesteam_x(SteamState S);
EXTERN double freesteam_mu(SteamState S);
EXTERN double freesteam_k(SteamState S);
//CL: getting SteamState for two given properties e.g. pressure and temperatur
EXTERN SteamState freesteam_set_pv(double,double);
EXTERN SteamState freesteam_set_pu(double,double);
EXTERN SteamState freesteam_set_pT(double,double);
EXTERN SteamState freesteam_set_ph(double,double);
//CL: getting region of the SteamState
EXTERN int freesteam_region(SteamState);
//CL: transport properties
EXTERN double freesteam_mu_rhoT(double,double);
EXTERN double freesteam_k_rhoT(double,double);
//CL: Region 1 --> see region1.h (freesteam)
EXTERN double freesteam_region1_v_pT(double,double);
EXTERN double freesteam_region1_h_pT(double,double);
EXTERN double freesteam_region1_kappaT_pT(double,double);
EXTERN double freesteam_region1_alphav_pT(double,double);
EXTERN double freesteam_region1_cp_pT(double,double);
EXTERN double freesteam_region1_u_pT(double,double);
EXTERN double freesteam_region1_s_pT(double,double);
EXTERN double freesteam_region1_cv_pT(double,double);
//CL: Region 2 --> see region2.h (freesteam)
EXTERN double freesteam_region2_v_pT(double,double);
EXTERN double freesteam_region2_u_pT(double,double);
EXTERN double freesteam_region2_s_pT(double,double);
EXTERN double freesteam_region2_h_pT(double,double);
EXTERN double freesteam_region2_cp_pT(double,double);
EXTERN double freesteam_region2_cv_pT(double,double);
EXTERN double freesteam_region2_alphav_pT(double,double);
EXTERN double freesteam_region2_kappaT_pT(double,double);
//CL: Region 3 --> see region3.h (freesteam)
EXTERN double freesteam_region3_p_rhoT(double,double);
EXTERN double freesteam_region3_u_rhoT(double,double);
EXTERN double freesteam_region3_s_rhoT(double,double);
EXTERN double freesteam_region3_h_rhoT(double,double);
EXTERN double freesteam_region3_cp_rhoT(double,double);
EXTERN double freesteam_region3_cv_rhoT(double,double);
EXTERN double freesteam_region3_alphap_rhoT(double,double);
EXTERN double freesteam_region3_betap_rhoT(double,double);
//CL: Region 4 --> see region4.h (freesteam)
EXTERN double freesteam_region4_psat_T(double);
EXTERN double freesteam_region4_Tsat_p(double);
EXTERN double freesteam_region4_rhof_T(double);
EXTERN double freesteam_region4_rhog_T(double);
EXTERN double freesteam_region4_v_Tx(double,double);
EXTERN double freesteam_region4_u_Tx(double,double);
EXTERN double freesteam_region4_h_Tx(double,double);
EXTERN double freesteam_region4_s_Tx(double,double);
EXTERN double freesteam_region4_cp_Tx(double,double);
EXTERN double freesteam_region4_cv_Tx(double,double);
EXTERN double freesteam_region4_dpsatdT_T(double);
namespace Foam
{
//CL: Functions to caluculate all fluid properties
void calculateProperties_h
(
SteamState S,
scalar &rho,
scalar &h,
scalar &T,
scalar &p,
scalar &psi,
scalar &drhodh,
scalar &mu,
scalar &alpha,
scalar &x
);
//CL: This functions returns all (minimal) needed propeties (p,T,h,rho,psi,drhodh,mu and alpha) for given p and T
void calculateProperties_pT
(
scalar &p,
scalar &T,
scalar &h,
scalar &rho,
scalar &psi,
scalar &drhodh,
scalar &mu,
scalar &alpha
);
//CL: This function returns the same values as the function above for given p and T
//CL: Additionally, the vapor mass fraction x is return
//CL: NOTE: This function is only included to have the possibility to update x at the fixedValue (Temperature) BC
//CL: can only return x=0 and x=1 because it is not possible to describe the vapour dome with p and T
void calculateProperties_pT
(
scalar &p,
scalar &T,
scalar &h,
scalar &rho,
scalar &psi,
scalar &drhodh,
scalar &mu,
scalar &alpha,
scalar &x
);
//CL: This functions returns all (minimal) needed properties (p,T,h,rho,psi,drhodh,mu and alpha) for given p and h
void calculateProperties_ph
(
scalar &p,
scalar &h,
scalar &T,
scalar &rho,
scalar &psi,
scalar &drhodh,
scalar &mu,
scalar &alpha
);
//CL: This function returns the same values as the function above for given p and h
//CL: Additionally, the vapor mass fraction x is return
void calculateProperties_ph
(
scalar &p,
scalar &h,
scalar &T,
scalar &rho,
scalar &psi,
scalar &drhodh,
scalar &mu,
scalar &alpha,
scalar &x
);
//CL: Return density for given pT or ph;
scalar rho_pT(scalar p,scalar T);
scalar rho_ph(scalar p,scalar h);
//CL: Return cp for given pT or ph;
scalar cp_pT(scalar p,scalar T);
scalar cp_ph(scalar p,scalar h);
//CL: Return cv for given pT or ph;
scalar cv_pT(scalar p,scalar T);
scalar cv_ph(scalar p,scalar h);
//CL: Return enthalpy for given pT;
scalar h_pT(scalar p,scalar T);
//CL: Return temperature for given ph;
scalar T_ph(scalar p,scalar T);
//CL: Return psiH=(drho/dp)_h=constant for given pT or ph;
scalar psiH_pT(scalar p,scalar T);
scalar psiH_ph(scalar p,scalar h);
scalar psiH(SteamState S);
//CL: Return drhodh=(drho/dh)_p=constant for given pT or ph;
scalar drhodh_pT(scalar p,scalar T);
scalar drhodh_ph(scalar p,scalar h);
scalar drhodh(SteamState S);
}
#endif //IAPWSIF97_C_

62
src/thermophysicalModels/externalMedia/IAPWS_Waterproperties/IAPWSThermo/IAPWSThermo.C
View file

@ -25,7 +25,7 @@ License
Author
Christian Lucas
Institut für Thermodynamik
Technische Universität Braunschweig
Technische Universität Braunschweig
Germany
\*---------------------------------------------------------------------------*/
@ -45,21 +45,21 @@ void Foam::IAPWSThermo::calculate()
scalarField& drhodhCells = this->drhodh_.internalField();
scalarField& muCells = this->mu_.internalField();
scalarField& alphaCells = this->alpha_.internalField();
//CL: Updating all cell properties
//CL: loop through all cells
forAll(TCells, celli)
{
{
//CL: see IAPWAS-IF97.H
calculateProperties_ph
(
pCells[celli],
hCells[celli],
TCells[celli],
rhoCells[celli],
psiCells[celli],
drhodhCells[celli],
muCells[celli],
pCells[celli],
hCells[celli],
TCells[celli],
rhoCells[celli],
psiCells[celli],
drhodhCells[celli],
muCells[celli],
alphaCells[celli]
);
}
@ -84,13 +84,13 @@ void Foam::IAPWSThermo::calculate()
//CL: see IAPWAS-IF97.H
calculateProperties_pT
(
pp[facei],
pT[facei],
ph[facei],
prho[facei],
ppsi[facei],
pdrhodh[facei],
pmu[facei],
pp[facei],
pT[facei],
ph[facei],
prho[facei],
ppsi[facei],
pdrhodh[facei],
pmu[facei],
palpha[facei]
);
}
@ -103,13 +103,13 @@ void Foam::IAPWSThermo::calculate()
//CL: see IAPWAS-IF97.H
calculateProperties_ph
(
pp[facei],
ph[facei],
pT[facei],
prho[facei],
ppsi[facei],
pdrhodh[facei],
pmu[facei],
pp[facei],
ph[facei],
pT[facei],
prho[facei],
ppsi[facei],
pdrhodh[facei],
pmu[facei],
palpha[facei]
);
}
@ -143,7 +143,7 @@ Foam::IAPWSThermo::IAPWSThermo
dimensionSet(0, 2, -2, 0, 0),
this->hBoundaryTypes()
),
rho_
(
IOobject
@ -186,18 +186,18 @@ Foam::IAPWSThermo::IAPWSThermo
forAll(h_.boundaryField(), patchi)
{
h_.boundaryField()[patchi] ==
h(this->T_.boundaryField()[patchi], patchi);
h(this->T_.boundaryField()[patchi], patchi);
}
forAll(rhoCells, celli)
{
rhoCells[celli] = rho_pT(pCells[celli],TCells[celli]);
rhoCells[celli] = rho_pT(pCells[celli],TCells[celli]);
}
forAll(rho_.boundaryField(), patchi)
{
rho_.boundaryField()[patchi] ==
rho(this->p_.boundaryField()[patchi] ,this->h_.boundaryField()[patchi], patchi);
rho(this->p_.boundaryField()[patchi] ,this->h_.boundaryField()[patchi], patchi);
}
hBoundaryCorrection(h_);
@ -256,7 +256,7 @@ Foam::tmp<Foam::scalarField> Foam::IAPWSThermo::h
}
return th;
}
}
Foam::tmp<Foam::scalarField> Foam::IAPWSThermo::h
@ -273,7 +273,7 @@ Foam::tmp<Foam::scalarField> Foam::IAPWSThermo::h
forAll(T, facei)
{
h[facei] = h_pT(pp[facei], T[facei]);
h[facei] = h_pT(pp[facei], T[facei]);
}
return th;
@ -365,7 +365,7 @@ Foam::tmp<Foam::volScalarField> Foam::IAPWSThermo::Cp() const
}
//CL: Returns an updated field for rho
//CL: Returns an updated field for rho
Foam::tmp<Foam::volScalarField> Foam::IAPWSThermo::rho() const
{
const fvMesh& mesh = this->p_.mesh();

10
src/thermophysicalModels/externalMedia/IAPWS_Waterproperties/IAPWSThermo/IAPWSThermo.H
View file

@ -27,12 +27,12 @@ Class
Description:
Waterproperties based on the IAPWS 97 tables
The water properties are caluclated using freeSteam (http://freesteam.sourceforge.net/example.php)
Waterproperties based on the IAPWS 97 tables
The water properties are caluclated using freeSteam (http://freesteam.sourceforge.net/example.php)
General paper decribing the water tables:
"Revised Release on the IAPWS Industrial Formulation 1997 for the Thermodynamic Properties of Water and Steam"
"Revised Release on the IAPWS Industrial Formulation 1997 for the Thermodynamic Properties of Water and Steam"
SourceFiles
IAPWSThermo.C
@ -40,7 +40,7 @@ SourceFiles
Author
Christian Lucas
Institut für Thermodynamik
Technische Universität Braunschweig
Technische Universität Braunschweig
Germany

4
src/thermophysicalModels/externalMedia/IAPWS_Waterproperties/IAPWSThermo/IAPWSThermos.C
View file

@ -25,7 +25,7 @@ License
Author
Christian Lucas
Institut für Thermodynamik
Technische Universität Braunschweig
Technische Universität Braunschweig
Germany
\*---------------------------------------------------------------------------*/
@ -34,7 +34,7 @@ Germany
#include "makeBasicPsiThermo.H"
#include "IAPWSThermo.H"
// including dummy classes --> this classes do nothing
// including dummy classes --> this classes do nothing
// except satisfy the template structure
//#include "dummyEqnOfState.H"
//#include "dummyThermo.H"

2
src/thermophysicalModels/externalMedia/IAPWS_Waterproperties/Make/files
View file

@ -1,4 +1,4 @@
IAPWSThermo/IAPWS-IF97.C
IAPWSThermo/IAPWSThermos.C
IAPWSThermo/IAPWSThermos.C
LIB = $(FOAM_LIBBIN)/libIAPWSThermo

2
src/thermophysicalModels/externalMedia/IAPWS_Waterproperties/Make/options
View file

@ -1,7 +1,7 @@
EXE_INC = \
-I$(LIB_SRC)/finiteVolume/lnInclude \
-I$(LIB_SRC)/thermophysicalModels/specie/lnInclude \
-I$(LIB_SRC)/thermophysicalModels/basic/lnInclude
-I$(LIB_SRC)/thermophysicalModels/basic/lnInclude
LIB_LIBS = \
-lfiniteVolume\

2
src/thermophysicalModels/externalMedia/README
View file

@ -1,4 +1,4 @@
CL:
CL:
In this folder, thermodynamic classes similar to "hPsiThermo" are provided. These classes use external fluid property libraries e.g. freesteam for water.The classes are not directly linked into the basicThermophysicalModels dynamic library (see folder OpenFOAM-x.x/src/thermophysicalModels/basic). To use the code, add the dynamic library (the one you want to use from this folder + the one from the external library) in the controldict of your case

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

@ -28,7 +28,7 @@ Description
Author
Christian Lucas
Institut für Thermodynamik
Technische Universität Braunschweig
Technische Universität Braunschweig
Germany
\*---------------------------------------------------------------------------*/
@ -49,7 +49,7 @@ aungierRedlichKwong::aungierRedlichKwong(Istream& is)
pcrit_(readScalar(is)),
Tcrit_(readScalar(is)),
azentricFactor_(readScalar(is)),
rhocrit_(readScalar(is)),
rhocrit_(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_),
@ -88,7 +88,7 @@ aungierRedlichKwong::aungierRedlichKwong(const dictionary& dict)
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),
TSave(0.0),
b2_(pow(b_,2)),
b3_(pow(b_,3)),
b4_(pow(b_,4)),
@ -119,7 +119,7 @@ void Foam::aungierRedlichKwong::write(Ostream& os) const
Ostream& operator<<(Ostream& os, const aungierRedlichKwong& ark)
{
os << static_cast<const specie&>(ark)<< token::SPACE
os << static_cast<const specie&>(ark)<< token::SPACE
<< ark.pcrit_ << tab<< ark.Tcrit_<< tab<<ark.azentricFactor_<< tab<<ark.rhocrit_;
os.check("Ostream& operator<<(Ostream& os, const aungierRedlichKwong& st)");

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

@ -40,7 +40,7 @@ SourceFiles
Author
Christian Lucas
Institut für Thermodynamik
Technische Universität Braunschweig
Technische Universität Braunschweig
Germany
\*---------------------------------------------------------------------------*/
@ -70,7 +70,7 @@ protected:
// Private data
scalar pcrit_;
scalar Tcrit_;
scalar Tcrit_;
scalar azentricFactor_;
scalar rhocrit_;
@ -105,10 +105,10 @@ protected:
//Protected functions
//CL: function updates the coefficients (aSave, daSave, d2aSave)
inline void updateModelCoefficients(const scalar T) const;
inline void updateModelCoefficients(const scalar T) const;
public:
// Constructors
@ -170,46 +170,46 @@ public:
//CL: Equation of state
inline scalar p(const scalar rho, const scalar T) const;
//CL: 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;
inline scalar dvdT(const scalar rho, const scalar T) const;
inline scalar dvdp(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 isobarExpCoef(const scalar rho, const scalar T) const;
inline scalar isothermalCompressiblity
(
const scalar rho,
const scalar T
) const;
) const;
//CL: Used for cv
inline scalar integral_d2pdT2_dv
(
const scalar rho,
const scalar T
) const ;
) const ;
//CL: second order derivatives, not Used At The Moment
inline scalar d2pdv2(const scalar rho, const scalar T) const;
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, scalar T) const;
inline scalar d2pdvdT(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;
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;
//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]
//- Return density [kg/m^3]
// rho0 is the starting point of the newton solver used to calculate rho
inline scalar rho
(
@ -217,8 +217,8 @@ public:
const scalar T,
const scalar rho0
) const;
inline scalar rho(const scalar p, const scalar T) 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;
@ -227,7 +227,7 @@ public:
inline scalar Z
(
const scalar p,
const scalar T,
const scalar T,
const scalar rho0
) const;

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

@ -25,7 +25,7 @@ License
Author
Christian Lucas
Institut für Thermodynamik
Technische Universität Braunschweig
Technische Universität Braunschweig
Germany
\*---------------------------------------------------------------------------*/
@ -56,7 +56,7 @@ inline aungierRedlichKwong::aungierRedlichKwong(const word& name, const aungierR
pcrit_(pg.pcrit_),
Tcrit_(pg.Tcrit_),
azentricFactor_(pg.azentricFactor_),
rhocrit_(pg.rhocrit_),
rhocrit_(pg.rhocrit_),
a0_(pg.a0_),
b_(pg.b_),
c_(pg.c_),
@ -108,15 +108,15 @@ inline scalar aungierRedlichKwong::a(const scalar T)const
{
//CL: check if a has already been calculated for this temperature
if(TSave==T)
{
return aSave;
{
return aSave;
}
//CL: If not, recalculate a(T), dadT(T) and d2adT2(T)
else
{
updateModelCoefficients(T);
return aSave;
}
}
}
@ -126,14 +126,14 @@ inline scalar aungierRedlichKwong::dadT(const scalar T)const
// check if a has already been calculated for this temperature
if(TSave==T)
{
return daSave;
return daSave;
}
//CL: If not, recalculate a(T), dadT(T) and d2adT2(T)
else
{
updateModelCoefficients(T);
return daSave;
}
}
}
@ -143,14 +143,14 @@ inline scalar aungierRedlichKwong::d2adT2(const scalar T)const
// check if a has already been calculated for this temperature
if(TSave==T)
{
return d2aSave;
return d2aSave;
}
//CL: If not, recalculate a(T), dadT(T) and d2adT2(T)
else
{
updateModelCoefficients(T);
return d2aSave;
}
}
}
//Aungier Redlich Kwong factors
@ -202,7 +202,7 @@ inline scalar aungierRedlichKwong::c2()const
}
//returns the pressure for a given density and temperature
//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;
@ -211,7 +211,7 @@ inline scalar aungierRedlichKwong::p(const scalar rho,const scalar T) const
}
//Real deviative dp/dv at constant temperature
//Real deviative dp/dv at constant temperature
//(molar values)
inline scalar aungierRedlichKwong::dpdv(const scalar rho, const scalar T) const
{
@ -224,7 +224,7 @@ inline scalar aungierRedlichKwong::dpdv(const scalar rho, const scalar T) const
}
//Real deviative dp/dT at constant molar volume
//Real deviative dp/dT at constant molar volume
//(molar values)
inline scalar aungierRedlichKwong::dpdT(const scalar rho, const scalar T) const
{
@ -233,7 +233,7 @@ inline scalar aungierRedlichKwong::dpdT(const scalar rho, const scalar T) const
}
//Real deviative dv/dT at constant pressure
//Real deviative dv/dT at constant pressure
//using implicit differentiation
//(molar values)
inline scalar aungierRedlichKwong::dvdT(const scalar rho,const scalar T) const
@ -242,9 +242,9 @@ inline scalar aungierRedlichKwong::dvdT(const scalar rho,const scalar T) const
}
//Real deviative dv/dp at constant temperature
//Real deviative dv/dp at constant temperature
//(molar values)
inline scalar aungierRedlichKwong::dvdp(const scalar rho,const scalar T) const
inline scalar aungierRedlichKwong::dvdp(const scalar rho,const scalar T) const
{
return 1/this->dpdv(rho,T);
}
@ -266,7 +266,7 @@ inline scalar aungierRedlichKwong::integral_p_dv
//needed to calculate the entropy
//(molar values)
inline scalar aungierRedlichKwong::integral_dpdT_dv
(
(
const scalar rho,
const scalar T
) const
@ -300,10 +300,10 @@ inline scalar aungierRedlichKwong::d2pdv2(const scalar rho,const scalar T) const
//(molar values)
//using second order implicit differentiation
//using second order implicit differentiation
inline scalar aungierRedlichKwong::d2vdT2(const scalar rho, const scalar T) const
{
return
return
-(
pow(this->dpdT(rho,T),2)*this->d2pdv2(rho,T)
+ pow(this->dpdv(rho,T),2)*this->d2pdT2(rho,T)
@ -325,13 +325,13 @@ inline scalar aungierRedlichKwong::d2pdvdT(const scalar rho, const scalar T) con
}
// the result of this intergal is needed for the nasa based cp polynomial
// the result of this intergal is needed for the nasa based cp polynomial
//(molar values)
inline scalar aungierRedlichKwong::integral_d2pdT2_dv
(
const scalar rho,
const scalar T
) const
) const
{
scalar Vm = this->W()/rho;
return d2adT2(T)*log(b()+Vm)/b()-d2adT2(T)*log(Vm)/b();
@ -376,7 +376,7 @@ inline scalar aungierRedlichKwong::rho
do
{
molarVolumePrevIteration= molarVolume;
label i=0;
do
{
@ -387,7 +387,7 @@ inline scalar aungierRedlichKwong::rho
)/pow(2,i);
i++;
if (i>8)
if (i>8)
{
//CL: using bisection methode as backup,
//CL: solution must be between rho=0.001 to rho=1500;
@ -402,7 +402,7 @@ inline scalar aungierRedlichKwong::rho
if ((f2 < 0 && f3 > 0) || (f2 > 0 && f3 < 0))
{
rho1=rho3;
}
}
else if ((f1 < 0 && f3 > 0)||(f1 > 0 && f3 < 0))
{
rho2=rho3;
@ -416,7 +416,7 @@ inline scalar aungierRedlichKwong::rho
{
molarVolume=this->W()/rho3;
molarVolumePrevIteration=this->W()/rho3;
break;
break;
}
else
{

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

@ -1,5 +1,5 @@
CL: Real gas mixture classes using mixture models.
CL: other models are possible
CL: other models are possible
CL: Mixtures based on a pseudo critical point approach

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

@ -28,7 +28,7 @@ Description
Author
Christian Lucas
Institut für Thermodynamik
Technische Universität Braunschweig
Technische Universität Braunschweig
Germany
\*---------------------------------------------------------------------------*/

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

@ -28,8 +28,8 @@ Class
Description
Mixture Aungier Redlich Kwong equation of state.
Mixture based on a pseudo critical point approach.
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
@ -42,7 +42,7 @@ SourceFiles
Author
Christian Lucas
Institut für Thermodynamik
Technische Universität Braunschweig
Technische Universität Braunschweig
Germany
\*---------------------------------------------------------------------------*/
@ -73,27 +73,27 @@ 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,
const aungierRedlichKwong& sp,
scalar Tcrit,
scalar azentricFactor,
scalar Vcrit,

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

@ -25,7 +25,7 @@ License
Author
Christian Lucas
Institut für Thermodynamik
Technische Universität Braunschweig
Technische Universität Braunschweig
Germany
\*---------------------------------------------------------------------------*/
@ -51,14 +51,14 @@ inline mixtureAungierRedlichKwong::mixtureAungierRedlichKwong
scalar rhoMax
)
:
aungierRedlichKwong(sp),
aungierRedlichKwong(sp),
Vcrit_(Vcrit),
Zcrit_(Zcrit)
{
//CL: Saving critical data
Tcrit_=Tcrit;
pcrit_=Zcrit*this->RR*Tcrit/Vcrit;
rhocrit_=this->W()/Vcrit_;
pcrit_=Zcrit*this->RR*Tcrit/Vcrit;
rhocrit_=this->W()/Vcrit_;
azentricFactor_=azentricFactor;
//CL: calculating the aungier redlich kwong coefficience
@ -71,7 +71,7 @@ inline mixtureAungierRedlichKwong::mixtureAungierRedlichKwong
rhoMin_=rhoMin;
rhoMax_=rhoMax;
//CL:
//CL:
b2_=pow(b_,2);
b3_=pow(b_,3);
b4_=pow(b_,4);
@ -128,7 +128,7 @@ inline void mixtureAungierRedlichKwong::operator+=(const mixtureAungierRedlichKw
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_;
@ -142,7 +142,7 @@ inline void mixtureAungierRedlichKwong::operator+=(const mixtureAungierRedlichKw
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:
//CL:
b2_=pow(b_,2);
b3_=pow(b_,3);
b4_=pow(b_,4);
@ -195,11 +195,11 @@ inline mixtureAungierRedlichKwong operator*
return mixtureAungierRedlichKwong
(
s*static_cast<const aungierRedlichKwong&>(pr),
pr.Tcrit_,
pr.azentricFactor_,
pr.Vcrit_,
pr.Zcrit_,
pr.rhoMin_,
pr.Tcrit_,
pr.azentricFactor_,
pr.Vcrit_,
pr.Zcrit_,
pr.rhoMin_,
pr.rhoMax_
);
}

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

@ -28,7 +28,7 @@ Description
Author
Christian Lucas
Institut für Thermodynamik
Technische Universität Braunschweig
Technische Universität Braunschweig
Germany
\*---------------------------------------------------------------------------*/
@ -50,7 +50,7 @@ mixturePengRobinson::mixturePengRobinson(Istream& is)
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
@ -101,21 +101,21 @@ mixturePengRobinson::mixturePengRobinson(const dictionary& dict)
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
//CL: determine the size of the vector
for (i=3; i<=nCom_;i++)
{
{
sizeOfVector_=sizeOfVector_+(i-1);
}
//CL: setting the size
//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

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

@ -41,7 +41,7 @@ SourceFiles
Author
Christian Lucas
Institut für Thermodynamik
Technische Universität Braunschweig
Technische Universität Braunschweig
Germany
\*---------------------------------------------------------------------------*/
@ -77,34 +77,34 @@ 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;
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: 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;
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_;
@ -121,7 +121,7 @@ public:
label numOfComp,
DynamicList<scalar> weigths,
DynamicList<mixturePengRobinson*> mixtureComponents,
scalar a0,
scalar a0,
scalar b,
scalar Tcrit,
scalar n,
@ -161,29 +161,29 @@ public:
//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
//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 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
(
const scalar rho,
const scalar T
) const;
) const;
// Used for cv
inline scalar integral_d2pdT2_dv
@ -193,7 +193,7 @@ public:
) const;
//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;
@ -218,7 +218,7 @@ public:
inline scalar rho(const scalar p, const scalar T) const;
//- Return compressibility drho/dp at T=constant [s^2/m^2]
//- Return compressibility drho/dp at T=constant [s^2/m^2]
inline scalar psi(const scalar rho, const scalar T) const;
//- Return compression factor []

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

@ -25,7 +25,7 @@ License
Author
Christian Lucas
Institut für Thermodynamik
Technische Universität Braunschweig
Technische Universität Braunschweig
Germany
\*---------------------------------------------------------------------------*/
@ -40,7 +40,7 @@ namespace Foam
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
// Construct from components
// CL: needed for operator*
// CL: needed for operator*
inline mixturePengRobinson::mixturePengRobinson
(
const pengRobinson& pr,
@ -55,7 +55,7 @@ inline mixturePengRobinson::mixturePengRobinson
scalar rhoMin,
scalar rhoMax
)
:
:
pengRobinson(pr),
singleComponent(1),
weigths(weigths),
@ -64,7 +64,7 @@ inline mixturePengRobinson::mixturePengRobinson
{
a0_=a0;
b_=b;
Tcrit_=Tcrit;
Tcrit_=Tcrit;
n_=n;
rhostd_=rhostd;
rhoMin_=rhoMin;
@ -82,12 +82,12 @@ inline mixturePengRobinson::mixturePengRobinson
scalar rhoMin,
scalar rhoMax
)
:
:
pengRobinson(pr),
singleComponent(0),
weigths(weigths),
mixtureComponents(mixtureComponents),
numOfComp(numOfComp)
numOfComp(numOfComp)
{
TSave=0.0;
rhoMin_=rhoMin;
@ -120,23 +120,23 @@ inline autoPtr<mixturePengRobinson> mixturePengRobinson::New(Istream& is)
//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
{
// 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: getting a(T), dadT(T) and d2adT2(T) for the mixture
//CL: using van der waals mixing rule
for (j=0;j<numOfComp;j++)
{
@ -144,22 +144,22 @@ inline void mixturePengRobinson::updateModelCoefficients(const scalar T) const
//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)
// 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
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]);
}
d2aSave=d2aSave+d2adT2mix(T,i,j)*(1-mixtureComponents[0]->realMixtureCorrCoef_[k]);
}
}
else
@ -167,17 +167,17 @@ inline void mixturePengRobinson::updateModelCoefficients(const scalar T) const
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: 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);
@ -229,7 +229,7 @@ inline scalar mixturePengRobinson::d2adT2mix(const scalar T, const label i, cons
//returns the pressure for a given density and temperature
//returns the pressure for a given density and temperature
inline scalar mixturePengRobinson::p(const scalar rho, const scalar T) const
{
updateModelCoefficients(T);
@ -237,7 +237,7 @@ inline scalar mixturePengRobinson::p(const scalar rho, const scalar T) const
}
//Real deviative dp/dv at constant temperature
//Real deviative dp/dv at constant temperature
//(molar values)
inline scalar mixturePengRobinson::dpdv(const scalar rho, const scalar T) const
{
@ -246,7 +246,7 @@ inline scalar mixturePengRobinson::dpdv(const scalar rho, const scalar T) const
}
//Real deviative dp/dT at constant molar volume
//Real deviative dp/dT at constant molar volume
//(molar values)
inline scalar mixturePengRobinson::dpdT(const scalar rho, const scalar T) const
{
@ -255,7 +255,7 @@ inline scalar mixturePengRobinson::dpdT(const scalar rho, const scalar T) const
}
//Real deviative dv/dT at constant pressure
//Real deviative dv/dT at constant pressure
//(molar values)
inline scalar mixturePengRobinson::dvdT(const scalar rho, const scalar T) const
{
@ -264,7 +264,7 @@ inline scalar mixturePengRobinson::dvdT(const scalar rho, const scalar T) const
}
//Real deviative dv/dp at constant temperature
//Real deviative dv/dp at constant temperature
//(molar values)
inline scalar mixturePengRobinson::dvdp(const scalar rho, const scalar T) const
{
@ -280,7 +280,7 @@ inline scalar mixturePengRobinson::integral_p_dv
const scalar rho,
const scalar T
) const
{
{
updateModelCoefficients(T);
return pengRobinson::integral_p_dv(rho,T);
}
@ -304,7 +304,7 @@ inline scalar mixturePengRobinson::d2pdT2(const scalar rho, const scalar T) cons
{
updateModelCoefficients(T);
return pengRobinson::d2pdT2(rho,T);
}
}
//(molar values)
@ -331,7 +331,7 @@ inline scalar mixturePengRobinson::d2pdvdT(const scalar rho, const scalar T) con
}
// the result of this intergal is needed for the nasa based cp polynomial
// the result of this intergal is needed for the nasa based cp polynomial
//(molar values)
inline scalar mixturePengRobinson::integral_d2pdT2_dv
(
@ -370,7 +370,7 @@ inline scalar mixturePengRobinson::isothermalCompressiblity
}
//- Return density [kg/m^3]on
//- Return density [kg/m^3]on
inline scalar mixturePengRobinson::rho
(
const scalar p,
@ -383,7 +383,7 @@ inline scalar mixturePengRobinson::rho
}
//- Return density [kg/m^3]on
//- Return density [kg/m^3]on
inline scalar mixturePengRobinson::rho(const scalar p, const scalar T) const
{
updateModelCoefficients(T);
@ -422,9 +422,9 @@ inline void mixturePengRobinson::operator+=(const mixturePengRobinson& pr)
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
//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];
@ -465,7 +465,7 @@ inline mixturePengRobinson operator+
mixtureComponents.setSize(2*pr1.numOfComp);
}
//CL: Getting the new weigths and mixtureComponents lists,
//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];
@ -484,7 +484,7 @@ inline mixturePengRobinson operator*
//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();
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_);

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

@ -28,7 +28,7 @@ Description
Author
Christian Lucas
Institut für Thermodynamik
Technische Universität Braunschweig
Technische Universität Braunschweig
Germany
\*---------------------------------------------------------------------------*/
@ -49,7 +49,7 @@ mixtureRedlichKwong::mixtureRedlichKwong(Istream& 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
@ -98,22 +98,22 @@ mixtureRedlichKwong::mixtureRedlichKwong(const dictionary& dict)
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
//CL: determine the size of the vector
for (i=3; i<=nCom_;i++)
{
{
sizeOfVector_=sizeOfVector_+(i-1);
}
//CL: setting the size
//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

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

@ -41,7 +41,7 @@ SourceFiles
Author
Christian Lucas
Institut für Thermodynamik
Technische Universität Braunschweig
Technische Universität Braunschweig
Germany
\*---------------------------------------------------------------------------*/
@ -72,23 +72,23 @@ class mixtureRedlichKwong
:
public redlichKwong
{
// Private data
// 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;
//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;
@ -118,7 +118,7 @@ public:
label numOfComp,
DynamicList<scalar> weigths,
DynamicList<mixtureRedlichKwong*> mixtureComponents,
scalar a,
scalar a,
scalar b,
scalar rhostd,
scalar rhoMin,

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

@ -25,7 +25,7 @@ License
Author
Christian Lucas
Institut für Thermodynamik
Technische Universität Braunschweig
Technische Universität Braunschweig
Germany
\*---------------------------------------------------------------------------*/
@ -40,7 +40,7 @@ namespace Foam
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
// Construct from components
// CL: needed for operator*
// CL: needed for operator*
inline mixtureRedlichKwong::mixtureRedlichKwong
(
const redlichKwong& rK,
@ -53,7 +53,7 @@ inline mixtureRedlichKwong::mixtureRedlichKwong
scalar rhoMin,
scalar rhoMax
)
:
:
redlichKwong(rK),
weigths(weigths),
mixtureComponents(mixtureComponents),
@ -82,7 +82,7 @@ inline mixtureRedlichKwong::mixtureRedlichKwong
scalar rhoMin,
scalar rhoMax
)
:
:
redlichKwong(rK),
weigths(weigths),
mixtureComponents(mixtureComponents),
@ -126,15 +126,15 @@ inline autoPtr<mixtureRedlichKwong> mixtureRedlichKwong::New(Istream& is)
//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: getting a for the mixture
//CL: using van der waals mixing rule
for (j=0;j<numOfComp;j++)
{
@ -157,14 +157,14 @@ inline void mixtureRedlichKwong::updateModelCoefficients() const
//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
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);
@ -182,8 +182,8 @@ inline void mixtureRedlichKwong::operator+=(const mixtureRedlichKwong& rK)
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
//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];
@ -223,7 +223,7 @@ inline mixtureRedlichKwong operator+
mixtureComponents.setSize(2*rK1.numOfComp);
}
//CL: Getting the new weigths and mixtureComponents lists,
//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];
@ -241,9 +241,9 @@ inline mixtureRedlichKwong operator*
{
//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();
weigths[rK.numOfComp-1]=s*rK.nMoles();
return mixtureRedlichKwong(s*static_cast<const redlichKwong&>(rK), rK.numOfComp, weigths,
return mixtureRedlichKwong(s*static_cast<const redlichKwong&>(rK), rK.numOfComp, weigths,
rK.mixtureComponents,rK.a_,rK.b_,rK.rhostd_, rK.rhoMin_, rK.rhoMax_);
}

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

@ -28,7 +28,7 @@ Description
Author
Christian Lucas
Institut für Thermodynamik
Technische Universität Braunschweig
Technische Universität Braunschweig
Germany
\*---------------------------------------------------------------------------*/
@ -50,7 +50,7 @@ mixtureSoaveRedlichKwong::mixtureSoaveRedlichKwong(Istream& is)
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
@ -72,7 +72,7 @@ mixtureSoaveRedlichKwong::mixtureSoaveRedlichKwong(const dictionary& 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
@ -101,21 +101,21 @@ mixtureSoaveRedlichKwong::mixtureSoaveRedlichKwong(const dictionary& dict)
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
//CL: determine the size of the vector
for (i=3; i<=nCom_;i++)
{
{
sizeOfVector_=sizeOfVector_+(i-1);
}
//CL: setting the size
//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

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

@ -41,7 +41,7 @@ SourceFiles
Author
Christian Lucas
Institut für Thermodynamik
Technische Universität Braunschweig
Technische Universität Braunschweig
Germany
\*---------------------------------------------------------------------------*/
@ -77,33 +77,33 @@ 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;
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: 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_;
@ -120,7 +120,7 @@ public:
label numOfComp,
DynamicList<scalar> weigths,
DynamicList<mixtureSoaveRedlichKwong*> mixtureComponents,
scalar a0,
scalar a0,
scalar b,
scalar Tcrit,
scalar n,
@ -160,29 +160,29 @@ public:
//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
//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 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
(
const scalar rho,
const scalar T
) const;
) const;
// Used for cv
inline scalar integral_d2pdT2_dv
@ -192,7 +192,7 @@ public:
) const;
//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;
@ -217,7 +217,7 @@ public:
inline scalar rho(const scalar p, const scalar T) const;
//- Return compressibility drho/dp at T=constant [s^2/m^2]
//- Return compressibility drho/dp at T=constant [s^2/m^2]
inline scalar psi(const scalar rho, const scalar T) const;
//- Return compression factor []

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

@ -25,7 +25,7 @@ License
Author
Christian Lucas
Institut für Thermodynamik
Technische Universität Braunschweig
Technische Universität Braunschweig
Germany
\*---------------------------------------------------------------------------*/
@ -40,7 +40,7 @@ namespace Foam
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
// Construct from components
// CL: needed for operator*
// CL: needed for operator*
inline mixtureSoaveRedlichKwong::mixtureSoaveRedlichKwong
(
const soaveRedlichKwong& srk,
@ -55,7 +55,7 @@ inline mixtureSoaveRedlichKwong::mixtureSoaveRedlichKwong
scalar rhoMin,
scalar rhoMax
)
:
:
soaveRedlichKwong(srk),
singleComponent(1),
weigths(weigths),
@ -64,7 +64,7 @@ inline mixtureSoaveRedlichKwong::mixtureSoaveRedlichKwong
{
a0_=a0;
b_=b;
Tcrit_=Tcrit;
Tcrit_=Tcrit;
n_=n;
rhostd_=rhostd;
rhoMin_=rhoMin;
@ -83,7 +83,7 @@ inline mixtureSoaveRedlichKwong::mixtureSoaveRedlichKwong
scalar rhoMin,
scalar rhoMax
)
:
:
soaveRedlichKwong(srk),
singleComponent(0),
weigths(weigths),
@ -123,23 +123,23 @@ inline autoPtr<mixtureSoaveRedlichKwong> mixtureSoaveRedlichKwong::New(Istream&
//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
{
// 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: getting a(T), dadT(T) and d2adT2(T) for the mixture
//CL: using van der waals mixing rule
for (j=0;j<numOfComp;j++)
{
@ -147,22 +147,22 @@ inline void mixtureSoaveRedlichKwong::updateModelCoefficients(const scalar T) c
//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)
// 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
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]);
}
d2aSave=d2aSave+d2adT2mix(T,i,j)*(1-mixtureComponents[0]->realMixtureCorrCoef_[k]);
}
}
else
@ -170,17 +170,17 @@ inline void mixtureSoaveRedlichKwong::updateModelCoefficients(const scalar T) c
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: 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);
@ -228,7 +228,7 @@ inline scalar mixtureSoaveRedlichKwong::d2adT2mix(const scalar T, const label i,
);
}
//returns the pressure for a given density and temperature
//returns the pressure for a given density and temperature
inline scalar mixtureSoaveRedlichKwong::p(const scalar rho, const scalar T) const
{
updateModelCoefficients(T);
@ -236,7 +236,7 @@ inline scalar mixtureSoaveRedlichKwong::p(const scalar rho, const scalar T) cons
}
//Real deviative dp/dv at constant temperature
//Real deviative dp/dv at constant temperature
//(molar values)
inline scalar mixtureSoaveRedlichKwong::dpdv(const scalar rho, const scalar T) const
{
@ -245,7 +245,7 @@ inline scalar mixtureSoaveRedlichKwong::dpdv(const scalar rho, const scalar T) c
}
//Real deviative dp/dT at constant molar volume
//Real deviative dp/dT at constant molar volume
//(molar values)
inline scalar mixtureSoaveRedlichKwong::dpdT(const scalar rho, const scalar T) const
{
@ -254,7 +254,7 @@ inline scalar mixtureSoaveRedlichKwong::dpdT(const scalar rho, const scalar T) c
}
//Real deviative dv/dT at constant pressure
//Real deviative dv/dT at constant pressure
//(molar values)
inline scalar mixtureSoaveRedlichKwong::dvdT(const scalar rho, const scalar T) const
{
@ -263,7 +263,7 @@ inline scalar mixtureSoaveRedlichKwong::dvdT(const scalar rho, const scalar T) c
}
//Real deviative dv/dp at constant temperature
//Real deviative dv/dp at constant temperature
//(molar values)
inline scalar mixtureSoaveRedlichKwong::dvdp(const scalar rho, const scalar T) const
{
@ -279,7 +279,7 @@ inline scalar mixtureSoaveRedlichKwong::integral_p_dv
const scalar rho,
const scalar T
) const
{
{
updateModelCoefficients(T);
return soaveRedlichKwong::integral_p_dv(rho,T);
}
@ -303,7 +303,7 @@ inline scalar mixtureSoaveRedlichKwong::d2pdT2(const scalar rho, const scalar T)
{
updateModelCoefficients(T);
return soaveRedlichKwong::d2pdT2(rho,T);
}
}
//(molar values)
@ -330,7 +330,7 @@ inline scalar mixtureSoaveRedlichKwong::d2pdvdT(const scalar rho, const scalar T
}
// the result of this intergal is needed for the nasa based cp polynomial
// the result of this intergal is needed for the nasa based cp polynomial
//(molar values)
inline scalar mixtureSoaveRedlichKwong::integral_d2pdT2_dv
(
@ -369,7 +369,7 @@ inline scalar mixtureSoaveRedlichKwong::isothermalCompressiblity
}
//- Return density [kg/m^3]on
//- Return density [kg/m^3]on
inline scalar mixtureSoaveRedlichKwong::rho
(
const scalar p,
@ -383,7 +383,7 @@ inline scalar mixtureSoaveRedlichKwong::rho
}
//- Return density [kg/m^3]on
//- Return density [kg/m^3]on
inline scalar mixtureSoaveRedlichKwong::rho(const scalar p, const scalar T) const
{
updateModelCoefficients(T);
@ -423,8 +423,8 @@ inline void mixtureSoaveRedlichKwong::operator+=(const mixtureSoaveRedlichKwong&
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
//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];
@ -466,7 +466,7 @@ inline mixtureSoaveRedlichKwong operator+
mixtureComponents.setSize(2*srk1.numOfComp);
}
//CL: Getting the new weigths and mixtureComponents lists,
//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];
@ -485,7 +485,7 @@ inline mixtureSoaveRedlichKwong operator*
//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();
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_);

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

@ -28,7 +28,7 @@ Description
Author
Christian Lucas
Institut für Thermodynamik
Technische Universität Braunschweig
Technische Universität Braunschweig
Germany
\*---------------------------------------------------------------------------*/
@ -81,14 +81,14 @@ pengRobinson::pengRobinson(const dictionary& dict)
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_),
b_(0.077796*this->RR()*Tcrit_/pcrit_),
n_(0.37464+1.54226*azentricFactor_-0.26992*pow(azentricFactor_,2)),
TSave(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())))
rhostd_(this->rho(Pstd,Tstd,Pstd/(Tstd*this->R())))
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
@ -111,7 +111,7 @@ void Foam::pengRobinson::write(Ostream& os) const
Ostream& operator<<(Ostream& os, const pengRobinson& pr)
{
os << static_cast<const specie&>(pr)<< token::SPACE
os << static_cast<const specie&>(pr)<< token::SPACE
<< pr.pcrit_ << tab<< pr.Tcrit_<< tab << pr.azentricFactor_;
os.check("Ostream& operator<<(Ostream& os, const pengRobinson& st)");

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

@ -42,7 +42,7 @@ SourceFiles
Author
Christian Lucas
Institut für Thermodynamik
Technische Universität Braunschweig
Technische Universität Braunschweig
Germany
\*---------------------------------------------------------------------------*/
@ -72,11 +72,11 @@ class pengRobinson
protected:
// Protected data
scalar pcrit_;
scalar Tcrit_;
scalar Tcrit_;
scalar azentricFactor_;
//-Peng Robinson factors
scalar n_;
scalar n_;
scalar a0_;
mutable scalar b_;
@ -104,10 +104,10 @@ protected:
//Protected functions
//CL: function updates the coefficients (aSave, daSave, d2aSave)
inline void updateModelCoefficients(const scalar T) const;
inline void updateModelCoefficients(const scalar T) const;
public:
// Constructors
@ -163,33 +163,33 @@ public:
//CL: Equation of state
inline scalar p(const scalar rho, const scalar T) const;
//CL: 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;
inline scalar dvdT(const scalar rho,const scalar T) const;
inline scalar dvdp(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;
) const;
inline scalar isothermalCompressiblity
(
const scalar rho,
const scalar T
) const;
) const;
//CL: Used for cv
inline scalar integral_d2pdT2_dv
(
(
const scalar rho,
const scalar T
) const ;
) const ;
//CL: second order derivatives, not Used At The Moment
inline scalar d2pdv2(const scalar rho,const scalar T) const;
@ -202,11 +202,11 @@ public:
//CL: Used for internal Energy
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;
//- Return density [kg/m^3]
//- Return density [kg/m^3]
// rho0 is the starting point of the newton solver used to calculate rho
inline scalar rho
(

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

@ -25,7 +25,7 @@ License
Author
Christian Lucas
Institut für Thermodynamik
Technische Universität Braunschweig
Technische Universität Braunschweig
Germany
\*---------------------------------------------------------------------------*/
@ -42,11 +42,11 @@ namespace Foam
// Construct from components
inline pengRobinson::pengRobinson
(
const specie& sp
const specie& sp
)
:
specie(sp),
TSave(0)
TSave(0)
{}
@ -107,15 +107,15 @@ inline scalar pengRobinson::a(const scalar T)const
{
//CL: check if a has already been calculated for this temperature
if(TSave==T)
{
return aSave;
{
return aSave;
}
//CL: If not, recalculate a(T), dadT(T) and d2adT2(T)
else
{
updateModelCoefficients(T);
return aSave;
}
}
}
@ -125,14 +125,14 @@ inline scalar pengRobinson::dadT(const scalar T)const
// check if a has already been calculated for this temperature
if(TSave==T)
{
return daSave;
return daSave;
}
//CL: If not, recalculate a(T), dadT(T) and d2adT2(T)
else
{
updateModelCoefficients(T);
return daSave;
}
}
}
@ -142,14 +142,14 @@ inline scalar pengRobinson::d2adT2(const scalar T)const
// check if a has already been calculated for this temperature
if(TSave==T)
{
return d2aSave;
return d2aSave;
}
//CL: If not, recalculate a(T), dadT(T) and d2adT2(T)
else
{
updateModelCoefficients(T);
return d2aSave;
}
}
}
@ -191,16 +191,16 @@ inline scalar pengRobinson::b5()const
return b5_;
}
//returns the pressure for a given density and temperature
//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());
}
//Real deviative dp/dv at constant temperature
//Real deviative dp/dv at constant temperature
//(molar values)
inline scalar pengRobinson::dpdv(const scalar rho, const scalar T) const
{
@ -223,7 +223,7 @@ inline scalar pengRobinson::dpdv(const scalar rho, const scalar T) const
}
//Real deviative dp/dT at constant molar volume
//Real deviative dp/dT at constant molar volume
//(molar values)
inline scalar pengRobinson::dpdT(const scalar rho, const scalar T) const
{
@ -233,7 +233,7 @@ inline scalar pengRobinson::dpdT(const scalar rho, const scalar T) const
}
//Real deviative dv/dT at constant pressure
//Real deviative dv/dT at constant pressure
//by using implicit differentiation
//(molar values)
inline scalar pengRobinson::dvdT(const scalar rho,const scalar T) const
@ -243,7 +243,7 @@ inline scalar pengRobinson::dvdT(const scalar rho,const scalar T) const
//(molar values)
inline scalar pengRobinson::dvdp(const scalar rho,const scalar T) const
inline scalar pengRobinson::dvdp(const scalar rho,const scalar T) const
{
return 1/this->dpdv(rho,T);
}
@ -306,10 +306,10 @@ inline scalar pengRobinson::d2pdv2(const scalar rho,const scalar T) const
//(molar values)
//using second order implicit differentiation
//using second order implicit differentiation
inline scalar pengRobinson::d2vdT2(const scalar rho, const scalar T) const
{
return
return
-(
pow(this->dpdT(rho,T),2)*this->d2pdv2(rho,T)
+ pow(this->dpdv(rho,T),2)*this->d2pdT2(rho,T)
@ -325,7 +325,7 @@ 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;
return(
2*dadT(T)*
(
@ -341,9 +341,9 @@ inline scalar pengRobinson::d2pdvdT(const scalar rho, const scalar T) const
}
// the result of this intergal is needed for the nasa based cp polynomial
// 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;
@ -364,7 +364,7 @@ inline scalar pengRobinson::isobarExpCoef(const scalar rho,const scalar T) cons
inline scalar pengRobinson::isothermalCompressiblity(const scalar rho,const scalar T) const
{
return this->isobarExpCoef(rho, T)/this->dpdT(rho, T);
//CL: also possible
//CL: also possible
//CL: return -this->dvdp(rho,T)*rho/this->W();
}
@ -401,7 +401,7 @@ inline scalar pengRobinson::rho(
)/pow(2,i);
i++;
if (i>8)
if (i>8)
{
//CL: using bisection methode as backup,
//CL: solution must be between rho=0.001 to rho=1500;
@ -416,7 +416,7 @@ inline scalar pengRobinson::rho(
if ((f2 < 0 && f3 > 0) || (f2 > 0 && f3 < 0))
{
rho1=rho3;
}
}
else if ((f1 < 0 && f3 > 0)||(f1 > 0 && f3 < 0))
{
rho2=rho3;
@ -430,7 +430,7 @@ inline scalar pengRobinson::rho(
{
molarVolume=this->W()/rho3;
molarVolumePrevIteration=this->W()/rho3;
break;
break;
}
else
{
@ -460,15 +460,15 @@ inline scalar pengRobinson::rho(
}
//- Return density [kg/m^3]on
//- Return density [kg/m^3]on
inline scalar pengRobinson::rho(const scalar p,const scalar T) const
{
// using perfect gas equation as starting point
return rho(p,T,p/(this->R()*T));
return rho(p,T,p/(this->R()*T));
}
//- Return compressibility drho/dp at T=constant [s^2/m^2]
//- 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();

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

@ -28,7 +28,7 @@ Description
Author
Christian Lucas
Institut für Thermodynamik
Technische Universität Braunschweig
Technische Universität Braunschweig
Germany
\*---------------------------------------------------------------------------*/
@ -58,7 +58,7 @@ redlichKwong::redlichKwong(Istream& is)
rhoMin_(1e-3),
// 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
@ -81,7 +81,7 @@ redlichKwong::redlichKwong(const dictionary& dict)
b3_(pow(b_,3)),
b5_(pow(b_,5)),
// Starting GUESS for the density by ideal gas law
rhostd_(this->rho(Pstd, Tstd, Pstd/(Tstd*this->R())))
rhostd_(this->rho(Pstd, Tstd, Pstd/(Tstd*this->R())))
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
@ -104,7 +104,7 @@ void Foam::redlichKwong::write(Ostream& os) const
Ostream& operator<<(Ostream& os, const redlichKwong& rk)
{
os << static_cast<const specie&>(rk)<< token::SPACE
os << static_cast<const specie&>(rk)<< token::SPACE
<< rk.pcrit_ << tab<< rk.Tcrit_;
os.check("Ostream& operator<<(Ostream& os, const redlichKwong& st)");

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

@ -35,7 +35,7 @@ SourceFiles
Author
Christian Lucas
Institut für Thermodynamik
Technische Universität Braunschweig
Technische Universität Braunschweig
Germany
\*---------------------------------------------------------------------------*/
@ -65,7 +65,7 @@ protected:
// Protected data
scalar pcrit_;
scalar Tcrit_;
scalar Tcrit_;
//-Redlich Kwong factors
mutable scalar a_;
@ -84,7 +84,7 @@ protected:
mutable scalar rhostd_;
public:
// Constructors
//- Construct from components
@ -127,22 +127,22 @@ public:
inline scalar b5()const;
//first order derivatives
//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 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
(
const scalar rho,
const scalar T
) const;
) const;
// Used for cv
inline scalar integral_d2pdT2_dv
@ -152,7 +152,7 @@ public:
) const;
//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;
@ -177,7 +177,7 @@ public:
inline scalar rho(const scalar p, const scalar T) const;
//- Return compressibility drho/dp at T=constant [s^2/m^2]
//- Return compressibility drho/dp at T=constant [s^2/m^2]
inline scalar psi(const scalar rho, const scalar T) const;
//- Return compression factor []
@ -198,7 +198,7 @@ public:
// Member operators
inline void operator+=(const redlichKwong&);
/*
/*
inline void operator-=(const redlichKwong&);
inline void operator*=(const scalar);
*/

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

@ -25,7 +25,7 @@ License
Author
Christian Lucas
Institut für Thermodynamik
Technische Universität Braunschweig
Technische Universität Braunschweig
Germany
\*---------------------------------------------------------------------------*/
@ -114,15 +114,15 @@ inline scalar redlichKwong::b5()const
//returns the pressure for a given density and temperature
//returns the pressure for a given density and temperature
inline scalar redlichKwong::p(const scalar rho, const scalar T) const
{
scalar Vm = this->W()/rho;
scalar Vm = this->W()/rho;
return this->RR()*T/(Vm - b_) - a_/(sqrt(T)*Vm*(Vm + b_));
}
//Real deviative dp/dv at constant temperature
//Real deviative dp/dv at constant temperature
//(molar values)
inline scalar redlichKwong::dpdv(const scalar rho, const scalar T) const
{
@ -135,7 +135,7 @@ inline scalar redlichKwong::dpdv(const scalar rho, const scalar T) const
}
//Real deviative dp/dT at constant molar volume
//Real deviative dp/dT at constant molar volume
//(molar values)
inline scalar redlichKwong::dpdT(const scalar rho, const scalar T) const
{
@ -144,7 +144,7 @@ inline scalar redlichKwong::dpdT(const scalar rho, const scalar T) const
}
//Real deviative dv/dT at constant pressure
//Real deviative dv/dT at constant pressure
//using implicit differentiation
//(molar values)
inline scalar redlichKwong::dvdT(const scalar rho, const scalar T) const
@ -153,9 +153,9 @@ inline scalar redlichKwong::dvdT(const scalar rho, const scalar T) const
}
//Real deviative dv/dp at constant temperature
//Real deviative dv/dp at constant temperature
//(molar values)
inline scalar redlichKwong::dvdp(const scalar rho, const scalar T) const
inline scalar redlichKwong::dvdp(const scalar rho, const scalar T) const
{
return 1/this->dpdv(rho,T);
}
@ -224,10 +224,10 @@ 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
{
return
return
-(
pow(this->dpdT(rho,T),2)*this->d2pdv2(rho,T)
+ pow(this->dpdv(rho,T),2)*this->d2pdT2(rho,T)
@ -242,7 +242,7 @@ inline scalar redlichKwong::d2pdvdT(const scalar rho, const scalar T) const
{
scalar Vm = this->W()/rho;
scalar Vm2 = Vm*Vm;
return
-(0.5*(
a_*(b3() - 3*b_*Vm2 + 2*pow(Vm,3))
@ -252,13 +252,13 @@ inline scalar redlichKwong::d2pdvdT(const scalar rho, const scalar T) const
}
// the result of this intergal is needed for the nasa based cp polynomial
// the result of this intergal is needed for the nasa based cp polynomial
//(molar values)
inline scalar redlichKwong::integral_d2pdT2_dv
(
const scalar rho,
const scalar T
) const
) const
{
scalar Vm = this->W()/rho;
return 0.75*a_*log(b_ + Vm)/(pow(T,2.5)*b_)
@ -299,7 +299,7 @@ inline scalar redlichKwong::rho
const scalar rho0
) const
{
scalar molarVolumePrevIteration;
scalar molarVolume;
label iter=0;
@ -313,7 +313,7 @@ inline scalar redlichKwong::rho
do
{
molarVolumePrevIteration= molarVolume;
label i=0;
do
{
@ -325,7 +325,7 @@ inline scalar redlichKwong::rho
)/pow(2,i);
i++;
if (i>8)
if (i>8)
{
//CL: using bisection methode as backup,
//CL: solution must be between rhoMin_ to rhoMax
@ -339,7 +339,7 @@ inline scalar redlichKwong::rho
if ((f2 < 0 && f3 > 0) || (f2 > 0 && f3 < 0))
{
rho1=rho3;
}
}
else if ((f1 < 0 && f3 > 0)||(f1 > 0 && f3 < 0))
{
rho2=rho3;
@ -353,7 +353,7 @@ inline scalar redlichKwong::rho
{
molarVolume=this->W()/rho3;
molarVolumePrevIteration=this->W()/rho3;
break;
break;
}
else
{
@ -383,7 +383,7 @@ inline scalar redlichKwong::rho
}
//- Return density [kg/m^3]on
//- Return density [kg/m^3]on
inline scalar redlichKwong::rho(const scalar p, const scalar T) const
{
// using perfect gas equation as starting point

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

@ -28,7 +28,7 @@ Description
Author
Christian Lucas
Institut für Thermodynamik
Technische Universität Braunschweig
Technische Universität Braunschweig
Germany
\*---------------------------------------------------------------------------*/
@ -82,7 +82,7 @@ soaveRedlichKwong::soaveRedlichKwong(const dictionary& dict)
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),
TSave(0.0),
b2_(pow(b_,2)),
b3_(pow(b_,3)),
b5_(pow(b_,5)),
@ -110,7 +110,7 @@ void Foam::soaveRedlichKwong::write(Ostream& os) const
Ostream& operator<<(Ostream& os, const soaveRedlichKwong& srk)
{
os << static_cast<const specie&>(srk)<< token::SPACE
os << static_cast<const specie&>(srk)<< token::SPACE
<< srk.pcrit_ << tab<< srk.Tcrit_<<tab<<srk.azentricFactor_;
os.check("Ostream& operator<<(Ostream& os, const soaveRedlichKwong& st)");

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

@ -40,7 +40,7 @@ SourceFiles
Author
Christian Lucas
Institut für Thermodynamik
Technische Universität Braunschweig
Technische Universität Braunschweig
Germany
\*---------------------------------------------------------------------------*/
@ -70,11 +70,11 @@ class soaveRedlichKwong
protected:
// Protected data
scalar pcrit_;
scalar Tcrit_;
scalar Tcrit_;
scalar azentricFactor_;
//-Soave Redlich Kwong
scalar n_;
scalar n_;
scalar a0_;
mutable scalar b_;
@ -102,10 +102,10 @@ protected:
//Protected functions
//CL: function updates the coefficients (aSave, daSave, d2aSave)
inline void updateModelCoefficients(const scalar T) const;
inline void updateModelCoefficients(const scalar T) const;
public:
// Constructors
@ -129,7 +129,7 @@ public:
// Selector from Istream
inline static autoPtr<soaveRedlichKwong> New(Istream& is);
// Member functions
inline scalar rhostd()const;
@ -159,46 +159,46 @@ public:
//CL: Equation of state
inline scalar p(const scalar rho, const scalar T) const;
//CL: 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;
inline scalar dvdT(const scalar rho,const scalar T) const;
inline scalar dvdp(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;
) const;
inline scalar isothermalCompressiblity
(
const scalar rho,
const scalar T
) const;
) 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;
inline scalar d2pdv2(const scalar rho,const scalar T) const;
inline scalar d2pdT2(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 d2pdvdT(const scalar rho,const scalar T) const;
inline scalar d2vdT2(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;
//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]
//- Return density [kg/m^3]
// rho0 is the starting point of the newton solver used to calculate rho
inline scalar rho
(
@ -209,7 +209,7 @@ public:
inline scalar rho(const scalar p,const scalar T) const;
//- Return compressibility drho/dp at T=constant [s^2/m^2]
//- Return compressibility drho/dp at T=constant [s^2/m^2]
inline scalar psi(const scalar rho, const scalar T) const;
//- Return compression factor []
@ -232,7 +232,7 @@ public:
inline void operator+=(const soaveRedlichKwong&);
/*
/*
inline void operator-=(const soaveRedlichKwong&);
inline void operator*=(const scalar);

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

@ -25,7 +25,7 @@ License
Author
Christian Lucas
Institut für Thermodynamik
Technische Universität Braunschweig
Technische Universität Braunschweig
Germany
@ -43,11 +43,11 @@ namespace Foam
// Construct from components
inline soaveRedlichKwong::soaveRedlichKwong
(
const specie& sp
const specie& sp
)
:
specie(sp),
TSave(0)
TSave(0)
{}
@ -107,15 +107,15 @@ inline scalar soaveRedlichKwong::a(const scalar T)const
{
//CL: check if a has already been calculated for this temperature
if(TSave==T)
{
return aSave;
{
return aSave;
}
//CL: If not, recalculate a(T), dadT(T) and d2adT2(T)
else
{
updateModelCoefficients(T);
return aSave;
}
}
}
@ -125,14 +125,14 @@ inline scalar soaveRedlichKwong::dadT(const scalar T)const
// check if a has already been calculated for this temperature
if(TSave==T)
{
return daSave;
return daSave;
}
//CL: If not, recalculate a(T), dadT(T) and d2adT2(T)
else
{
updateModelCoefficients(T);
return daSave;
}
}
}
@ -142,14 +142,14 @@ inline scalar soaveRedlichKwong::d2adT2(const scalar T)const
// check if a has already been calculated for this temperature
if(TSave==T)
{
return d2aSave;
return d2aSave;
}
//CL: If not, recalculate a(T), dadT(T) and d2adT2(T)
else
{
updateModelCoefficients(T);
return d2aSave;
}
}
}
@ -187,12 +187,12 @@ inline scalar soaveRedlichKwong::b5()const
}
//returns the pressure for a given density and temperature
//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
return
(
this->RR()*T/(Vm-b_)
-a(T)/(Vm*(Vm+b_))
@ -200,29 +200,29 @@ inline scalar soaveRedlichKwong::p(const scalar rho,const scalar T) const
}
//Real deviative dp/dv at constant temperature
//Real deviative dp/dv at constant temperature
//(molar values)
inline scalar soaveRedlichKwong::dpdv(const scalar rho, const scalar T) const
{
scalar Vm = this->W()/rho;
scalar Vm2 = Vm*Vm;
return
return
(
a(T)*(b3()-3*b_*Vm2+2*pow(Vm,3))
-this->RR()*T*Vm2*(b2()+2*b_*Vm+Vm2)
)
/(Vm2*pow(b_+Vm,2)*pow(b_-Vm,2));
/(Vm2*pow(b_+Vm,2)*pow(b_-Vm,2));
}
//Real deviative dp/dT at constant molar volume
//Real deviative dp/dT at constant molar volume
//(molar values)
inline scalar soaveRedlichKwong::dpdT(const scalar rho, const scalar T) const
{
scalar Vm = this->W()/rho;
return
return
(
this->RR()/(Vm-b_)
-dadT(T)/(Vm*(Vm+b_))
@ -239,20 +239,20 @@ inline scalar soaveRedlichKwong::dvdT(const scalar rho,const scalar T) const
}
//Real deviative dv/dp at constant temperature
//Real deviative dv/dp at constant temperature
//(molar values)
inline scalar soaveRedlichKwong::dvdp(const scalar rho,const scalar T) const
inline scalar soaveRedlichKwong::dvdp(const scalar rho,const scalar T) const
{
return 1/this->dpdv(rho,T);
}
//needed to calculate the internal energy
//needed to calculate the internal energy
//(molar values)
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_;
}
@ -276,14 +276,14 @@ inline scalar soaveRedlichKwong::d2pdT2(const scalar rho,const scalar T) const
}
//(molar values)
//(molar values)
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;
return
return
2*
(
a(T)*
@ -295,15 +295,15 @@ inline scalar soaveRedlichKwong::d2pdv2(const scalar rho,const scalar T) const
b3()+3*b2()*Vm+3*b_*Vm2+Vm3
)
)
/(Vm3*pow(b_+Vm,3)*pow(Vm-b_,3));
/(Vm3*pow(b_+Vm,3)*pow(Vm-b_,3));
}
//(molar values)
// using second Order implicit differentiation
// using second Order implicit differentiation
inline scalar soaveRedlichKwong::d2vdT2(const scalar rho, const scalar T) const
{
return
return
-(
pow(this->dpdT(rho,T),2)*this->d2pdv2(rho,T)
+ pow(this->dpdv(rho,T),2)*this->d2pdT2(rho,T)
@ -319,19 +319,19 @@ inline scalar soaveRedlichKwong::d2pdvdT(const scalar rho, const scalar T) const
scalar Vm = this->W()/rho;
scalar Vm2 = Vm*Vm;
return
return
(
dadT(T)*(b3()-3*b_*Vm2+2*pow(Vm,3))
-this->RR()*Vm2*(b2()+2*b_*Vm+Vm2)
)
/(Vm2*pow(b_+Vm,2)*pow(b_-Vm,2));
/(Vm2*pow(b_+Vm,2)*pow(b_-Vm,2));
}
// the result of this intergal is needed for the nasa based cp polynomial
// 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_;
@ -361,7 +361,7 @@ inline scalar soaveRedlichKwong::rho(
const scalar rho0
) const
{
scalar molarVolumePrevIteration;
scalar molarVolume;
label iter=0;
@ -375,7 +375,7 @@ inline scalar soaveRedlichKwong::rho(
do
{
molarVolumePrevIteration= molarVolume;
label i=0;
do
{
@ -386,7 +386,7 @@ inline scalar soaveRedlichKwong::rho(
)/pow(2,i);
i++;
if (i>8)
if (i>8)
{
//CL: using bisection methode as backup,
//CL: solution must be between rho=0.001 to rho=1500;
@ -401,11 +401,11 @@ inline scalar soaveRedlichKwong::rho(
if ((f2 < 0 && f3 > 0) || (f2 > 0 && f3 < 0))
{
rho1=rho3;
}
}
else if ((f1 < 0 && f3 > 0)||(f1 > 0 && f3 < 0))
{
rho2=rho3;
}
}
else
{
rho2=(rho2 + rho3)/2;
@ -415,7 +415,7 @@ inline scalar soaveRedlichKwong::rho(
{
molarVolume=this->W()/rho3;
molarVolumePrevIteration=this->W()/rho3;
break;
break;
}
else
{
@ -444,7 +444,7 @@ inline scalar soaveRedlichKwong::rho(
return this->W()/molarVolume;
}
//- Return density [kg/m^3]on
//- Return density [kg/m^3]on
inline scalar soaveRedlichKwong::rho(const scalar p,const scalar T) const
{
//CL: using perfect gas equation as starting point

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

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd |
\\ / A nd |
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -25,7 +25,7 @@ License
Author
Christian Lucas
Institut für Thermodynamik
Technische Universität Braunschweig
Technische Universität Braunschweig
Germany
\*---------------------------------------------------------------------------*/

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

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd |
\\ / A nd |
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -26,13 +26,13 @@ Class
Foam::nasaHeatCapacityPolynomial
Description
Nasa Heat Capacity Polynomial for perfect Gas (7.order polynominal) --> freely available for many fluids
--> paper title: NASA Glenn Coefficients for Calculating Thermodynamic Properties of Individual Species
Nasa Heat Capacity Polynomial for perfect Gas (7.order polynominal) --> freely available for many fluids
--> paper title: NASA Glenn Coefficients for Calculating Thermodynamic Properties of Individual Species
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
@ -42,7 +42,7 @@ SourceFiles
Author
Christian Lucas
Institut für Thermodynamik
Technische Universität Braunschweig
Technische Universität Braunschweig
Germany
\*---------------------------------------------------------------------------*/
@ -174,7 +174,7 @@ public:
inline static autoPtr<nasaHeatCapacityPolynomial> New(Istream& is);
// Member Functions
//- perfect Gas Enthalpy [J/kmol]
inline scalar h0(const scalar T) const;

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

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd |
\\ / A nd |
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -25,7 +25,7 @@ License
Author
Christian Lucas
Institut für Thermodynamik
Technische Universität Braunschweig
Technische Universität Braunschweig
Germany
\*---------------------------------------------------------------------------*/
@ -62,7 +62,7 @@ inline Foam::nasaHeatCapacityPolynomial<equationOfState>::nasaHeatCapacityPolyno
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()))
cp_std(this->cp_nonLimited(this->rhostd(),this->Tstd()))
{}
//- Construct from components
@ -97,7 +97,7 @@ inline Foam::nasaHeatCapacityPolynomial<equationOfState>::nasaHeatCapacityPolyno
s0_std(s0_std_),
integral_p_dv_std(integral_p_dv_std_),
integral_dpdT_dv_std(integral_dpdT_dv_std_),
cp_std(cp_std_)
cp_std(cp_std_)
{}
@ -149,14 +149,14 @@ Foam::nasaHeatCapacityPolynomial<equationOfState>::New(Istream& is)
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
//used to calculate the internal energy
//perfect gas enthalpy
//perfect gas enthalpy
template<class equationOfState>
inline Foam::scalar Foam::nasaHeatCapacityPolynomial<equationOfState>::h0
(
const scalar T
) const
{
return
return
this->RR()*T*
(
-this->a1_*pow(T,-2)
@ -171,7 +171,7 @@ inline Foam::scalar Foam::nasaHeatCapacityPolynomial<equationOfState>::h0
//used to calculate the internal energy
//perfect gas internal energy
//perfect gas internal energy
template<class equationOfState>
inline Foam::scalar Foam::nasaHeatCapacityPolynomial<equationOfState>::e0
(
@ -181,7 +181,7 @@ inline Foam::scalar Foam::nasaHeatCapacityPolynomial<equationOfState>::e0
return this->h0(T) - this->RR()*T;
}
// used to calculate the entropy
// perfect gas entropy
// perfect gas entropy
template<class equationOfState>
inline Foam::scalar Foam::nasaHeatCapacityPolynomial<equationOfState>::s0
(
@ -189,7 +189,7 @@ inline Foam::scalar Foam::nasaHeatCapacityPolynomial<equationOfState>::s0
) const
{
return this->RR()*
(
(
this->a1_*(-1)/(2*pow(T,2))
-this->a2_/T+this->a3_*log(T)
+this->a4_*T
@ -241,14 +241,14 @@ template<class equationOfState>
inline Foam::scalar Foam::nasaHeatCapacityPolynomial<equationOfState>::cp
(
const scalar rho,
const scalar T
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
// 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
// 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,
@ -258,7 +258,7 @@ inline Foam::scalar Foam::nasaHeatCapacityPolynomial<equationOfState>::cp
-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),
@ -267,10 +267,10 @@ template<class equationOfState>
inline Foam::scalar Foam::nasaHeatCapacityPolynomial<equationOfState>::cp_nonLimited
(
const scalar rho,
const scalar T
const scalar T
) const
{
return fabs(this->cv(rho,T)-T*pow((this->dpdT(rho, T)),2)/min(this->dpdv(rho, T),-1));
return fabs(this->cv(rho,T)-T*pow((this->dpdT(rho, T)),2)/min(this->dpdv(rho, T),-1));
}
@ -281,7 +281,7 @@ 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);
@ -293,9 +293,9 @@ inline Foam::scalar Foam::nasaHeatCapacityPolynomial<equationOfState>::cv
//function to calculate real gas enthalpy
template<class equationOfState>
inline Foam::scalar Foam::nasaHeatCapacityPolynomial<equationOfState>::h
(
(
const scalar rho,
const scalar T
const scalar T
) const
{
return this->e(rho,T)+this->p(rho,T)/rho*this->W()-this->Pstd()/this->rhostd()*this->W();
@ -307,12 +307,12 @@ inline Foam::scalar Foam::nasaHeatCapacityPolynomial<equationOfState>::h
// equation: du= cv0 dT +[T*dp/dT -p]dv
template<class equationOfState>
inline Foam::scalar Foam::nasaHeatCapacityPolynomial<equationOfState>::e
(
(
const scalar rho,
const scalar T
const scalar T
) const
{
return
return
(
-this->Tstd()*integral_dpdT_dv_std
+integral_p_dv_std
@ -332,7 +332,7 @@ inline Foam::scalar Foam::nasaHeatCapacityPolynomial<equationOfState>::s
(
const scalar rho,
const scalar T
) const
{
return -integral_dpdT_dv_std
@ -422,7 +422,7 @@ inline Foam::nasaHeatCapacityPolynomial<equationOfState> Foam::operator+
np1.nMoles()/eofs.nMoles()*np1.a5_
+ np2.nMoles()/eofs.nMoles()*np2.a5_,
np1.nMoles()/eofs.nMoles()*np1.a6_
+ np2.nMoles()/eofs.nMoles()*np2.a6_,
+ np2.nMoles()/eofs.nMoles()*np2.a6_,
np1.nMoles()/eofs.nMoles()*np1.a7_
+ np2.nMoles()/eofs.nMoles()*np2.a7_
);
@ -456,7 +456,7 @@ inline Foam::nasaHeatCapacityPolynomial<equationOfState> Foam::operator-
np1.nMoles()/eofs.nMoles()*np1.a5_
- np2.nMoles()/eofs.nMoles()*np2.a5_,
np1.nMoles()/eofs.nMoles()*np1.a6_
- np2.nMoles()/eofs.nMoles()*np2.a6_,
- np2.nMoles()/eofs.nMoles()*np2.a6_,
np1.nMoles()/eofs.nMoles()*np1.a7_
- np2.nMoles()/eofs.nMoles()*np2.a7_
);
@ -470,8 +470,8 @@ inline Foam::nasaHeatCapacityPolynomial<equationOfState> Foam::operator*
const nasaHeatCapacityPolynomial<equationOfState>& np
)
{
//CL: values at STD don't need to be recalculated,
//CL: therefore, providing the values in the constructor
//CL: values at STD don't need to be recalculated,
//CL: therefore, providing the values in the constructor
return nasaHeatCapacityPolynomial<equationOfState>
(
s*static_cast<const equationOfState&>(np),

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

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd |
\\ / A nd |
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -25,7 +25,7 @@ License
Author
Christian Lucas
Institut für Thermodynamik
Technische Universität Braunschweig
Technische Universität Braunschweig
Germany
\*---------------------------------------------------------------------------*/

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

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd |
\\ / A nd |
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -27,9 +27,9 @@ Class
Description
This is an subversion of the specieThermo function. While the "orginal" specieThermo function
This is an subversion of the specieThermo function. While the "orginal" specieThermo function
is only valid for perfect Gases, this function is valid for real gases.
Basic thermodynamics type based on the use of fitting functions for
cp, h, s obtained from the template argument type thermo. All other
properties are derived from these primitive functions.
@ -43,7 +43,7 @@ SourceFiles
Author
Christian Lucas
Institut für Thermodynamik
Technische Universität Braunschweig
Technische Universität Braunschweig
Germany
@ -157,7 +157,7 @@ public:
// Calculate and return derived properties
// (These functions need not provided in derived types)
//CL: isentropic expansion factor "gamma" (heat capacity ratio for perfect gas)
//CL: isentropic expansion factor "gamma" (heat capacity ratio for perfect gas)
inline scalar gamma(const scalar T, const scalar rho) const;
//- Sensible internal energy [J/kmol]
@ -201,16 +201,16 @@ public:
//CL: Other variables
//- Return compressibility drho/dp at h=constant [s^2/m^2]
//- Return compressibility drho/dp at h=constant [s^2/m^2]
inline scalar psiH(const scalar rho, const scalar T) const;
//- Return compressibility drho/dp at e=constant [s^2/m^2]
//- Return compressibility drho/dp at e=constant [s^2/m^2]
inline scalar psiE(const scalar rho, const scalar T) const;
//- Return compressibility drho/dH at p=constant
//- Return compressibility drho/dH at p=constant
inline scalar drhodH(const scalar rho, const scalar T) const;
//- Return compressibility drho/dE at p=constant
//- Return compressibility drho/dE at p=constant
inline scalar drhodE(const scalar rho, const scalar T) const;
// Energy->temperature inversion functions

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

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd |
\\ / A nd |
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -25,7 +25,7 @@ License
Author
Christian Lucas
Institut für Thermodynamik
Technische Universität Braunschweig
Technische Universität Braunschweig
Germany
@ -75,18 +75,18 @@ inline void Foam::realGasSpecieThermo<thermo>::T
i=0;
do
{
//CL: using a stabilizing newton solver
{
//CL: using a stabilizing newton solver
//CL: if the solve is diverging, the step is reduced until the solver converges
Tnew = Test - ((this->*F)(rho,Test) - f)/(this->*dFdT)(rho,Test)/(pow(2,i));
Tnew = Test - ((this->*F)(rho,Test) - f)/(this->*dFdT)(rho,Test)/(pow(2,i));
i++;
}while
(
(i<20)
&&
((
mag((this->*F)(rho,Tnew) - f)
>
mag((this->*F)(rho,Tnew) - f)
>
mag((this->*F)(rho,Test) - f)
))
);
@ -104,7 +104,7 @@ inline void Foam::realGasSpecieThermo<thermo>::T
}
} while
//CL: both fields must converge
(
(
(mag(mag(Tnew) - mag(Test)) > Ttol)
||
(mag(mag(rho) - mag(rhoOld)) > rhotol)
@ -206,12 +206,12 @@ inline Foam::scalar Foam::realGasSpecieThermo<thermo>::psiH
(
const scalar rho,
const scalar T
) const
{
scalar beta=this->isobarExpCoef(rho,T);
return
return
-(
(T*beta*beta-beta)/this->Cp(rho,T)
-this->isothermalCompressiblity(rho,T)*rho
@ -225,14 +225,14 @@ inline Foam::scalar Foam::realGasSpecieThermo<thermo>::psiE
(
const scalar rho,
const scalar T
) const
{
scalar V = 1/rho;
scalar cp=this->Cp(rho,T);
scalar beta=this->isobarExpCoef(rho,T);
return
return
-(
(
T*pow(beta,2)*V
@ -246,27 +246,27 @@ inline Foam::scalar Foam::realGasSpecieThermo<thermo>::psiE
);
}
//CL:- Returns drho/dH at p=constant
//CL:- Returns drho/dH at p=constant
//CL:- using Bridgeman's Table
template<class thermo>
inline Foam::scalar Foam::realGasSpecieThermo<thermo>::drhodH
(
const scalar rho,
const scalar T
) const
{
return -(rho*this->isobarExpCoef(rho,T))/this->Cp(rho,T);
}
//CL:- Returns drho/dE at p=constant
//CL:- Returns drho/dE at p=constant
//CL:- using Bridgeman's Table
template<class thermo>
inline Foam::scalar Foam::realGasSpecieThermo<thermo>::drhodE
(
const scalar rho,
const scalar T
) const
{
scalar beta=this->isobarExpCoef(rho,T);