# Peng Robinson Equation of state

 Register Blogs Members List Search Today's Posts Mark Forums Read

 November 18, 2021, 04:48 Peng Robinson Equation of state #1 Senior Member   Join Date: Dec 2019 Posts: 215 Rep Power: 7 Hi everyone, is anyone familiar with the implementation of the Peng Robinson EQS? It is implemented differently than the basic equation. What confuses me first is that in PengRobinsonGas.C a compressibility factor is calculated: Code: ```template Foam::PengRobinsonGas::PengRobinsonGas ( const dictionary& dict ) : Specie(dict), Tc_(readScalar(dict.subDict("equationOfState").lookup("Tc"))), Vc_(readScalar(dict.subDict("equationOfState").lookup("Vc"))), Zc_(1.0), Pc_(readScalar(dict.subDict("equationOfState").lookup("Pc"))), omega_(readScalar(dict.subDict("equationOfState").lookup("omega"))) { Zc_ = Pc_*Vc_/(RR*Tc_); }``` I do not understand why it is needed in the first place. I always thought that for the Peng Robinson equation, if you assume that Pc is infinity, the equation becomes equal to the perfect gas law. Here, it is the opposite. But in PengRobinsonGasl.H there is something similar to the original equation implemented: Code: ```// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // template inline Foam::scalar Foam::PengRobinsonGas::rho ( scalar p, scalar T ) const { scalar Z = this->Z(p, T); return p/(Z*this->R()*T); } template inline Foam::scalar Foam::PengRobinsonGas::h(scalar p, scalar T) const { scalar Pr = p/Pc_; scalar Tr = T/Tc_; scalar B = 0.07780*Pr/Tr; scalar kappa = 0.37464 + 1.54226*omega_ - 0.26992*sqr(omega_); scalar alpha = sqr(1 + kappa*(1 - sqrt(Tr))); scalar Z = this->Z(p, T); return RR*Tc_ *( Tr*(Z - 1) - 2.078*(1 + kappa)*sqrt(alpha) *log((Z + 2.414*B)/(Z - 0.414*B)) ); } template inline Foam::scalar Foam::PengRobinsonGas::cp(scalar p, scalar T) const { scalar Tr = T/Tc_; scalar a = 0.45724*sqr(RR*Tc_)/Pc_; scalar b = 0.07780*RR*Tc_/Pc_; scalar kappa = 0.37464 + 1.54226*omega_ - 0.26992*sqr(omega_); scalar alpha = sqr(1 + kappa*(1 - sqrt(Tr))); scalar A = a*alpha*p/sqr(RR*T); scalar B = b*p/(RR*T); scalar Z = this->Z(p, T); scalar ap = kappa*a*(kappa/Tc_ - (1 + kappa)/sqrt(T*Tc_)); scalar app = kappa*a*(1 + kappa)/(2*sqrt(pow3(T)*Tc_)); scalar M = (sqr(Z) + 2*B*Z - sqr(B))/(Z - B); scalar N = ap*B/(b*RR); const scalar root2 = sqrt(2.0); return app*(T/(2*root2*b))*log((Z + (root2 + 1)*B)/(Z - (root2 - 1)*B)) + RR*sqr(M - N)/(sqr(M) - 2*A*(Z + B)) - RR; } template inline Foam::scalar Foam::PengRobinsonGas::s ( scalar p, scalar T ) const { scalar Pr = p/Pc_; scalar Tr = T/Tc_; scalar B = 0.07780*Pr/Tr; scalar kappa = 0.37464 + 1.54226*omega_ - 0.26992*sqr(omega_); scalar Z = this->Z(p, T); return - RR*log(p/Pstd) + RR *( log(Z - B) - 2.078*kappa*((1 + kappa)/sqrt(Tr) - kappa) *log((Z + 2.414*B)/(Z - 0.414*B)) ); } template inline Foam::scalar Foam::PengRobinsonGas::psi ( scalar p, scalar T ) const { scalar Z = this->Z(p, T); return 1.0/(Z*this->R()*T); } template inline Foam::scalar Foam::PengRobinsonGas::Z ( scalar p, scalar T ) const { scalar Tr = T/Tc_; scalar a = 0.45724*sqr(RR*Tc_)/Pc_; scalar b = 0.07780*RR*Tc_/Pc_; scalar kappa = 0.37464 + 1.54226*omega_ - 0.26992*sqr(omega_); scalar alpha = sqr(1 + kappa*(1 - sqrt(Tr))); scalar A = a*alpha*p/sqr(RR*T); scalar B = b*p/(RR*T); scalar a2 = B - 1; scalar a1 = A - 2*B - 3*sqr(B); scalar a0 = -A*B + sqr(B) + pow3(B); scalar Q = (3*a1 - a2*a2)/9.0; scalar Rl = (9*a2*a1 - 27*a0 - 2*a2*a2*a2)/54.0; scalar Q3 = Q*Q*Q; scalar D = Q3 + Rl*Rl; scalar root = -1; if (D <= 0) { scalar th = ::acos(Rl/sqrt(-Q3)); scalar qm = 2*sqrt(-Q); scalar r1 = qm*cos(th/3.0) - a2/3.0; scalar r2 = qm*cos((th + 2*constant::mathematical::pi)/3.0) - a2/3.0; scalar r3 = qm*cos((th + 4*constant::mathematical::pi)/3.0) - a2/3.0; root = max(r1, max(r2, r3)); } else { // One root is real scalar D05 = sqrt(D); scalar S = pow(Rl + D05, 1.0/3.0); scalar Tl = 0; if (D05 > Rl) { Tl = -pow(mag(Rl - D05), 1.0/3.0); } else { Tl = pow(Rl - D05, 1.0/3.0); } root = S + Tl - a2/3.0; } return root; } template inline Foam::scalar Foam::PengRobinsonGas::cpMcv ( scalar p, scalar T ) const { scalar Tr = T/Tc_; scalar a = 0.45724*sqr(RR*Tc_)/Pc_; scalar b = 0.07780*RR*Tc_/Pc_; scalar kappa = 0.37464 + 1.54226*omega_ - 0.26992*sqr(omega_); scalar alpha = sqr(1 + kappa*(1 - sqrt(Tr))); scalar A = alpha*a*p/sqr(RR*T); scalar B = b*p/(RR*T); scalar Z = this->Z(p, T); scalar ap = kappa*a*(kappa/Tc_ - (1 + kappa)/sqrt(T*Tc_)); scalar M = (sqr(Z) + 2*B*Z - sqr(B))/(Z - B); scalar N = ap*B/(b*RR); return RR*sqr(M - N)/(sqr(M) - 2*A*(Z + B)); }``` I hope some can clearify this to me