5: Equation of State
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 \label{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), Z L is very close to zero for either a sub-cooled or saturated liquid, and Z v 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 E quation 2-2 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 Z v >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 E quation ( \ref{5-1} ) is re-written as
\[Z=\dfrac{PV}{RT} \label{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: \(Z_V\) and \(Z_L\). 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 \label{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 Z c . There are very few fluids that have such a high Z c - 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 Z c of 0.307 at the critical point. While this Z c value is still higher than that of chlorosilanes, the accuracy of Z v and Z L values are adequate for fugacity estimation , through 90% of critical pressure (P c ). Since commercial applications only rarely approach 90% of P c 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} \label{5-4} \]
with
- \(a=\dfrac{0.45724R^2T_{c}^2}{p_{c}} \nonumber\)
- \(b=\dfrac{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\) when ω < or = 0.49
- \(\kappa=0.379642+1.48503\omega-0.164423\omega^2 +0.016666\omega^3 \nonumber\) when ω>0.49
- \(T_{r}=\dfrac{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=\dfrac{\alpha a p}{R^2T^2} \nonumber\) and \(B=\dfrac{bp}{RT} \nonumber\)
\[Z^3 - (1-B)Z^2+(A-2B-3B^2)Z -(AB-B^2-B^3) = 0 \label{5-5} \]
At a given temperature, \(T\), and for fluid physical properties T c , P c , 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 E quation ( \ref{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 Z v ; the smallest root is the value of \(Z_L\); and the third root is discarded as having no physical meaning.
Then knowing \(Z_v\) and \(Z_L\), 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 Z v and Z L , 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) . While the values for saturated vapor-phase densities appear to be quite reasonable up to 80% of critical pressure, the estimated values for saturated liquid-phase densities are generally found to be rather poor. To get an idea of this inaccuracy, compare the ambient temperature liquid density value to the EOS value for saturation pressure equivalent to ambient temperature. For a good work-around to this problem, see the attachment "Saturated Liquid Densities.pdf". For a good work-around to the inaccuracy of duplicating the latent heat of vaporization at near-ambient pressures, see the attachment "Heat of Vaporization.pdf".
A cautionary statement is in order regarding the use of an EOS (whether Van der Waals, SRK, P-R or others) with liquid mixtures. In some hydrocarbon applications, VLE convergence calculations attempt to combine the critical properties of a mixture's components, as a "psuedo-critical" mixture property, and use these "psuedo-criticals" in EOS calculations along with various mixing models. While there has been some success in this approach for certain hydrocarbons ( mostly in the oil and gas industry), such approach rarely gives good results with inorganic polar or mildly polar mixtures, such as chlorosilanes. Various papers have been issued over the recent years promoting such approach, suggesting the use of "binary interaction parameters" to estimate an EOS for a mixture. There are simply too many differences in polarity, acentric factor, and individual fluid critical properties for this to work. Therefore, the reader is discouraged from such practice. Sometimes calculations must be worked out, checked along the way, and results validated - rather than taking an "easy way out". Certainly with computer automation, there is no reason not to do the full calculation procedure.
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 y 1 , y 2 and y 3
\[y^3 +py^2 + qy + r =0 \label{5-6} \]
may be transformed to the format , \(x^3+ax+b =0 \) by substituting for \(y\) the value, \(x - \frac{p}{3} \nonumber\). Where \(a= \frac{1}{3}(3q-p^2) \) and \(b=\frac{1}{27}(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.
For saturated liquids and vapors, the last case is always true so 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 x 1 , x 2 , and x 3 will have the following values:
- \(x_1= 2 \sqrt{\frac{-a}{3}} \times cos(\frac{\phi}{3}) \) and so cubic root \(y_1= -\frac{p}{3} +x_1 \nonumber\)
- \(x_2 =2 \sqrt{\frac{-a}{3}} \times cos(\frac{\phi}{3}+120º) \) and cubic root \(y_2= -\frac{p}{3} +x_2 \nonumber\)
- \(x_3=2 \sqrt{\frac{-a}{3}} \times 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 (which is frequently needed) to evaluate Z for superheated vapors. But for liquids and vapors near saturation, this alternate solution 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 still 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 E quation ( \ref{5-5} ) above, the cubic equation’s \(Z^3\) term is unity, the \(Z^2\) is \([- (1-B)] = B-1\), the \(Z\) term is \([A‑2B‑3B^3]\) and the constant is \([- (AB-B^2-B^3)]\). Then the solution of that cubic equation, using E quation ( \ref{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; and the quantity \(\frac{b^2}{4} + \frac{a^3}{27} \nonumber\) will always be negative for liquids and vapors at or near saturation: so there will be three unequal roots. Note that for superheated vapors or vapors above critical pressure (aka gases), the quantity \(\frac{b^2}{4} + \frac{a^3}{27} \) can be >0, in which case there is only one real root for Z, and an alternate algebraic solution method will be needed ( i.e., the trigonometric method breaks down) . In such case, \(Z\) will be >1.
To use the trigonometric method to get \(Z_v\) and \(Z_L\), first solve for the angle ø (in radians) of E quation ( \ref{5-6} ) :
\(\phi = cos^{-1}(-b/(2\sqrt{-a^3/27}) \nonumber\) and so the three roots, y 1 , y 2 , and y 3 will 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 Z L , with the intermediate valued root being discarded as meaningless.
Since the math in Part V 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 Z v , Z L , and of ΔZ = Z v - Z L . From the values of Z v and Z L determine 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 .
Solution
From Part III , Table 3-3, the pertinent property values for TCS are: MW = 135.452; T c =479.15°K; P c =41.15 atm ; and ω=0.2090. T b , V c and Z c 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 E quation ( \ref{5-5} ) above for the P-R EOS and plugging in the values for T, P, T r , T c , P c , 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.68518\times (1-0.7243^{0.5}))^2 =1.2145 @: 347.05°K \nonumber\)
\(A=1.2145 \times 1.7177E7 \times 3.50/(82.057^2 \times 347.05^2) \nonumber \)=0.09003 @: 347.05°K, 3.50 atm
\(B=(74.336 \times 3.50)/(82.057 \times 437.05)= 9.136E-3 \nonumber \) @: 347.05°K, 3.50 atm
Then the cubic equation to solve (using E quation ( \ref{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.09003 -2\times 0.00913161 - 3\times 0.009136^2= 0.071516 \nonumber\)
- \(r=-(0.09003 \times 0.009136 -0.009136^2 -0.009136^3)=-0.00073819 \nonumber\)
so for the trigonometric solution a and b are:
\(a=(3 \times 0.071516 - (-0.99086)^2) /3 = -.25575 \nonumber\)
\(b=(2(-0.99086)^3 -9\times (-0.99086)\times 0.071516 +27\times (-0.00073819))/27 =-0.049179 \nonumber\)
b 2 /4 +a 3 /27 = -1.49326E-5, so three real roots for Z are confirmed
\(\phi = cos^{-1}(-b/(2\sqrt{-a^3/27}) \nonumber\)= 0.15588 radians and so the three roots are:
- \(Z_1= -p/3 +2\sqrt{-a/3} \times cos(\phi /3)= 0.91345\nonumber\)
- \(Z_2= -p/3 +2\sqrt{-a/3} \times cos(\phi /3 + 2\pi/3) = 0.012439 \nonumber\)
- \(Z_3= -p/3 +2\sqrt{-a/3} \times cos(\phi /3 + 4\pi/3) = 0.064968 \nonumber\)
Z v is the largest root = 0.91345 and Z L is the smallest root = 0.01244 , with the Z 3 root being discarded as meaningless. Z v and Z L will be used in Part VI to give estimates of fugacity coefficients, which partially account for departures from ideal mixture behavior.
ΔZ = Z v - Z L = 0.91345-0.01244 = 0.90101 (note that this is less than the unity value of Part II , Equation 2-3 and 2-4 , which led to simplified VP equations).
If E quation ( \ref{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.91345 \times 82.057 \times 347.05) \nonumber\) = 1.3455E-4 g-mole/cc . In mass units = (135.45 x 1.3455E-4) = 1.822E-2 g/cc
Predicted liquid density @3.50 atm ,\(73.9°C = 3.50/(0.012439 \times 82.057 \times 347.05) \nonumber\) = 9.8804E-3 g-mole/cc . In mass units = (135.45 x 9.8804E-3) = 1.338 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 8%, 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 ±10% on saturated phase density. Comparable TCS saturated vapor phase density data does not seem to exist, so no accuracy estimates can be made. See also the attached file "Saturated Liquid Densities", which can alleviate the EOS' inaccuracy on liquid densities.
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: T r = 0.7243; A=9.9796; B 1 =0.38395; B 2 = -0.11827; B 3 =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 ΔH v in energy units of calories/gram-mole;
\[\Delta H_{\nu}/(0.90101 \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 ΔH v is predicted at = 4,114 cal/g-mole. Multiplying by the MW, the latent heat of vaporization is predicted as 30.4 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 indicated a value closer to 50 cal/gram (but these were still estimates, not direct measurements, and are possibly in error as well).
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 IV , Equation 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.