pr 状态方程的另一种求解方法
-
template<class Specie> inline Foam::scalar Foam::PengRobinson<Specie>::rho ( scalar p, scalar T ) const { scalar b = 0.07780*this->R()*Tc_/Pc_; scalar a = 0.45724*sqr(this->R())*sqr(Tc_)/Pc_; scalar Tr = T/Tc_; scalar alpha = sqr ( 1.0 + (0.37464 + 1.54226*omega_- 0.26992*sqr(omega_)) * (1.0 - sqrt(Tr)) ); scalar c = (alpha*a);//a(T) scalar V0 = b; do { d = sqr(V0) + 2*b*V0 - sqr(b); V = b + (this->R()*T - p)*(d/c); V0 = V; } while ( fabs(V - V0) < 1e-6 ); return (this->W()/V0); }