[shock77]

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:

template<class Specie>
Foam::PengRobinsonGas<Specie>::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:

// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //

template<class Specie>
inline Foam::scalar Foam::PengRobinsonGas<Specie>::rho
(
    scalar p,
    scalar T
) const
{
    scalar Z = this->Z(p, T);
    return p/(Z*this->R()*T);
}


template<class Specie>
inline Foam::scalar Foam::PengRobinsonGas<Specie>::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<class Specie>
inline Foam::scalar Foam::PengRobinsonGas<Specie>::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<class Specie>
inline Foam::scalar Foam::PengRobinsonGas<Specie>::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<class Specie>
inline Foam::scalar Foam::PengRobinsonGas<Specie>::psi
(
    scalar p,
    scalar T
) const
{
    scalar Z = this->Z(p, T);

    return 1.0/(Z*this->R()*T);
}


template<class Specie>
inline Foam::scalar Foam::PengRobinsonGas<Specie>::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<class Specie>
inline Foam::scalar Foam::PengRobinsonGas<Specie>::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