CFD中文网

我这样写不对啊？ 不知道哪位可以给看看毛病在哪@administrators @李东岳

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);
}

b = 0.0778RTc*Pc Tc Pc 临界参数
V 是摩尔体积
p t 是温度压力


pr方程可以看做 v 的一元三次函数
通过 b p t 构造迭代函数求解 v
有谁可以指点一下怎么迭代求解V的吗