Skip to main content
[ "article:topic", "showtoc:no", "authorname:lcoleman", "Equation of State", "Peng-Robinson Equation of State" ]
Engineering LibreTexts

5: Equation of State

  • Page ID
    • 0.jpg
    • Contributed by Larry Coleman
    • Engineering Consultant at Consultant on Demand

    Distillation Science (a blend of Chemistry and Chemical Engineering)

    This is Part V, Equation of State of a ten-part series of technical articles on Distillation Science, as is currently practiced on an industrial level. See also Part I, Overview for introductory comments, the scope of the article series, and nomenclature. Part V, Equation of State explains why an Equation of State (EOS) is needed to evaluate the pure-component physical properties that are used in distillation science, reviews the various options and makes a recommendation. Also included are techniques for solving the EOS cubic equations via spreadsheet (i.e., MS Excel) and some work-arounds near the critical point. This article uses the critical constants and acentric factor of Part III. The EOS solutions of this article will feed values to Part VI, Fugacity.

    The general PVT relationship for a fluid is

    \(PV=ZRT \nonumber\)            Equation 5-1

    where the fluid could be a sub-cooled liquid, a saturated liquid (i.e, at its boiling point), a saturated vapor (i.e., at its dewpoint), or a superheated vapor (i.e., a gas). In the ideal situation ( ambient pressure and temperature), ZL is very close to zero for either a sub-cooled or saturated liquid, and  Zv is close to  unity for either a saturated or superheated vapor. The Ideal Gas Law has \(Z_V=1\), so \(ΔZ=(Z_v-Z_L) ≈ 1\nonumber\)

    \(ΔZ\) is an important parameter to Equation 2-1 and 4-1 in previous articles Part II and Part IV. For calculations in this and other Parts, temperature and pressure are in absolute units (°K and atmospheres), and volume in molar units (e.g., cc/mole). In that case, the gas constant \(R\) has units of 82.057 atm-cc/mole-K.

    Van der Waals EOS and more recent improvements

    For a real fluid, \(Z\) is never 0  and almost never =1, but rather a value between 0 and 1, except for gases at high temperatures and pressures (where Zv >1). While the Ideal Gas Law is a nice introductory concept, it is rarely used in distillation science. Instead, a more accurate expression is needed, and so Equation 5-1 is re-written as 

    \(Z=\dfrac{PV}{RT} \nonumber\)               Equation 5-2

    and different models are used to give results: \(Z\) as a function of \(T\) & \(P\). The oldest is the Van der Waals EOS from the 1873 based on a model that molecules are hard bodies that take up space and have inter-particle forces. Note that an Equation of State covers all non-solid phases, so for saturated systems (i.e., fluids at their boiling/condensation point) it will have two real solutions: one for the vapor and one for the liquid: Zand ZL. For a non-saturated system (i.e., a super-heated vapor or sub-cooled liquid) it will have only one real solution.

    Written in terms of Z, the Van der Waals EOS is

    \(Z=\left(P+ \frac{a}{V^2}\right) \times (V-b)/RT \nonumber\)               Equation 5-3

    where a and b are constants derived from just the fluid’s critical properties. While a good general approximation at lower pressures (better than the Ideal Gas Law), this expression becomes less accurate as the critical point is approached, and it always yields 0.375 as Zc. There are very few fluids that have such a high Zc - most are closer to 0.30. Solutions for saturated phase values, \(Z_v\) and \(Z_L\), require manipulation of a cubic equation and solving its roots. However inaccurate it may be, it does allow estimation of the saturated vapor and liquid molar densities: \(\rho_{\nu}=1/V_{\nu} \nonumber \) and \(\rho_{L}=1/V_{L} \nonumber \), as well as \(\Delta Z= Z_{v}-Z_{L} \nonumber \).

    In Part IV, Equation 4‑1 showed how the latent heat of vaporization, \(ΔH_v\), is calculated from the differential Clausius-Clapyeron equation, once a good vapor pressure relationship is known and ΔZ is closely determined. And in Part VI, it will be seen how fugacity coefficients are calculated from an EOS, to partially account for some departures to ideal behavior in mixtures of fluids. 

    In 1949, the Redlich-Kwong EOS improved on Van der Waals, and in 1972 Soave made a modification to yield the S-R-K EOS, by adding in the effect of the acentric factor, ω. Further improvement was done in a general-application EOS by Peng-Robinson in 1976. Since then, there have been many additional EOS proposed for specialized applications which give better results, but require additional data the additional constants. A full history of various EOS’s is given in Wikipedia.

    For purposes of practical chlorosilane distillation design of fugacity (see Part VI), the Peng-Robinson (P-R) EOS appears to give satisfactory results. P-R does not appear to be accurate enough to extend to closely estimating saturated phase densities of chlorosilanes, or latent heats. If enough data were available, it might be possible to determine which of the newer, more specific EOS forms is accurate enough for those derivative properties (i.e., superior to Peng-Robinson). Unfortunately, that data does not yet exist (although a few key data points show the failings of P-R).

    Regardless of the fluid parameters chosen, the P-R EOS yields a Zc of 0.307 at the critical point. While this Zcvalue is still higher than that of chlorosilanes, the accuracy of Zand Zvalues are adequate for fugacity estimation, through 90% of critical pressure (Pc). Since commercial applications only rarely approach 90% of Pc for reasons of process and equipment design, the P-R EOS is an acceptable compromise. It must be pointed out that the P-R EOS should not be used in any calculations that involve differential or integral calculus manipulation of Z, since erroneous (or impossible) values result. 

    Peng-Robinson Equation of State

    \(p= \frac{RT}{V_{m}-b} - \frac{a \alpha}{V_{m}^2 +2bV_{m} -b^2} \nonumber\)                  Equation 5-4


    • \(a=\frac{0.45724R^2T_{c}^2}{p_{c}} \nonumber\)
    • \(b=\frac{0.07780RT_{c}}{p_{c}} \nonumber\)
    • \(\alpha= (1+\kappa(1-T_{r}^{0.5})^2 \nonumber\)
    • \(\kappa=0.37464+1.54226\omega-0.26992\omega^2 \nonumber\)
    • \(T_{r}=\frac{T}{T_{c}} \nonumber\) noting that \(\omega,\kappa, \alpha\) are dimensionless

    In more compact and dimensionless form, the Peng-Robinson EOS can be re-stated as a cubic in \(Z\) using

    \(A=\frac{\alpha a p}{R^2T^2} \nonumber\)   and   \(B=\frac{bp}{RT} \nonumber\)

    \(Z^3 - (1-B)Z^2+(A-2B-3B^2)Z -(AB-B^2-B^3) = 0 \nonumber\)                  Equation 5-5

    At a given temperature, T, and for fluid physical properties Tc, Pc, and ω (from Part III), and where P is the vapor pressure calculated for that fluid by the VP equation of Part IV, dimensionless parameters “A” and “B” are determined for Equation 5-5. This sets up a cubic equation in Z, which has three real roots (as opposed to being imaginary). The largest root is the value of Zv; the smallest root is the value of ZL; and the third root is discarded as having no physical meaning.

    Then knowing Zv and ZL, the value of ΔZ can be calculated, from which the latent heat of vaporization (at that T & P) could be calculated using Equation 4-1. Also from Zv and ZL, the saturated phase densities can be calculated as estimates. But again, the P-R EOS will not hold up under this much mathematical manipulation to give any more than rough estimates of these distillation design properties; and so it is recommended to only use the P-R EOS to estimate fugacities (unless rough estimates of latent heat and phase densities are acceptable) . The P-R EOS is offered to lay the groundwork for future improvement.

    One stumbling block to many is that solving the EOS equation for Z means solving a cubic equation. Currently there are no MS Excel functions for solving the roots of a cubic equation. However a ready solution is found in an older CRC “Standard Math Tables” book, which can be solved with a spreadsheet or programmed into a macro.

    Cubic Equations 

    Method for solving for the roots of a cubic equation  in the variable "y", with roots y1, y2 and y3                        Equation 5-6

    \(y^3 +py^2 + qy + r =0 \nonumber\) may be transformed to the format , \(x^3+ax+b =0 \)

    by substituting for \(y\) the value, \(x - \frac{p}{3} \nonumber\). Here \(a= \frac{1}{3}(3q-p^2) \) and \(b=\frac{1}/{2\tau}(2p^3 - 9pq + 27r \)

    If \(p\), \(q\), and \(r\) are real, then

    • If \(\frac{b^2}{4} +\frac{a^3}{27} > 0 \nonumber\), there will be one real root and two conjugate imaginary roots.
    • If \(\frac{b^2}{4} +\frac{a^3}{27} = 0 \nonumber\), there will be three real roots of which at least two are equal.
    • If \(\frac{b^2}{4} +\frac{a^3}{27} < 0 \nonumber\), there will be three real roots and unequal roots.

    In the last case, a trigonometric solution is useful for solving EOS equations. Compute the value of the angle \(\phi \nonumber\) in the expression 

    \(cos\phi = \frac{-b}{2} \div \sqrt{(\frac{-a^3}{27})}\)

    Then the three transformed roots x1, x2, and x3 will have the following values:

    \(x_1= 2 \sqrt{\frac{-a}{3}cos(\frac{\phi}{3})} \)        and so  cubic root    \(y_1= -\frac{p}{3} +x_1 \nonumber\)

    \(x_2 =2 \sqrt{\frac{-a}{3}cos(\frac{\phi}{3}+120º)} \)      and cubic root     \(y_2= -\frac{p}{3} +x_2 \nonumber\)

    \(x_3=2 \sqrt{\frac{-a}{3}cos(\frac{\phi}{3}+240º)}\)      and  cubic root   \(y_3= -\frac{p}{3} +x_3 \nonumber\)

    There is an alternate algebraic solution to solving a cubic equation, but it frequently breaks down close to 80-85% of critical pressure using MS Excel, due to software problems. Hence, the recommendation of the more stable trigonometric method, even though its math is seemingly odd. The trigonometric root-solving technique can break down somewhere near 95% of critical, using MS Excel because of numerical precision; but since the P-R EOS is only valid to 90% of critical, that is not a problem. Note that there are other cubic root solving algorithms, but these also require even higher levels of numerical precision and are often found to break down (e.g., division by zero or requiring math with imaginary numbers).

    When doing the above calculations, note that cos-1(ø) is only valid over the range between -1 and 1, returning values of ø between 1 to 0 radians, respectively. When doing automated calculations (via a spreadsheet like MS Excel), it is a good idea to check to see if the quantity \({-b/(2\sqrt{-a^3/27})} \nonumber\) is between -1 and 1, in order to avoid an error message. 

    In Equation 5-5 above, the cubic equation’s Zterm is unity, the Zis [- (1-B)] = B-1, the Z term is [A‑2B‑3B3] and the constant is [- (AB-B2-B3)]. Then the solution of that cubic equation, using Equation 5-6:

    • \(p=B-1 \nonumber\) 
    • \(q=A-2B-3B^2 \nonumber\) 
    • \(r= -(AB-B^2-B^3) \nonumber\)
    • \(a=(3q-p^2)/3 \nonumber\)
    • \(b=(2p^3 -9pq +27r)/27 \nonumber\)

    Intermediate parameters “p” and “r” will usually be negative, and parameter “q” will be positive; the quantity \(\frac{b^2}{4} + \frac{a^3}{27} \nonumber\) will always be negative: so there will be three unequal roots. 

    To use the trigonometric method, first solve for the angle ø (in radians) of Equation 5-6:

    \(\phi = cos^{-1}(-b/(2\sqrt{-a^3/27}) \nonumber\) and so the three roots, y1, y2, and ywill be:

    \(y_1= -p/3 +2\sqrt{-a/3} \times cos(\phi /3) \nonumber\), \(y_2= -p/3 +2\sqrt{-a/3} \times cos(\phi /3 + 2\pi/3) \nonumber\) and \(y_3= -p/3 +2\sqrt{-a/3} \times cos(\phi /3 + 4\pi/3) \nonumber\)

    The largest valued root is \(Z_v\), and the smallest valued root is ZL, with the intermediate valued root being discarded as meaningless. 

    Since the math in Part of both the EOS and the cubic solution might seem daunting, an example is provided with calculations.

    Example \(\PageIndex{1}\)

    Continuing the example of Part IV, of TCS vaporizing at 3.50 atm, and 347.05°K = 73.9°C, determine the values of Zv, ZL, and of ΔZ = Z- ZL. From the values of Zand Zdetermine the saturated vapor and liquid densities and then from the value of ΔZ determine the value of the latent heat of vaporization at 3.50 atm


    From Part III, Table 3-3, the pertinent property values for TCS are: MW = 135.452; Tc=479.15°K; Pc=41.15 atm ; and ω=0.2090. Tb , Vc and Zc are not needed for this example. 

    For use in the P-R EOS, at the examples temperature of 347.05°K, T r= 347.05/479.15 = 0.7243.

    Using Equation 5-5 above for the P-R EOS and plugging in the values for T, P, Tr, Tc, Pc, and ω:

    \(a=0.45724 \times 82.057^2 \times 479.15^2/41.15 = 1.7177E7 \nonumber \)

     \(b=0.07780 \times 82.057 \times 479.15/41.15 = 74.336 \nonumber \)

     \(\kappa=0.37464 + 1.54226 \times 0.02090 - 0.26992 \times 0.2090^2=.68518 \nonumber \)

    \(\alpha =(1+0.68516\times (1-0.7243^{0.5}))^2 =1.7575 @: 347.05°K \nonumber\)

    \(A=1.7575 \times 1.7177E7 \times 3.50/(82.057^2 \times 347.05^2)=0.13029 @: 347.05°K, 3.50 atm

    \(B=(74.336 \times 3.50)/(82.057 \times 437.05)= 9.136E-3 @: 347.05°K, 3.50 atm

    Then the cubic equation to solve (using Equation 5-6) is \(Z^3 +pZ^2 +qZ +r=0 \nonumber\), with terms:

    • \(p=B-1 \nonumber\) 
    • \(q=A-2B-3B^2 \nonumber\) 
    • \(r= -(AB-B^2-B^3) \nonumber\)
    • \(a=(3q-p^2)/3 \nonumber\)
    • \(b=(2p^3 -9pq +27r)/27 \nonumber\)
    • \(p=0.0091361-1 = -0.99086 \nonumber\)
    • \(q= 0.13029 -2\times 0.00913161 - 3\times 0.009136^2= 0.11177 \nonumber\)
    • \(r=-(0.13029 \times 0.009136 -0.009136^2 -0.009136^3)=-0.001106 \nonumber\)

    so for the trigonometric solution a and b are:

    \(a=(3 \times 0.11177 - (-0.99086)^2) /3 = -.21550 \nonumber\)

    \(b=(2(-0.99086)^3 -9\times (-0.99086)\times 0.11177 +27\times (-0.0001106))/27 =0.036254 \nonumber\)

    b2/4 +a3/27 = -4.2106E-5, so three real roots for Z are confirmed

    \(\phi = cos^{-1}(-b/(2\sqrt{-a^3/27}) \nonumber\)= 0.34376 radians and so the three roots are:

    •  \(Z_1= -p/3 +2\sqrt{-a/3} \times cos(\phi /3)= 0.86281\nonumber\)
    • \(Z_2= -p/3 +2\sqrt{-a/3} \times cos(\phi /3 + 2\pi/3) = 0.010947 \nonumber\) 
    • \(Z_3= -p/3 +2\sqrt{-a/3} \times cos(\phi /3 + 4\pi/3) = 0.11710 \nonumber\)

    Zis the largest root = 0.86281 and ZLis the smallest root = 0.01095 , with the Z3 root being discarded as meaningless. Zand ZL will be used in Part VI to give estimates of fugacity coefficients, which partially account for departures from ideal mixture behavior.

    ΔZ = Z- ZL= 0.86281-0.01095 = 0.85186 (note that this is less than the unity value of the Ideal Gas Law of Equation 5-1, and the unity assumption of Part II, Equation 2-1 which led to simplified VP equations).

    If Equation 5-1 is re-written as \(\frac{1}{V} = \frac{P}{ZRT} \nonumber\) that gives the density in molar units,

    and then multiply by molecular weight to get density in more familiar mass units.

    Predicted vapor density @ 3.50 atm, \(73.9°C = 3.50/(0.86281 \times 82.057 \times 347.05) \nonumber\) = 1.4244E-4 g-mole/cc . In mass units = (135.45 x 1.4244E-4) = 1.929E-2 g/cc

    Liquid density @3.50 atm,\(73.9°C = 3.50/(0.010947 \times 82.057 \times 347.05) \nonumber\) = 1.1227E-2 g-mole/cc . In mass units = (135.45 x 1.1227E-2) = 1.521 g/cc

    Looking at some TCS saturated liquid phase density data, it can be seen that the predicted saturated liquid density by the P-R EOS is high by about 10%, with a somewhat greater inaccuracy for sub-cooled liquid TCS. So unless there is better data available, the P-R EOS should not be thought of as better than ±15% on saturated phase density. Comparable TCS saturated vapor phase density data does not seem to exist, so no accuracy estimates can be made.

    In order to estimate the latent heat of vaporization at these conditions, it is necessary to go back to Part IV, Equation 4.5 and the example values: Tr= 0.7243; A=9.9796; B1=0.38395; B2= -0.11827; B3=0.063708; c=2.5687; n=5.7950; and k=0.0446.

    \[\psi=\frac{\Delta H_{\nu}}{\Delta Z RT} = A(1-B_{1}T_{r} +B_{2}T_{r}^2-B_{3}T_{r}^3 +c(T_{r}^n -k)) \nonumber\]

    and substituting in the values from above, but using a value for R=1.9859 to get ΔHv in energy units of calories/gram-mole;

    \[\Delta H_{\nu}/(0.85186 \times 1.9859 \times 347.05)=9.9796 \times (1-0.38395\times 0.7243 -0.11827 \times 0.7243^2 -0.063708 \times 0.7243^3) + 2.5687(0.7243{5.7950}-0.0446) \nonumber\]

    So ΔHv is predicted at = 3,890 cal/g-mole. Multiplying by the MW, the latent heat of vaporization is predicted as 28.7 cal/gram at 3.50 atm and 72.9°C. Better estimates of the TCS latent heat at these conditions by previous work done under NASA contract showed a value closer to 50 cal/gram ( but these were estimates only, not direct measurements).

    This illustrates the potential power of using an EOS with a thermodynamically consistent VP equation, but also shows the degree of inaccuracy of pairing the P-R EOS with the VP derivative of Part IVEquation 4.5. In this case, it is unknown how much of the inaccuracy is caused by the EOS, and how much caused by VP equation’s derivative. 

    As will be seen in Part VI, the best use of the P-R EOS is just to estimate fugacity coefficients, which are relative quantities rather than absolute.