Equations of state: Corrected coefficient mixing

and removed it from equations of state for which it is too inaccurate to be
useful, in particular from perfectFluid which has been superseded by rPolynomial
which is more accurate and naturally support coefficient mixing.
This commit is contained in:
Henry Weller
2019-08-31 12:04:57 +01:00
parent 4817971e13
commit c8d4346a54
21 changed files with 127 additions and 407 deletions

View File

@ -32,6 +32,9 @@ Description
rho = rho0*(1 - beta*(T - T0))
\endverbatim
Coefficient mixing is very inaccurate and not supported,
so this equation of state is not applicable to mixtures.
SourceFiles
BoussinesqI.H
Boussinesq.C
@ -172,7 +175,7 @@ public:
//- Return entropy [J/kg/K]
inline scalar S(const scalar p, const scalar T) const;
//- Return compressibility rho/p [s^2/m^2]
//- Return compressibility [s^2/m^2]
inline scalar psi(scalar p, scalar T) const;
//- Return compression factor []

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Copyright (C) 2015-2018 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2015-2019 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -187,18 +187,7 @@ inline void Foam::Boussinesq<Specie>::operator+=
const Boussinesq<Specie>& b
)
{
scalar Y1 = this->Y();
Specie::operator+=(b);
if (mag(this->Y()) > small)
{
Y1 /= this->Y();
const scalar Y2 = b.Y()/this->Y();
rho0_ = Y1*rho0_ + Y2*b.rho0_;
T0_ = Y1*T0_ + Y2*b.T0_;
beta_ = Y1*beta_ + Y2*b.beta_;
}
NotImplemented;
}
@ -218,31 +207,8 @@ inline Foam::Boussinesq<Specie> Foam::operator+
const Boussinesq<Specie>& b2
)
{
Specie sp(static_cast<const Specie&>(b1) + static_cast<const Specie&>(b2));
if (mag(sp.Y()) < small)
{
return Boussinesq<Specie>
(
sp,
b1.rho0_,
b1.T0_,
b1.beta_
);
}
else
{
const scalar Y1 = b1.Y()/sp.Y();
const scalar Y2 = b2.Y()/sp.Y();
return Boussinesq<Specie>
(
sp,
Y1*b1.rho0_ + Y2*b2.rho0_,
Y1*b1.T0_ + Y2*b2.T0_,
Y1*b1.beta_ + Y2*b2.beta_
);
}
NotImplemented;
return b1;
}
@ -270,18 +236,8 @@ inline Foam::Boussinesq<Specie> Foam::operator==
const Boussinesq<Specie>& b2
)
{
Specie sp(static_cast<const Specie&>(b1) == static_cast<const Specie&>(b2));
const scalar Y1 = b1.Y()/sp.Y();
const scalar Y2 = b2.Y()/sp.Y();
return Boussinesq<Specie>
(
sp,
Y2*b2.rho0_ - Y1*b1.rho0_,
Y2*b2.T0_ - Y1*b1.T0_,
Y2*b2.beta_ - Y1*b1.beta_
);
NotImplemented;
return b1;
}

View File

@ -170,7 +170,7 @@ public:
//- Return entropy [J/kg/K]
inline scalar S(const scalar p, const scalar T) const;
//- Return compressibility rho/p [s^2/m^2]
//- Return compressibility [s^2/m^2]
inline scalar psi(scalar p, scalar T) const;
//- Return compression factor []

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Copyright (C) 2014-2018 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2014-2019 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -338,19 +338,19 @@ inline void Foam::PengRobinsonGas<Specie>::operator+=
const PengRobinsonGas<Specie>& pg
)
{
scalar Y1 = this->Y();
scalar X1 = this->Y()/this->W();
Specie::operator+=(pg);
if (mag(this->Y()) > small)
{
Y1 /= this->Y();
const scalar Y2 = pg.Y()/this->Y();
X1 *= this->W()/this->Y();
const scalar X2 = this->W()*pg.Y()/(pg.W()*this->Y());
Tc_ = Y1*Tc_ + Y2*pg.Tc_;
Vc_ = Y1*Vc_ + Y2*pg.Vc_;
Zc_ = Y1*Zc_ + Y2*pg.Zc_;
Tc_ = X1*Tc_ + X2*pg.Tc_;
Vc_ = X1*Vc_ + X2*pg.Vc_;
Zc_ = X1*Zc_ + X2*pg.Zc_;
Pc_ = RR*Zc_*Tc_/Vc_;
omega_ = Y1*omega_ + Y2*pg.omega_;
omega_ = X1*omega_ + X2*pg.omega_;
}
}
@ -392,12 +392,12 @@ Foam::PengRobinsonGas<Specie> Foam::operator+
}
else
{
const scalar Y1 = pg1.Y()/sp.Y();
const scalar Y2 = pg2.Y()/sp.Y();
const scalar X1 = sp.W()*pg1.Y()/(pg1.W()*sp.Y());
const scalar X2 = sp.W()*pg2.Y()/(pg2.W()*sp.Y());
const scalar Tc = Y1*pg1.Tc_ + Y2*pg2.Tc_;
const scalar Vc = Y1*pg1.Vc_ + Y2*pg2.Vc_;
const scalar Zc = Y1*pg1.Zc_ + Y2*pg2.Zc_;
const scalar Tc = X1*pg1.Tc_ + X2*pg2.Tc_;
const scalar Vc = X1*pg1.Vc_ + X2*pg2.Vc_;
const scalar Zc = X1*pg1.Zc_ + X2*pg2.Zc_;
return PengRobinsonGas<Specie>
(
@ -406,7 +406,7 @@ Foam::PengRobinsonGas<Specie> Foam::operator+
Vc,
Zc,
RR*Zc*Tc/Vc,
Y1*pg1.omega_ + Y2*pg2.omega_
X1*pg1.omega_ + X2*pg2.omega_
);
}
}
@ -438,27 +438,14 @@ Foam::PengRobinsonGas<Specie> Foam::operator==
const PengRobinsonGas<Specie>& pg2
)
{
Specie sp
(
static_cast<const Specie&>(pg1)
== static_cast<const Specie&>(pg2)
);
const scalar Y1 = pg1.Y()/sp.Y();
const scalar Y2 = pg2.Y()/sp.Y();
const scalar Tc = Y2*pg2.Tc_ - Y1*pg1.Tc_;
const scalar Vc = Y2*pg2.Vc_ - Y1*pg1.Vc_;
const scalar Zc = Y2*pg2.Zc_ - Y1*pg1.Zc_;
return PengRobinsonGas<Specie>
(
sp,
Tc,
Vc,
Zc,
RR*Zc*Tc/Vc,
Y2*pg2.omega_ - Y1*pg1.omega_
static_cast<const Specie&>(pg1) == static_cast<const Specie&>(pg2),
NaN,
NaN,
NaN,
NaN,
NaN
);
}

View File

@ -27,6 +27,9 @@ Class
Description
Adiabatic perfect fluid equation of state.
Coefficient mixing is very inaccurate and not supported,
so this equation of state is not applicable to mixtures.
SourceFiles
adiabaticPerfectFluidI.H
adiabaticPerfectFluid.C
@ -169,7 +172,7 @@ public:
//- Return entropy [J/kg/K]
inline scalar S(const scalar p, const scalar T) const;
//- Return compressibility rho/p [s^2/m^2]
//- Return compressibility [s^2/m^2]
inline scalar psi(scalar p, scalar T) const;
//- Return compression factor []

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Copyright (C) 2013-2018 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2013-2019 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -199,19 +199,7 @@ inline void Foam::adiabaticPerfectFluid<Specie>::operator+=
const adiabaticPerfectFluid<Specie>& pf
)
{
scalar Y1 = this->Y();
Specie::operator+=(pf);
if (mag(this->Y()) > small)
{
Y1 /= this->Y();
const scalar Y2 = pf.Y()/this->Y();
p0_ = Y1*p0_ + Y2*pf.p0_;
rho0_ = Y1*rho0_ + Y2*pf.rho0_;
gamma_ = Y1*gamma_ + Y2*pf.gamma_;
B_ = Y1*B_ + Y2*pf.B_;
}
NotImplemented;
}
@ -231,37 +219,8 @@ inline Foam::adiabaticPerfectFluid<Specie> Foam::operator+
const adiabaticPerfectFluid<Specie>& pf2
)
{
Specie sp
(
static_cast<const Specie&>(pf1)
+ static_cast<const Specie&>(pf2)
);
if (mag(sp.Y()) < small)
{
return adiabaticPerfectFluid<Specie>
(
sp,
pf1.p0_,
pf1.rho0_,
pf1.gamma_,
pf1.B_
);
}
else
{
const scalar Y1 = pf1.Y()/sp.Y();
const scalar Y2 = pf2.Y()/sp.Y();
return adiabaticPerfectFluid<Specie>
(
sp,
Y1*pf1.p0_ + Y2*pf2.p0_,
Y1*pf1.rho0_ + Y2*pf2.rho0_,
Y1*pf1.gamma_ + Y2*pf2.gamma_,
Y1*pf1.B_ + Y2*pf2.B_
);
}
NotImplemented;
return pf1;
}
@ -290,23 +249,8 @@ inline Foam::adiabaticPerfectFluid<Specie> Foam::operator==
const adiabaticPerfectFluid<Specie>& pf2
)
{
Specie sp
(
static_cast<const Specie&>(pf1)
== static_cast<const Specie&>(pf2)
);
const scalar Y1 = pf1.Y()/sp.Y();
const scalar Y2 = pf2.Y()/sp.Y();
return adiabaticPerfectFluid<Specie>
(
sp,
Y2*pf2.p0_ - Y1*pf1.p0_,
Y2*pf2.rho0_ - Y1*pf1.rho0_,
Y2*pf2.gamma_ - Y1*pf1.gamma_,
Y2*pf2.B_ - Y1*pf1.B_
);
NotImplemented;
return pf1;
}

View File

@ -28,6 +28,9 @@ Description
Incompressible, polynomial form of equation of state, using a polynomial
function for density.
Coefficient mixing is very inaccurate and not supported,
so this equation of state is not applicable to mixtures.
Usage
\table
Property | Description
@ -176,7 +179,7 @@ public:
//- Return entropy [J/kg/K]
inline scalar S(const scalar p, const scalar T) const;
//- Return compressibility rho/p [s^2/m^2]
//- Return compressibility [s^2/m^2]
inline scalar psi(scalar p, scalar T) const;
//- Return compression factor []

View File

@ -184,16 +184,7 @@ inline void Foam::icoPolynomial<Specie, PolySize>::operator+=
const icoPolynomial<Specie, PolySize>& ip
)
{
scalar Y1 = this->Y();
Specie::operator+=(ip);
if (mag(this->Y()) > small)
{
Y1 /= this->Y();
const scalar Y2 = ip.Y()/this->Y();
rhoCoeffs_ = Y1*rhoCoeffs_ + Y2*ip.rhoCoeffs_;
}
NotImplemented;
}
@ -213,31 +204,8 @@ Foam::icoPolynomial<Specie, PolySize> Foam::operator+
const icoPolynomial<Specie, PolySize>& ip2
)
{
Specie sp
(
static_cast<const Specie&>(ip1)
+ static_cast<const Specie&>(ip2)
);
if (mag(sp.Y()) < small)
{
return icoPolynomial<Specie, PolySize>
(
sp,
ip1.rhoCoeffs_
);
}
else
{
const scalar Y1 = ip1.Y()/sp.Y();
const scalar Y2 = ip2.Y()/sp.Y();
return icoPolynomial<Specie, PolySize>
(
sp,
Y1*ip1.rhoCoeffs_ + Y2*ip2.rhoCoeffs_
);
}
NotImplemented;
return ip1;
}
@ -263,20 +231,8 @@ Foam::icoPolynomial<Specie, PolySize> Foam::operator==
const icoPolynomial<Specie, PolySize>& ip2
)
{
Specie sp
(
static_cast<const Specie&>(ip1)
== static_cast<const Specie&>(ip2)
);
const scalar Y1 = ip1.Y()/sp.Y();
const scalar Y2 = ip2.Y()/sp.Y();
return icoPolynomial<Specie, PolySize>
(
sp,
Y2*ip2.rhoCoeffs_ - Y1*ip1.rhoCoeffs_
);
NotImplemented;
return ip1;
}

View File

@ -157,7 +157,7 @@ public:
//- Return entropy [J/kg/K]
inline scalar S(const scalar p, const scalar T) const;
//- Return compressibility rho/p [s^2/m^2]
//- Return compressibility [s^2/m^2]
inline scalar psi(scalar p, scalar T) const;
//- Return compression factor []

View File

@ -185,16 +185,17 @@ inline void Foam::incompressiblePerfectGas<Specie>::operator+=
const incompressiblePerfectGas<Specie>& ipg
)
{
scalar Y1 = this->Y();
Specie::operator+=(ipg);
if (mag(this->Y()) > small)
if (notEqual(pRef_, ipg.pRef_))
{
Y1 /= this->Y();
const scalar Y2 = ipg.Y()/this->Y();
pRef_ = Y1*pRef_ + Y2*ipg.pRef_;
FatalErrorInFunction
<< "pRef " << pRef_ << " for "
<< (this->name().size() ? this->name() : "others")
<< " != " << ipg.pRef_ << " for "
<< (ipg.name().size() ? ipg.name() : "others")
<< exit(FatalError);
}
Specie::operator+=(ipg);
}
@ -214,31 +215,21 @@ inline Foam::incompressiblePerfectGas<Specie> Foam::operator+
const incompressiblePerfectGas<Specie>& ipg2
)
{
Specie sp
if (notEqual(ipg1.pRef_, ipg2.pRef_))
{
FatalErrorInFunction
<< "pRef " << ipg1.pRef_ << " for "
<< (ipg1.name().size() ? ipg1.name() : "others")
<< " != " << ipg2.pRef_ << " for "
<< (ipg2.name().size() ? ipg2.name() : "others")
<< exit(FatalError);
}
return incompressiblePerfectGas<Specie>
(
static_cast<const Specie&>(ipg1)
+ static_cast<const Specie&>(ipg2)
static_cast<const Specie&>(ipg1) + static_cast<const Specie&>(ipg2),
ipg1.pRef_
);
if (mag(sp.Y()) < small)
{
return incompressiblePerfectGas<Specie>
(
sp,
ipg1.pRef_
);
}
else
{
const scalar Y1 = ipg1.Y()/sp.Y();
const scalar Y2 = ipg2.Y()/sp.Y();
return incompressiblePerfectGas<Specie>
(
sp,
Y1*ipg1.pRef_ + Y2*ipg2.pRef_
);
}
}
@ -264,19 +255,10 @@ inline Foam::incompressiblePerfectGas<Specie> Foam::operator==
const incompressiblePerfectGas<Specie>& ipg2
)
{
Specie sp
(
static_cast<const Specie&>(ipg1)
== static_cast<const Specie&>(ipg2)
);
const scalar Y1 = ipg1.Y()/sp.Y();
const scalar Y2 = ipg2.Y()/sp.Y();
return incompressiblePerfectGas<Specie>
(
sp,
Y2*ipg2.pRef_ - Y1*ipg1.pRef_
static_cast<const Specie&>(ipg1) == static_cast<const Specie&>(ipg2),
NaN
);
}

View File

@ -31,6 +31,9 @@ Description
rho = rho0 + psi*p
\endverbatim
Coefficient mixing is very inaccurate and not supported,
so this equation of state is not applicable to mixtures.
SourceFiles
linearI.H
linear.C
@ -158,7 +161,7 @@ public:
//- Return entropy [J/kg/K]
inline scalar S(const scalar p, const scalar T) const;
//- Return compressibility rho/p [s^2/m^2]
//- Return compressibility [s^2/m^2]
inline scalar psi(scalar p, scalar T) const;
//- Return compression factor []

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Copyright (C) 2013-2018 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2013-2019 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -150,17 +150,7 @@ inline void Foam::linear<Specie>::operator+=
const linear<Specie>& pf
)
{
scalar Y1 = this->Y();
Specie::operator+=(pf);
if (mag(this->Y()) > small)
{
Y1 /= this->Y();
const scalar Y2 = pf.Y()/this->Y();
psi_ = Y1*psi_ + Y2*pf.psi_;
rho0_ = Y1*rho0_ + Y2*pf.rho0_;
}
NotImplemented;
}
@ -180,33 +170,8 @@ inline Foam::linear<Specie> Foam::operator+
const linear<Specie>& pf2
)
{
Specie sp
(
static_cast<const Specie&>(pf1)
+ static_cast<const Specie&>(pf2)
);
if (mag(sp.Y()) < small)
{
return linear<Specie>
(
sp,
pf1.psi_,
pf1.rho0_
);
}
else
{
const scalar Y1 = pf1.Y()/sp.Y();
const scalar Y2 = pf2.Y()/sp.Y();
return linear<Specie>
(
sp,
Y1*pf1.psi_ + Y2*pf2.psi_,
Y1*pf1.rho0_ + Y2*pf2.rho0_
);
}
NotImplemented;
return pf1;
}
@ -233,21 +198,8 @@ inline Foam::linear<Specie> Foam::operator==
const linear<Specie>& pf2
)
{
Specie sp
(
static_cast<const Specie&>(pf1)
== static_cast<const Specie&>(pf2)
);
const scalar Y1 = pf1.Y()/sp.Y();
const scalar Y2 = pf2.Y()/sp.Y();
return linear<Specie>
(
sp,
Y2*pf2.psi_ - Y1*pf1.psi_,
Y2*pf2.rho0_ - Y1*pf1.rho0_
);
NotImplemented;
return pf1;
}

View File

@ -25,7 +25,14 @@ Class
Foam::perfectFluid
Description
Perfect gas equation of state.
Simple extension of the perfect gas equation of state to liquids by
the addition of a constant density off-set.
Coefficient mixing is very inaccurate and not supported,
so this equation of state is not applicable to mixtures.
This equation of state is rather inaccurate and has been superseded by
rPolynomial which is much more accurate and supports mixtures.
SourceFiles
perfectFluidI.H
@ -157,7 +164,7 @@ public:
//- Return entropy [J/kg/K]
inline scalar S(const scalar p, const scalar T) const;
//- Return compressibility rho/p [s^2/m^2]
//- Return compressibility [s^2/m^2]
inline scalar psi(scalar p, scalar T) const;
//- Return compression factor []

View File

@ -160,17 +160,7 @@ inline void Foam::perfectFluid<Specie>::operator+=
const perfectFluid<Specie>& pf
)
{
scalar Y1 = this->Y();
Specie::operator+=(pf);
if (mag(this->Y()) > small)
{
Y1 /= this->Y();
const scalar Y2 = pf.Y()/this->Y();
R_ = 1.0/(Y1/R_ + Y2/pf.R_);
rho0_ = Y1*rho0_ + Y2*pf.rho0_;
}
NotImplemented;
}
@ -190,33 +180,8 @@ inline Foam::perfectFluid<Specie> Foam::operator+
const perfectFluid<Specie>& pf2
)
{
Specie sp
(
static_cast<const Specie&>(pf1)
+ static_cast<const Specie&>(pf2)
);
if (mag(sp.Y()) < small)
{
return perfectFluid<Specie>
(
sp,
pf1.R_,
pf1.rho0_
);
}
else
{
const scalar Y1 = pf1.Y()/sp.Y();
const scalar Y2 = pf2.Y()/sp.Y();
return perfectFluid<Specie>
(
sp,
1.0/(Y1/pf1.R_ + Y2/pf2.R_),
Y1*pf1.rho0_ + Y2*pf2.rho0_
);
}
NotImplemented;
return pf1;
}
@ -243,34 +208,8 @@ inline Foam::perfectFluid<Specie> Foam::operator==
const perfectFluid<Specie>& pf2
)
{
Specie sp
(
static_cast<const Specie&>(pf1)
== static_cast<const Specie&>(pf2)
);
if (mag(sp.Y()) < small)
{
return perfectFluid<Specie>
(
sp,
pf1.R_,
pf1.rho0_
);
}
else
{
const scalar Y1 = pf1.Y()/sp.Y();
const scalar Y2 = pf2.Y()/sp.Y();
const scalar oneByR = Y2/pf2.R_ - Y1/pf1.R_;
return perfectFluid<Specie>
(
sp,
mag(oneByR) < small ? great : 1/oneByR,
Y2*pf2.rho0_ - Y1*pf1.rho0_
);
}
NotImplemented;
return pf1;
}

View File

@ -141,7 +141,7 @@ public:
//- Return entropy [J/kg/K]
inline scalar S(const scalar p, const scalar T) const;
//- Return compressibility rho/p [s^2/m^2]
//- Return compressibility [s^2/m^2]
inline scalar psi(scalar p, scalar T) const;
//- Return compression factor []

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Copyright (C) 2012-2018 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2012-2019 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -53,9 +53,9 @@ void Foam::rhoConst<Specie>::write(Ostream& os) const
// * * * * * * * * * * * * * * * Ostream Operator * * * * * * * * * * * * * //
template<class Specie>
Foam::Ostream& Foam::operator<<(Ostream& os, const rhoConst<Specie>& ico)
Foam::Ostream& Foam::operator<<(Ostream& os, const rhoConst<Specie>& rc)
{
ico.write(os);
rc.write(os);
return os;
}

View File

@ -143,7 +143,7 @@ public:
//- Return entropy [J/kg/K]
inline scalar S(const scalar p, const scalar T) const;
//- Return compressibility rho/p [s^2/m^2]
//- Return compressibility [s^2/m^2]
inline scalar psi(scalar p, scalar T) const;
//- Return compression factor []

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Copyright (C) 2012-2018 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2012-2019 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -45,11 +45,11 @@ template<class Specie>
inline Foam::rhoConst<Specie>::rhoConst
(
const word& name,
const rhoConst<Specie>& ico
const rhoConst<Specie>& rc
)
:
Specie(name, ico),
rho_(ico.rho_)
Specie(name, rc),
rho_(rc.rho_)
{}
@ -129,17 +129,14 @@ inline Foam::scalar Foam::rhoConst<Specie>::CpMCv(scalar p, scalar T) const
// * * * * * * * * * * * * * * * Member Operators * * * * * * * * * * * * * //
template<class Specie>
inline void Foam::rhoConst<Specie>::operator+=(const rhoConst<Specie>& ico)
inline void Foam::rhoConst<Specie>::operator+=(const rhoConst<Specie>& rc)
{
scalar Y1 = this->Y();
Specie::operator+=(ico);
const scalar Y1 = this->Y();
Specie::operator+=(rc);
if (mag(this->Y()) > small)
{
Y1 /= this->Y();
const scalar Y2 = ico.Y()/this->Y();
rho_ = Y1*rho_ + Y2*ico.rho_;
rho_ = this->Y()/(Y1/rho_ + rc.Y()/rc.rho_);
}
}
@ -156,14 +153,14 @@ inline void Foam::rhoConst<Specie>::operator*=(const scalar s)
template<class Specie>
inline Foam::rhoConst<Specie> Foam::operator+
(
const rhoConst<Specie>& ico1,
const rhoConst<Specie>& ico2
const rhoConst<Specie>& rc1,
const rhoConst<Specie>& rc2
)
{
Specie sp
(
static_cast<const Specie&>(ico1)
+ static_cast<const Specie&>(ico2)
static_cast<const Specie&>(rc1)
+ static_cast<const Specie&>(rc2)
);
if (mag(sp.Y()) < small)
@ -171,18 +168,15 @@ inline Foam::rhoConst<Specie> Foam::operator+
return rhoConst<Specie>
(
sp,
ico1.rho_
rc1.rho_
);
}
else
{
const scalar Y1 = ico1.Y()/sp.Y();
const scalar Y2 = ico2.Y()/sp.Y();
return rhoConst<Specie>
(
sp,
Y1*ico1.rho_ + Y2*ico2.rho_
sp.Y()/(rc1.Y()/rc1.rho_ + rc2.Y()/rc2.rho_)
);
}
}
@ -192,33 +186,24 @@ template<class Specie>
inline Foam::rhoConst<Specie> Foam::operator*
(
const scalar s,
const rhoConst<Specie>& ico
const rhoConst<Specie>& rc
)
{
return rhoConst<Specie>(s*static_cast<const Specie&>(ico), ico.rho_);
return rhoConst<Specie>(s*static_cast<const Specie&>(rc), rc.rho_);
}
template<class Specie>
inline Foam::rhoConst<Specie> Foam::operator==
(
const rhoConst<Specie>& ico1,
const rhoConst<Specie>& ico2
const rhoConst<Specie>& rc1,
const rhoConst<Specie>& rc2
)
{
Specie sp
(
static_cast<const Specie&>(ico1)
== static_cast<const Specie&>(ico2)
);
const scalar Y1 = ico1.Y()/sp.Y();
const scalar Y2 = ico2.Y()/sp.Y();
return rhoConst<Specie>
(
sp,
Y2*ico2.rho_ - Y1*ico1.rho_
static_cast<const Specie&>(rc1) == static_cast<const Specie&>(rc2),
NaN
);
}

View File

@ -202,7 +202,7 @@ public:
// Fundamental equation of state properties
//- Liquid compressibility rho/p [s^2/m^2]
//- Liquid compressibility [s^2/m^2]
// Note: currently it is assumed the liquid is incompressible
inline scalar psi(scalar p, scalar T) const;

View File

@ -124,7 +124,7 @@ public:
//- Liquid density [kg/m^3]
virtual scalar rho(scalar p, scalar T) const = 0;
//- Liquid compressibility rho/p [s^2/m^2]
//- Liquid compressibility [s^2/m^2]
// Note: currently it is assumed the liquid is incompressible
virtual scalar psi(scalar p, scalar T) const = 0;

View File

@ -107,7 +107,7 @@ public:
//- Liquid density [kg/m^3]
inline scalar rho(scalar p, scalar T) const;
//- Liquid compressibility rho/p [s^2/m^2]
//- Liquid compressibility [s^2/m^2]
// Note: currently it is assumed the liquid is incompressible
inline scalar psi(scalar p, scalar T) const;