Information • Textbooks • Media • Resources edited by
Computer Bulletin Board
Steven D. Gammon University of Idaho Moscow, ID 83844
Virial Coefficients Using Different Equations of State Charles B. Wakefield* and Constance Phillips Department of Science, Mathematics and Computer Science, University of Texas of the Permian Basin, Odessa, TX 79762; *
[email protected] Inexpensive desktop computers and powerful computer algebra software packages have made it possible for students and faculty to look at problems unthinkable a few years past because of the intensive algebra involved. One such problem is the derivation of expressions for the first three virial coefficients when using different equations of state. Using Mathematica, expressions for the virial coefficients at constant volume, constant pressure, and the fugacity have been derived for several equations of state. Our purpose is to illustrate the actual Mathematica instructions that lead to the virial coefficients for the well-known van der Waals equation of state. The set of Mathematica instructions for this equation of state will serve as a model for other more sophisticated equations of state. We selected the Redlich–Kwong and the Gibbons–Laughton modification to the Redlich–Kwong equations of state. We will use the abbreviations “vdw” for the van der Waals equation (1), “RK” for the Redlich–Kwong equation (2), and “GLRK” for Gibbons–Laughton modification to the Redlich– Kwong equation (3) as suffixes for the different equations of state. The virial coefficients will be designated as usual with the above suffixes to differentiate between the equations of state. For example, BVvdw represents the second virial coefficient at constant volume for the van der Waals equation of state. CpRK is the third virial constant at constant pressure for the Redlich–Kwong equation of state. The constants for the different equations of state are also distinguished by the appropriate suffix. Therefore, the “a” in the van der Waals equation is written as “avdw”. We performed the derivation with the van der Waals equation of state. We thought that most physical chemistry students would be familiar with this equation and could therefore pay more attention to the instructions to the computer algebra program. Also, the output from the computer algebra system would look familiar, thus giving a sense of confidence in the program and the method. Development The well-known van der Waals equation for one mole is written as follows: P = RT – a 2 (1) V–b V To transform the equation into a form compatible with the virial equation expressed in terms of the volume, we must multiply by V/RT. This gives the following:
Zvdw = PV = V – a RT V – b RTV
(2)
This compares to the temperature dependent virial equation in terms of volume, which is:
Z = PV = 1 + RT
BV T V
+
CV T V
2
+
DV T V3
+…
(3)
The instructions to Mathematica are indicated in bold type (4). A semicolon at the end of a Mathematica instruction will suppress the output for that instruction. This was desirable for some instructions because the algebraic results were so long. The output will follow the bold typed instruction unless a semicolon was used to suppress the output. What appears in the bold type is exactly the instruction given to and then executed by Mathematica. The following is the instruction set for determining the virial coefficients and the fugacity for the van der Waals equation that will serve as a model for our other two equations of state. Zvdw[v_ ] = v/(v - bvdw) - avdw(RTv) V avdw RT v + bvdw + V seriesvdw = Normal[Series[Zvdw[v], {v, Infinity, 4}]] bvdw – avdw bvdw4 bvdw3 bvdw2 RT 1+ + + + 4 3 2 v v v v BVvdw[T_] := Coefficient[seriesvdw, 1/v] CVvdw[T_] := Coefficient[seriesvdw, 1/v^2] DVvdw[T_] := Coefficient[seriesvdw, 1/v^3] Print[“BVvdw(T) = “, BVvdw[T]] Print[“CVvdw(T) = “, CVvdw[T]] Print[“DVvdw(T) = “, DVvdw[T]]
BVvdw(T) = bvdw – avdw RT CVvdw(T) = bvdw2 DVvdw(T) = bvdw3
The BVvdw result is found in most physical chemistry textbooks (5); however, the terms for CVvdw and DVvdw have little meaning because they are not functions of the temperature but are constants for the van der Waals equation. They would appear as lines, which cut through the center of the virial data with no curvature, giving a kind of average value. This is not true for the more sophisticated equations of state. Next, we will derive the pressure-dependent coefficients from the volume-dependent coefficients. To do this, we solve the volume-dependent virial equation for the pressure P. Then, after substituting the results into the pressure-dependent virial
JChemEd.chem.wisc.edu • Vol. 77 No. 10 October 2000 • Journal of Chemical Education
1371
Information • Textbooks • Media • Resources
equations, we obtain the coefficients of the volume terms, all of which is done in the following segment of Mathematica code. This approach gives us the expressions of the coefficients in the virial equation for pressure in terms of the volume equation coefficients. While we carry the calculation only up to the Dp coefficient, the process could easily be extended to higher terms if needed. Though the terms above BVvdw are not useful, we include them here because the same method could be used for more sophisticated equations of state in the same manner as illustrated below. For several of the instructions below the output is excessive, it is suppressed with a semicolon at the end of the Mathematica instruction. The constant-pressure virial coefficient Bvdw should be recognizable to the physical chemical student. It is included in most physical chemistry texts. The other coefficients are much more algebra intensive but are found with ease using Mathematica.
RT ln
ln
1 + Bpp +
+
ln
sol1 = Solve[Coefficient[fun, 1/v] = = Bv, Bp]//Flatten sol2 = Solve[Coefficient[fun, 1/v^2] = = Cv, Cp]/.sol1//Flatten sol3 = Solve[Coefficient[fun, 1/v^3] = = Dv, Dp]/.sol1/.sol2//Flatten Bv Bp → RT
Cp →
Bv2 – Cv R2 T2
Dp →
2Bv Bv2 – Cv + BvCv – Dv R3 T3
Bpvdw[T_ ] = BVvdw[T]/(RT) Cpvdw[T_ ] = -(BVvdw[T]^2 - CVvdw[T])/(RT)^2//Expand; Dpvdw[T_ ] = -(-2BVvdw[T](BVvdw[T]^2 - CVvdw[T]) + BVvdw[T]CVvdw[T] - DVvdw[T]/(RT)^3//Expand; Print[“Bpvdw(T) = “, Bpvdw[T]] Print[“Cpvdw(T) = “, Cpvdw[T]] Print[“Dpvdw(T) = “, Dpvdw[T]]
Bpvdw(T) =
bvdw – avdw RT RT
Cpvdw(T) =
avdw2 bvdw + 2 avdw R4 T4 R3 T3
2 avdw3 6 avdw2 bvdw – 3 avdw bvdw2 Dpvdw(T) = – + R5 T5 R6 T6 R4 T4
V – RT dP P
f = P
P= P
P= 0
V – 1 dP RT P
f = P
P= P
(5)
P= 0
1 + B + C P + D P 2 + … – 1 dP (6) P P P P P
ln
f 2 3 = B P P + 1 CP P + 1 DP P + … 2 3 P
(7)
If we convert this to the exponential form we have f = P exp B P + 1 CP P 2 + 1 DP P 3 2 3
(8)
Now we can express f for the van der Waals equation as follows: fvdw = p Exp Bpvdw T p +
Cpvdw T 2 Dpvdw T 3 p + p 2 3
2 avdw3 6 avdw2 bvdw – 3 avdw bvdw2 E^ 1 p3 – + + 3 R5 T5 R6 T6 R4 T4 2 1 p2 – avdw + 2 avdw bvdw 4 4 2 R3 T3 R T
+
p bvdw – avdw RT RT
p
This same derivation process using Mathematica can be used for other equations of state. We will now apply the above model to the other two equations of state discussed in this paper. For these cases we will show only the pertinent output from the Mathematica instructions. Results for Other Equations of State In this section we will give the results for the other equations of state for the virial coefficients and the fugacity. These results were obtained using the same approach as for the van der Waals equation. The expression for DP for GLRK is too long to include here, but it can be derived with no difficulty using the method outlined above.
Redlich–Kwong (RK) The RK equation (2) is
Fugacity (vdw)
a P RK = RT – V–b V–b V T
Now we illustrate how to arrive at an expression for the fugacity assuming the gas behaves according to the van der Waals equation. In most physical chemistry texts we have the following relation for the fugacity:
R 2T c2.5
1372
(4)
Equation 6 is easily integrated to give
Dpp3
fun = Zp/.psub//Expand;
P= 0
Now we solve the virial equation in terms of pressure for V/RT and substitute into eq 5 to get
Zp = 1 + Bpp + Cpp^2 + Dpp^3 Cpp2
P= P
If we divide eq 4 by RT we have
psub = {p -> RT/v + BvRT/v^2 + CvRT/v^3 + DvRT/v^4}
p → DvRT + CvRT + BvRT + RT v v4 v3 v2
f = P
where a = 0.472
Pc
and b = 0.867
RT c Pc
Journal of Chemical Education • Vol. 77 No. 10 October 2000 • JChemEd.chem.wisc.edu
(9)
.
Information • Textbooks • Media • Resources
Figure 1. This plot, generated using Mathematica, shows that the results obtained by using the Mathematica-derived BvRK coefficients give reasonable results when compared to the data for (◊) H2, (䉭) CO2, and (ⵧ) CH4.
All the expressions that follow are determined in the same manner as was illustrated for the van der Waals example. First we give the virial form. bRK4 + 1+
aRK bRK2 aRK bRK3 bRK3 – 3/2 RT RT3/2 + + 4 3 v v aRK bRK bRK2 + bRK – aRK RT3/2 RT3/2 + v v2
Next we show the results for the virial coefficients at constant volume. BVRK(T) = bRK – aRK RT3/2 bRK CVRK(T) = bRK2 + aRK 3/2 RT
where a =
9KP c
, b=
KRTc 3P c
, and K = 2 0.3 – 1.
The primary difference between these two equations is the function of temperature. The Redlich–Kwong simply uses the inverse of the square root of the temperature, whereas the Gibbons–Laughton uses the temperature as a function of their α expression, where α = X(TR – 1) + Y(TR0.5 – 1). For this expression TR = T/Tc and Tc is the critical temperature of the gas. The values of the constants X and Y are considered to be “substance-dependent parameters” with Soave’s “m” factor, which is based on the Pitzer acentric factor ω2. First we give the virial form.
aRK bRK RT3/2
The following are the coefficients at constant pressure. BpRK(T) =
2 R Tc 2
T X + – 1.+ T aGLRK 1.+ – 1.+ Tc Tc bGLRK – RT 1 + v
2
DVRK(T) = bRK3 –
Figure 2. This plot shows that the results for BvGLRK coefficients also give reasonable results when compared to the data for (◊) H2, (䉭) CO2, and (ⵧ) CH4.
bGLRK2 +
bRK – aRK RT3/2
RT
aRK2 bRK CpRK(T) = – 4 5 + 3 aRK R T R3 T7/2
T X + – 1.+ T aGLRK bGLRK 1.+ – 1.+ Tc Tc RT v2
Y +
0.5
Y +
T X + – 1.+ T aGLRK bGLRK2 1.+ – 1.+ Tc Tc RT v3
0.5
T X + – 1.+ T aGLRK bGLRK3 1.+ – 1.+ Tc Tc bGLRK + RT v4
0.5
bGLRK3 –
7 aRK bRK2 9 aRK2 bRK 2 aRK3 – DpRK(T) = – + 5 6 15/2 R T R4T9/2 R 6T
0.5
Y +
Y
4
Finally, we have the expression for the fugacity for the Redlich–Kwong equation of state. 2 2 3 1 p3 – 2 aRK + 9 aRK bRK – 7 aRK bRK + 5 6 4 9/2 3 6 15/2 T R R T R T p bRK – aRK aRK2 RT3/2 1 3 aRK bRK 2 – E^ p + + RT 2 R4T5 R3 T7/2
Gibbons–Laughton (GLRK) The GLRK (3) is as follows: P GLRK = RT – a αT V–b V V+b
p
Shown next are the expressions for the virial coefficients at constant volume.
X – BVGLRK = bGLRK – 1.aGLRK + 1.aGLRK RT RT
(10)
aGLRK X + 1.aGLRK Y – RT RTc
T aGLRK Tc RT
JChemEd.chem.wisc.edu • Vol. 77 No. 10 October 2000 • Journal of Chemical Education
0.5
Y
1373
Information • Textbooks • Media • Resources
Conclusion
CVGLRK = bGLRK2 + 1.aGLRK bGLRK – 1.aGLRK bGLRK X + RT RT aGLRK bGLRK T Tc aGLRK bGLRK X – 1.aGLRK bGLRK Y + RT RT RTc
DVGLRK = bGLRK3 –
0.5
Y
1.aGLRK bGLRK2 + 1.aGLRK bGLRK2 X – RT RT
2 T aGLRK bGLRK2 X 1.aGLRK bGLRK2 Y aGLRK bGLRK Tc + – RTc RT RT
0.5
Y
We show here only the expression for BpGLRK—the expressions for CpGLRK and DpGLRK are excessively long. BpGLRK =
bGLRK– 1.aGLRK + 1.aGLRK X – aGLRK X + 1.aGLRK Y – RT RT RT RTc RT
aGLRK T Tc RT
0.5
Y
The expression for the fugacity for GLRK is also omitted because of its length. However, having expressions for the pressure coefficients we can find an expression for the fugacity as was carried out for the van der Waals equation. The instruction sent to Mathematica needed to retrieve the expression for the fugacity for GLRK is as follows: DpGLRK[T] p + p ] [BpGLRK[T]p+ CpGLRK[T] 2 3
fGLRK = p Exp
2
3
Expression Confirmation Graphics We have included two graphics (Figs. 1 and 2) to show that the expressions derived for the Redlich–Kwong and the Gibbons–Laughton equations of state return some reasonable results. The fit of the Cv data was about as good as the results shown for Bv (6 ).
1374
We see that Mathematica is a powerful tool that can be used to do the extensive algebra inherent in these types of derivations and calculations. As the speed and sophistication of computers and the computer algebra systems increase and their cost decreases, it will become increasingly important to teach our students how to utilize these powerful tools. With the click of a mouse we can retrieve algebraic expressions that could never be derived by hand. We can make calculations for different values of an independent variable using very long complicated expressions, and with the help of the plotting routines, we can instantly view the effect of making changes in an equation. These activities have in the past been reserved for the large computers that were out of the reach of most undergraduate students. This is no longer the case. Many students have their own computers with Mathematica. Our work also supports the claim that the two-parameter Redlich–Kwong equation fits the virial coefficient data better than the van der Waals equation (7). Literature Cited 1. Atkins, P. Physical Chemistry, 5th ed.; Freeman: New York, 1994; Chapter 1. 2. Shah, K. K.; Thodos, G. Ind. Eng. Chem. 1956, 57 (3), 30– 37. 3. Gibbons, R. M.; Laughton, A. P. J. Chem. Soc, Faraday Trans. 2 1984, 80, 10119–10138. 4. Wolfram, S. The Mathematica Book, 4th ed.; Cambridge University Press: New York, 1999. 5. Barrow, G. M. Physical Chemistry, 5th ed.; McGraw-Hill: New York, 1988; Chapter 2. 6. Dymond, J. H.; Smith, E. B. The Virial Coefficients of Pure Gases and Mixtures; Clarendon Press: Oxford, 1980. 7. Andrews, F. C. Thermodynamics; Wiley: New York, 1971; p 118.
Journal of Chemical Education • Vol. 77 No. 10 October 2000 • JChemEd.chem.wisc.edu
