From c8d4346a54ef4ef403401e0faadbbcffe9b9523e Mon Sep 17 00:00:00 2001 From: Henry Weller Date: Sat, 31 Aug 2019 12:04:57 +0100 Subject: [PATCH] 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. --- .../equationOfState/Boussinesq/Boussinesq.H | 5 +- .../equationOfState/Boussinesq/BoussinesqI.H | 56 ++------------- .../PengRobinsonGas/PengRobinsonGas.H | 2 +- .../PengRobinsonGas/PengRobinsonGasI.H | 53 ++++++-------- .../adiabaticPerfectFluid.H | 5 +- .../adiabaticPerfectFluidI.H | 68 ++---------------- .../icoPolynomial/icoPolynomial.H | 5 +- .../icoPolynomial/icoPolynomialI.H | 54 ++------------ .../incompressiblePerfectGas.H | 2 +- .../incompressiblePerfectGasI.H | 66 +++++++---------- .../specie/equationOfState/linear/linear.H | 5 +- .../specie/equationOfState/linear/linearI.H | 60 ++-------------- .../perfectFluid/perfectFluid.H | 11 ++- .../perfectFluid/perfectFluidI.H | 71 ++----------------- .../equationOfState/perfectGas/perfectGas.H | 2 +- .../equationOfState/rhoConst/rhoConst.C | 6 +- .../equationOfState/rhoConst/rhoConst.H | 2 +- .../equationOfState/rhoConst/rhoConstI.H | 55 ++++++-------- .../liquidProperties/liquidProperties.H | 2 +- .../thermophysicalProperties.H | 2 +- .../thermophysicalPropertiesSelector.H | 2 +- 21 files changed, 127 insertions(+), 407 deletions(-) diff --git a/src/thermophysicalModels/specie/equationOfState/Boussinesq/Boussinesq.H b/src/thermophysicalModels/specie/equationOfState/Boussinesq/Boussinesq.H index 02375780ac..076ad5ea03 100644 --- a/src/thermophysicalModels/specie/equationOfState/Boussinesq/Boussinesq.H +++ b/src/thermophysicalModels/specie/equationOfState/Boussinesq/Boussinesq.H @@ -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 [] diff --git a/src/thermophysicalModels/specie/equationOfState/Boussinesq/BoussinesqI.H b/src/thermophysicalModels/specie/equationOfState/Boussinesq/BoussinesqI.H index 74e3a581b3..eeefe7fbb9 100644 --- a/src/thermophysicalModels/specie/equationOfState/Boussinesq/BoussinesqI.H +++ b/src/thermophysicalModels/specie/equationOfState/Boussinesq/BoussinesqI.H @@ -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::operator+= const Boussinesq& 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 Foam::operator+ const Boussinesq& b2 ) { - Specie sp(static_cast(b1) + static_cast(b2)); - - if (mag(sp.Y()) < small) - { - return Boussinesq - ( - sp, - b1.rho0_, - b1.T0_, - b1.beta_ - ); - } - else - { - const scalar Y1 = b1.Y()/sp.Y(); - const scalar Y2 = b2.Y()/sp.Y(); - - return Boussinesq - ( - 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 Foam::operator== const Boussinesq& b2 ) { - Specie sp(static_cast(b1) == static_cast(b2)); - - const scalar Y1 = b1.Y()/sp.Y(); - const scalar Y2 = b2.Y()/sp.Y(); - - return Boussinesq - ( - sp, - Y2*b2.rho0_ - Y1*b1.rho0_, - Y2*b2.T0_ - Y1*b1.T0_, - Y2*b2.beta_ - Y1*b1.beta_ - ); + NotImplemented; + return b1; } diff --git a/src/thermophysicalModels/specie/equationOfState/PengRobinsonGas/PengRobinsonGas.H b/src/thermophysicalModels/specie/equationOfState/PengRobinsonGas/PengRobinsonGas.H index 29805f93fd..88a4ee1fd8 100644 --- a/src/thermophysicalModels/specie/equationOfState/PengRobinsonGas/PengRobinsonGas.H +++ b/src/thermophysicalModels/specie/equationOfState/PengRobinsonGas/PengRobinsonGas.H @@ -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 [] diff --git a/src/thermophysicalModels/specie/equationOfState/PengRobinsonGas/PengRobinsonGasI.H b/src/thermophysicalModels/specie/equationOfState/PengRobinsonGas/PengRobinsonGasI.H index dc16404ba2..9576e11500 100644 --- a/src/thermophysicalModels/specie/equationOfState/PengRobinsonGas/PengRobinsonGasI.H +++ b/src/thermophysicalModels/specie/equationOfState/PengRobinsonGas/PengRobinsonGasI.H @@ -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::operator+= const PengRobinsonGas& 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 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 ( @@ -406,7 +406,7 @@ Foam::PengRobinsonGas 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 Foam::operator== const PengRobinsonGas& pg2 ) { - Specie sp - ( - static_cast(pg1) - == static_cast(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 ( - sp, - Tc, - Vc, - Zc, - RR*Zc*Tc/Vc, - Y2*pg2.omega_ - Y1*pg1.omega_ + static_cast(pg1) == static_cast(pg2), + NaN, + NaN, + NaN, + NaN, + NaN ); } diff --git a/src/thermophysicalModels/specie/equationOfState/adiabaticPerfectFluid/adiabaticPerfectFluid.H b/src/thermophysicalModels/specie/equationOfState/adiabaticPerfectFluid/adiabaticPerfectFluid.H index c6ba0ea890..59ce65bf6a 100644 --- a/src/thermophysicalModels/specie/equationOfState/adiabaticPerfectFluid/adiabaticPerfectFluid.H +++ b/src/thermophysicalModels/specie/equationOfState/adiabaticPerfectFluid/adiabaticPerfectFluid.H @@ -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 [] diff --git a/src/thermophysicalModels/specie/equationOfState/adiabaticPerfectFluid/adiabaticPerfectFluidI.H b/src/thermophysicalModels/specie/equationOfState/adiabaticPerfectFluid/adiabaticPerfectFluidI.H index 8c838244fa..593649a43d 100644 --- a/src/thermophysicalModels/specie/equationOfState/adiabaticPerfectFluid/adiabaticPerfectFluidI.H +++ b/src/thermophysicalModels/specie/equationOfState/adiabaticPerfectFluid/adiabaticPerfectFluidI.H @@ -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::operator+= const adiabaticPerfectFluid& 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 Foam::operator+ const adiabaticPerfectFluid& pf2 ) { - Specie sp - ( - static_cast(pf1) - + static_cast(pf2) - ); - - if (mag(sp.Y()) < small) - { - return adiabaticPerfectFluid - ( - 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 - ( - 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 Foam::operator== const adiabaticPerfectFluid& pf2 ) { - Specie sp - ( - static_cast(pf1) - == static_cast(pf2) - ); - - const scalar Y1 = pf1.Y()/sp.Y(); - const scalar Y2 = pf2.Y()/sp.Y(); - - return adiabaticPerfectFluid - ( - 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; } diff --git a/src/thermophysicalModels/specie/equationOfState/icoPolynomial/icoPolynomial.H b/src/thermophysicalModels/specie/equationOfState/icoPolynomial/icoPolynomial.H index 41a5e887df..e910e1ea84 100644 --- a/src/thermophysicalModels/specie/equationOfState/icoPolynomial/icoPolynomial.H +++ b/src/thermophysicalModels/specie/equationOfState/icoPolynomial/icoPolynomial.H @@ -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 [] diff --git a/src/thermophysicalModels/specie/equationOfState/icoPolynomial/icoPolynomialI.H b/src/thermophysicalModels/specie/equationOfState/icoPolynomial/icoPolynomialI.H index 6a4304a91f..6c3c307a86 100644 --- a/src/thermophysicalModels/specie/equationOfState/icoPolynomial/icoPolynomialI.H +++ b/src/thermophysicalModels/specie/equationOfState/icoPolynomial/icoPolynomialI.H @@ -184,16 +184,7 @@ inline void Foam::icoPolynomial::operator+= const icoPolynomial& 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 Foam::operator+ const icoPolynomial& ip2 ) { - Specie sp - ( - static_cast(ip1) - + static_cast(ip2) - ); - - if (mag(sp.Y()) < small) - { - return icoPolynomial - ( - sp, - ip1.rhoCoeffs_ - ); - } - else - { - const scalar Y1 = ip1.Y()/sp.Y(); - const scalar Y2 = ip2.Y()/sp.Y(); - - return icoPolynomial - ( - sp, - Y1*ip1.rhoCoeffs_ + Y2*ip2.rhoCoeffs_ - ); - } + NotImplemented; + return ip1; } @@ -263,20 +231,8 @@ Foam::icoPolynomial Foam::operator== const icoPolynomial& ip2 ) { - Specie sp - ( - static_cast(ip1) - == static_cast(ip2) - ); - - const scalar Y1 = ip1.Y()/sp.Y(); - const scalar Y2 = ip2.Y()/sp.Y(); - - return icoPolynomial - ( - sp, - Y2*ip2.rhoCoeffs_ - Y1*ip1.rhoCoeffs_ - ); + NotImplemented; + return ip1; } diff --git a/src/thermophysicalModels/specie/equationOfState/incompressiblePerfectGas/incompressiblePerfectGas.H b/src/thermophysicalModels/specie/equationOfState/incompressiblePerfectGas/incompressiblePerfectGas.H index 6b30f00961..1321317b3c 100644 --- a/src/thermophysicalModels/specie/equationOfState/incompressiblePerfectGas/incompressiblePerfectGas.H +++ b/src/thermophysicalModels/specie/equationOfState/incompressiblePerfectGas/incompressiblePerfectGas.H @@ -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 [] diff --git a/src/thermophysicalModels/specie/equationOfState/incompressiblePerfectGas/incompressiblePerfectGasI.H b/src/thermophysicalModels/specie/equationOfState/incompressiblePerfectGas/incompressiblePerfectGasI.H index 2c48cd814f..4486de6272 100644 --- a/src/thermophysicalModels/specie/equationOfState/incompressiblePerfectGas/incompressiblePerfectGasI.H +++ b/src/thermophysicalModels/specie/equationOfState/incompressiblePerfectGas/incompressiblePerfectGasI.H @@ -185,16 +185,17 @@ inline void Foam::incompressiblePerfectGas::operator+= const incompressiblePerfectGas& 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 Foam::operator+ const incompressiblePerfectGas& 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 ( - static_cast(ipg1) - + static_cast(ipg2) + static_cast(ipg1) + static_cast(ipg2), + ipg1.pRef_ ); - - if (mag(sp.Y()) < small) - { - return incompressiblePerfectGas - ( - sp, - ipg1.pRef_ - ); - } - else - { - const scalar Y1 = ipg1.Y()/sp.Y(); - const scalar Y2 = ipg2.Y()/sp.Y(); - - return incompressiblePerfectGas - ( - sp, - Y1*ipg1.pRef_ + Y2*ipg2.pRef_ - ); - } } @@ -264,19 +255,10 @@ inline Foam::incompressiblePerfectGas Foam::operator== const incompressiblePerfectGas& ipg2 ) { - Specie sp - ( - static_cast(ipg1) - == static_cast(ipg2) - ); - - const scalar Y1 = ipg1.Y()/sp.Y(); - const scalar Y2 = ipg2.Y()/sp.Y(); - return incompressiblePerfectGas ( - sp, - Y2*ipg2.pRef_ - Y1*ipg1.pRef_ + static_cast(ipg1) == static_cast(ipg2), + NaN ); } diff --git a/src/thermophysicalModels/specie/equationOfState/linear/linear.H b/src/thermophysicalModels/specie/equationOfState/linear/linear.H index d49103c5b8..c8aac35688 100644 --- a/src/thermophysicalModels/specie/equationOfState/linear/linear.H +++ b/src/thermophysicalModels/specie/equationOfState/linear/linear.H @@ -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 [] diff --git a/src/thermophysicalModels/specie/equationOfState/linear/linearI.H b/src/thermophysicalModels/specie/equationOfState/linear/linearI.H index 6ec30b4793..ced54ac47f 100644 --- a/src/thermophysicalModels/specie/equationOfState/linear/linearI.H +++ b/src/thermophysicalModels/specie/equationOfState/linear/linearI.H @@ -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::operator+= const linear& 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 Foam::operator+ const linear& pf2 ) { - Specie sp - ( - static_cast(pf1) - + static_cast(pf2) - ); - - if (mag(sp.Y()) < small) - { - return linear - ( - sp, - pf1.psi_, - pf1.rho0_ - ); - } - else - { - const scalar Y1 = pf1.Y()/sp.Y(); - const scalar Y2 = pf2.Y()/sp.Y(); - - return linear - ( - sp, - Y1*pf1.psi_ + Y2*pf2.psi_, - Y1*pf1.rho0_ + Y2*pf2.rho0_ - ); - } + NotImplemented; + return pf1; } @@ -233,21 +198,8 @@ inline Foam::linear Foam::operator== const linear& pf2 ) { - Specie sp - ( - static_cast(pf1) - == static_cast(pf2) - ); - - const scalar Y1 = pf1.Y()/sp.Y(); - const scalar Y2 = pf2.Y()/sp.Y(); - - return linear - ( - sp, - Y2*pf2.psi_ - Y1*pf1.psi_, - Y2*pf2.rho0_ - Y1*pf1.rho0_ - ); + NotImplemented; + return pf1; } diff --git a/src/thermophysicalModels/specie/equationOfState/perfectFluid/perfectFluid.H b/src/thermophysicalModels/specie/equationOfState/perfectFluid/perfectFluid.H index 30eaf30ca4..495fd5f79c 100644 --- a/src/thermophysicalModels/specie/equationOfState/perfectFluid/perfectFluid.H +++ b/src/thermophysicalModels/specie/equationOfState/perfectFluid/perfectFluid.H @@ -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 [] diff --git a/src/thermophysicalModels/specie/equationOfState/perfectFluid/perfectFluidI.H b/src/thermophysicalModels/specie/equationOfState/perfectFluid/perfectFluidI.H index 9204872f40..933908c86b 100644 --- a/src/thermophysicalModels/specie/equationOfState/perfectFluid/perfectFluidI.H +++ b/src/thermophysicalModels/specie/equationOfState/perfectFluid/perfectFluidI.H @@ -160,17 +160,7 @@ inline void Foam::perfectFluid::operator+= const perfectFluid& 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 Foam::operator+ const perfectFluid& pf2 ) { - Specie sp - ( - static_cast(pf1) - + static_cast(pf2) - ); - - if (mag(sp.Y()) < small) - { - return perfectFluid - ( - sp, - pf1.R_, - pf1.rho0_ - ); - } - else - { - const scalar Y1 = pf1.Y()/sp.Y(); - const scalar Y2 = pf2.Y()/sp.Y(); - - return perfectFluid - ( - 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 Foam::operator== const perfectFluid& pf2 ) { - Specie sp - ( - static_cast(pf1) - == static_cast(pf2) - ); - - if (mag(sp.Y()) < small) - { - return perfectFluid - ( - 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 - ( - sp, - mag(oneByR) < small ? great : 1/oneByR, - Y2*pf2.rho0_ - Y1*pf1.rho0_ - ); - } + NotImplemented; + return pf1; } diff --git a/src/thermophysicalModels/specie/equationOfState/perfectGas/perfectGas.H b/src/thermophysicalModels/specie/equationOfState/perfectGas/perfectGas.H index dbbc8103bf..5cdd38f380 100644 --- a/src/thermophysicalModels/specie/equationOfState/perfectGas/perfectGas.H +++ b/src/thermophysicalModels/specie/equationOfState/perfectGas/perfectGas.H @@ -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 [] diff --git a/src/thermophysicalModels/specie/equationOfState/rhoConst/rhoConst.C b/src/thermophysicalModels/specie/equationOfState/rhoConst/rhoConst.C index 85b8937640..62d3f71691 100644 --- a/src/thermophysicalModels/specie/equationOfState/rhoConst/rhoConst.C +++ b/src/thermophysicalModels/specie/equationOfState/rhoConst/rhoConst.C @@ -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::write(Ostream& os) const // * * * * * * * * * * * * * * * Ostream Operator * * * * * * * * * * * * * // template -Foam::Ostream& Foam::operator<<(Ostream& os, const rhoConst& ico) +Foam::Ostream& Foam::operator<<(Ostream& os, const rhoConst& rc) { - ico.write(os); + rc.write(os); return os; } diff --git a/src/thermophysicalModels/specie/equationOfState/rhoConst/rhoConst.H b/src/thermophysicalModels/specie/equationOfState/rhoConst/rhoConst.H index 0c8af21904..3725802afd 100644 --- a/src/thermophysicalModels/specie/equationOfState/rhoConst/rhoConst.H +++ b/src/thermophysicalModels/specie/equationOfState/rhoConst/rhoConst.H @@ -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 [] diff --git a/src/thermophysicalModels/specie/equationOfState/rhoConst/rhoConstI.H b/src/thermophysicalModels/specie/equationOfState/rhoConst/rhoConstI.H index 3373614b2f..d32d08cdb4 100644 --- a/src/thermophysicalModels/specie/equationOfState/rhoConst/rhoConstI.H +++ b/src/thermophysicalModels/specie/equationOfState/rhoConst/rhoConstI.H @@ -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 inline Foam::rhoConst::rhoConst ( const word& name, - const rhoConst& ico + const rhoConst& rc ) : - Specie(name, ico), - rho_(ico.rho_) + Specie(name, rc), + rho_(rc.rho_) {} @@ -129,17 +129,14 @@ inline Foam::scalar Foam::rhoConst::CpMCv(scalar p, scalar T) const // * * * * * * * * * * * * * * * Member Operators * * * * * * * * * * * * * // template -inline void Foam::rhoConst::operator+=(const rhoConst& ico) +inline void Foam::rhoConst::operator+=(const rhoConst& 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::operator*=(const scalar s) template inline Foam::rhoConst Foam::operator+ ( - const rhoConst& ico1, - const rhoConst& ico2 + const rhoConst& rc1, + const rhoConst& rc2 ) { Specie sp ( - static_cast(ico1) - + static_cast(ico2) + static_cast(rc1) + + static_cast(rc2) ); if (mag(sp.Y()) < small) @@ -171,18 +168,15 @@ inline Foam::rhoConst Foam::operator+ return rhoConst ( sp, - ico1.rho_ + rc1.rho_ ); } else { - const scalar Y1 = ico1.Y()/sp.Y(); - const scalar Y2 = ico2.Y()/sp.Y(); - return rhoConst ( sp, - Y1*ico1.rho_ + Y2*ico2.rho_ + sp.Y()/(rc1.Y()/rc1.rho_ + rc2.Y()/rc2.rho_) ); } } @@ -192,33 +186,24 @@ template inline Foam::rhoConst Foam::operator* ( const scalar s, - const rhoConst& ico + const rhoConst& rc ) { - return rhoConst(s*static_cast(ico), ico.rho_); + return rhoConst(s*static_cast(rc), rc.rho_); } template inline Foam::rhoConst Foam::operator== ( - const rhoConst& ico1, - const rhoConst& ico2 + const rhoConst& rc1, + const rhoConst& rc2 ) { - Specie sp - ( - static_cast(ico1) - == static_cast(ico2) - ); - - const scalar Y1 = ico1.Y()/sp.Y(); - const scalar Y2 = ico2.Y()/sp.Y(); - return rhoConst ( - sp, - Y2*ico2.rho_ - Y1*ico1.rho_ + static_cast(rc1) == static_cast(rc2), + NaN ); } diff --git a/src/thermophysicalModels/thermophysicalProperties/liquidProperties/liquidProperties/liquidProperties.H b/src/thermophysicalModels/thermophysicalProperties/liquidProperties/liquidProperties/liquidProperties.H index f26e226883..016054d3ae 100644 --- a/src/thermophysicalModels/thermophysicalProperties/liquidProperties/liquidProperties/liquidProperties.H +++ b/src/thermophysicalModels/thermophysicalProperties/liquidProperties/liquidProperties/liquidProperties.H @@ -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; diff --git a/src/thermophysicalModels/thermophysicalProperties/thermophysicalProperties/thermophysicalProperties.H b/src/thermophysicalModels/thermophysicalProperties/thermophysicalProperties/thermophysicalProperties.H index 57912a6bb6..741cef0819 100644 --- a/src/thermophysicalModels/thermophysicalProperties/thermophysicalProperties/thermophysicalProperties.H +++ b/src/thermophysicalModels/thermophysicalProperties/thermophysicalProperties/thermophysicalProperties.H @@ -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; diff --git a/src/thermophysicalModels/thermophysicalProperties/thermophysicalPropertiesSelector/thermophysicalPropertiesSelector.H b/src/thermophysicalModels/thermophysicalProperties/thermophysicalPropertiesSelector/thermophysicalPropertiesSelector.H index 1842ca0e39..c606b6e121 100644 --- a/src/thermophysicalModels/thermophysicalProperties/thermophysicalPropertiesSelector/thermophysicalPropertiesSelector.H +++ b/src/thermophysicalModels/thermophysicalProperties/thermophysicalPropertiesSelector/thermophysicalPropertiesSelector.H @@ -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;