CFD Online Logo CFD Online URL
www.cfd-online.com
[Sponsors]
Home > Forums > Software User Forums > OpenFOAM > OpenFOAM Running, Solving & CFD

Peng Robinson Equation of state

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

Reply
 
LinkBack Thread Tools Search this Thread Display Modes
Old   November 18, 2021, 05:48
Default Peng Robinson Equation of state
  #1
Senior Member
 
Join Date: Dec 2019
Posts: 203
Rep Power: 6
shock77 is on a distinguished road
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<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:


Code:
// * * * * * * * * * * * * * * * 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
shock77 is offline   Reply With Quote

Reply

Thread Tools Search this Thread
Search this Thread:

Advanced Search
Display Modes

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off
Trackbacks are Off
Pingbacks are On
Refbacks are On


Similar Threads
Thread Thread Starter Forum Replies Last Post
Writing an equation of state FAKHREDDINE CFX 2 April 3, 2017 13:08
Simple piston movement in cylinder- fluid models arun1994 CFX 4 July 8, 2016 03:54
continuity equation was diverging in transient state stenber FLUENT 0 March 19, 2016 23:37
How can I add an equation of state for rho and cp ?? Canesin OpenFOAM Programming & Development 0 May 14, 2010 18:52
Real gas equation of state for gas mixture ahourri FLUENT 0 April 30, 2008 12:23


All times are GMT -4. The time now is 09:04.