Physics:Virial expansion

Short description: Series expansion of the equation of state for a many-particle system

The virial expansion is a model of thermodynamic equations of state. This model expresses the pressure P of a gas in local equilibrium as a power series of the density. This equation may be represented in terms of the compressibility factor, [math]\displaystyle{ Z }[/math], as [math]\displaystyle{ Z \equiv \frac{P}{RT\rho} = A + B\rho + C\rho^2 + \cdots }[/math] This equation was first proposed by Kamerlingh Onnes. The terms A, B, and C represent the virial coefficients. The leading coefficient, A, is defined as the constant value of 1, which enforces that the equation reduces to the ideal gas expression as the gas density approaches zero.

Second and third virial coefficients

The second, B, and third, C, virial coefficients have been studied extensively and tabulated for many fluids for more than a century. Two of the most extensive compilations are in the books by Dymond. and NIST Thermo Data Engine Database and its Web Thermo Tables. Tables of second and third virial coefficients of many fluids are included in these compilations.

The 2nd and 3rd virial coefficients of argon

The second and third virial coefficients as functions of temperature are shown in the following figure for argon. Reduced temperature and reduced virial coefficients, scaled by respective critical properties, are all dimensionless. Notice in the figure that the second virial coefficient decreases monotonically as temperature is lowered. However, the third virial coefficient has a bell shape. It increases as temperature is lowered to the critical temperature, then it passes through a peak and decreases rapidly to zero as temperature is lowered from the critical point to the triple point. It is physically unreasonable for it to decrease below the critical temperature. This is because the third virial coefficient theoretically represents the repulsive force among molecules, which is expected to increase at a lower temperature, as molecules are pressed together. The behaviors just described are similar to many other gases, as shown in the figure below.

The reason why the third virial coefficient decreases below the critical temperature is as follows. Generally, the PρT isotherms are measured in the gaseous phase. Below the critical temperature, the gaseous phase condenses and coexists with the liquid phase, and the PρT isotherm becomes flat. Saturation pressure does not change until the gas condenses completely to liquid, and then pressure rises as density increases. There is a large gap between pure gaseous phase and pure liquid phase where no useful pressure data, except saturation pressure, are available. If only data in the gaseous phase were analyzed, the third virial coefficient becomes very small, because the PρT isotherm is almost linear in the gaseous phase. However, if data points in the pure liquid phase are included, a second-order regression would give a large third virial coefficient. The third virial coefficient thus derived would increase monotonically as the temperature is lowered from the critical point to the triple point.

The expectation that the third virial coefficient is a monotonically increasing function of 1/T can be verified with equations of state which accurately predicted the PρT isotherms in the saturation region where gaseous and liquid phases coexist.

Casting equations of the state into virial form

Most equations of state can be reformulated and cast in virial equations to evaluate and compare their implicit second and third virial coefficients. The seminal van der Waals equation of state was proposed in 1873: [math]\displaystyle{ P = \frac{RT}{\left(v-b\right)} - \frac{a}{v^2} }[/math] where v = 1/ρ is molar volume. It can be rearranged by expanding 1/(v − b) into a Taylor series: [math]\displaystyle{ Z = 1 + \left(b-\frac{a}{RT}\right)\rho + b^2\rho^2 + b^3\rho^3 + \cdots }[/math]

In the van der Waals equation, the second virial coefficient has roughly the correct behavior, as it decreases monotonically when the temperature is lowered. The third and higher virial coefficients are independent of temperature, and are not correct, especially at low temperatures.

Almost all subsequent equations of state derived from the van der Waals equation, like those from Dieterici, Berthelot, Redlich-Kwong, Peng-Robinson, etc., suffer from the singularity introduced by 1/(v - b). Hence they cannot represent accurately the PρT isotherms at temperatures below critical temperature. Many of them produce adequate second virial coefficients, but most give incorrect third virial coefficients .

Other equations of state, started by Beattie-Bridgeman, however, are more closely related to virial equations, and show to be more accurate in representing behavior of fluids in both gaseous and liquid phases. They can be reformulated into virial equations of state, and compared with one another. The Beattie-Bridgeman equation of state, proposed in 1928, [math]\displaystyle{ p=\frac{RT}{v^2}\left(1-\frac{c}{vT^3}\right)(v+B)-\frac{A}{v^2} }[/math] where

• [math]\displaystyle{ A = A_0 \left(1 - \frac{a}{v} \right) }[/math]
• [math]\displaystyle{ B = B_0 \left(1 - \frac{b}{v} \right) }[/math]

can be rearranged: [math]\displaystyle{ Z=1 + \left(B_0 -\frac{A_0}{RT} - \frac{c}{T^3}\right) \rho - \left(B_0 b-\frac{A_0 a}{RT} + \frac{B_0 c}{T^3}\right) \rho^2 + \left(\frac{B_0 b c}{T^3}\right) \rho^3 }[/math] This equation of state represents well the second virial coefficient. However, the third virial coefficient has the wrong sign. Thus, it fails to represent isotherms close to and below the critical temperature.

The Benedict-Webb-Rubin equation of state of 1940 represents better isotherms below the critical temperature: [math]\displaystyle{ Z = 1 + \left(B_0 -\frac{A_0}{RT} - \frac{C_0}{RT^3}\right) \rho + \left(b-\frac{a}{RT}\right) \rho^2 + \left(\frac{\alpha a}{RT}\right) \rho^5 + \frac{c\rho^2}{RT^3}\left(1 + \gamma\rho^2\right)\exp\left(-\gamma\rho^2\right) }[/math]

More improvements are achieved by Starling in 1972: [math]\displaystyle{ Z = 1 + \left(B_0 -\frac{A_0}{RT} - \frac{C_0}{RT^3} + \frac{D_0}{RT^4} - \frac{E_0}{RT^5}\right) \rho + \left(b-\frac{a}{RT}-\frac{d}{RT^2}\right) \rho^2 + \alpha\left(\frac{a}{RT}+\frac{d}{RT^2}\right) \rho^5 + \frac{c\rho^2}{RT^3}\left(1 + \gamma\rho^2\right)\exp\left(-\gamma\rho^2\right) }[/math]

Following are plots of reduced second and third virial coefficients against reduced temperature according to Starling:

The 2nd and 3rd virial coefficients for 12 fluids

The exponential terms in the last two equations correct the third virial coefficient so that the isotherms in the liquid phase can be represented correctly. The exponential term converges rapidly as ρ increases, and if we took only the first two terms in its Taylor expansion series, [math]\displaystyle{ 1-\gamma\rho^2 }[/math], and multiply it with [math]\displaystyle{ 1 + \gamma\rho^2 }[/math], the result is [math]\displaystyle{ 1 - \gamma^2\rho^4 }[/math]. It thus contributes a [math]\displaystyle{ c / RT^3 }[/math] term to the third virial coefficient, and one term to the eighth virial coefficient, which can be ignored.

After the expansion of the exponential terms, the Benedict-Webb-Rubin and Starling equations of state have this form: [math]\displaystyle{ Z = 1 + b\rho_r + c\rho_r^2 + f\rho_r^5 }[/math] The fourth and fifth virial coefficients are zero. After the third virial term, the next significant term is the sixth virial coefficient. It seems that the first three virial terms dominate the compressibility factor for fluids, down to [math]\displaystyle{ 0.5T_c }[/math], and up to [math]\displaystyle{ 2.0\rho_c }[/math].

Re-analysis of the data by Starling

Re-analyzing the data reported by Starling, the virial coefficients are best represented by [math]\displaystyle{ b=b_0+\frac{b_1}{t_r}+\frac{b_2}{t_r^2}+\frac{b_3}{t_r^3} }[/math] [math]\displaystyle{ c=c_0+\frac{c_1}{t_r}+\frac{c_2}{t_r^2}+\frac{c_3}{t_r^3} }[/math] [math]\displaystyle{ f=f_0+\frac{f_1}{t_r} }[/math] b and c could be determined using simple second order regression analysis from experimental PρT isotherms. [math]\displaystyle{ b_0 - b_3 }[/math] and [math]\displaystyle{ c_0 - c_3 }[/math] could then be determined using third order regression analysis on b and c. [math]\displaystyle{ f_0 - f_1 }[/math] could then be determined by analyzing residues in compressibility factor after the first three virial terms are removed from the virial equation. The data by Starling are re-analyzed and the results are shown in the following table. These coefficients are all dimensionless, they are all scaled with critical molar volumes and critical temperature.

Fluid	[math]\displaystyle{ b_0 }[/math]	[math]\displaystyle{ b_1 }[/math]	[math]\displaystyle{ b_2 }[/math]	[math]\displaystyle{ b_3 }[/math]	[math]\displaystyle{ c_0 }[/math]	[math]\displaystyle{ c_1 }[/math]	[math]\displaystyle{ c_2 }[/math]	[math]\displaystyle{ c_3 }[/math]	[math]\displaystyle{ f_0 }[/math]	[math]\displaystyle{ f_1 }[/math]
Methane	0.440	−1.171	−0.236	−0.210	0.364	−0.275	−0.014	0.396	0.0319	1.71×10−03
Ethane	0.330	−0.806	−0.363	−0.378	0.553	−0.675	−0.038	0.680	0.0461	2.63×10−03
Propane	0.288	−0.706	−0.245	−0.575	0.532	−0.546	−0.308	0.843	0.0334	1.89×10−02
n-Butane	0.377	−0.916	−0.115	−0.610	0.547	−0.519	−0.347	0.871	0.0305	2.04×10−02
i-Butane	0.438	−1.051	−0.172	−0.401	0.483	−0.342	−0.021	0.538	0.0194	1.19×10−03
n-Pentane	0.481	−1.056	−0.166	−0.560	0.668	−0.720	−0.204	0.841	0.0411	1.17×10−02
i-Pentane	0.242	−0.674	−0.306	−0.520	0.815	−0.943	−0.194	0.868	0.0484	9.99×10−03
n-Hexane	0.435	−0.636	−0.358	−0.759	0.848	−1.275	−0.105	1.120	0.0604	4.98×10−03
n-Heptane	0.493	−0.798	−0.636	−0.428	0.589	−0.738	−0.017	0.814	0.0508	1.21×10−03
n-Octane	0.600	−0.744	−0.456	−0.763	0.174	−0.197	−0.272	0.919	0.0144	1.99×10−02
Nitrogen	0.502	−1.380	0.092	−0.333	0.400	−0.276	−0.027	0.322	0.0279	2.72×10−03
CO2	0.178	−0.044	−1.517	0.039	0.428	−0.422	−0.008	0.687	0.0490	9.52×10−04
H2S	0.191	−0.927	−0.078	−0.366	1.093	−1.227	−0.001	0.577	0.0578	8.37×10−05

Cubic virial equation of state

The three-term virial equation or a cubic virial equation of state [math]\displaystyle{ Z = 1+B\rho+C\rho^2 }[/math] has the simplicity of the Van der Waals equation of state without its singularity at v = b. Theoretically, the second virial coefficient represents bimolecular attraction forces, and the third virial term represents the repulsive forces among three molecules in close contact. Intuitively, we should expect B to become negative at low temperature, while C would remain positive to counterbalance the effect of B and pushes Z and hence the pressure to high values as ρ increases.

With this cubic virial equation, the coefficients B and C can be solved in closed form. Imposing the critical conditions: [math]\displaystyle{ \frac{\mathrm{d}P}{\mathrm{d}v}=0 \qquad \text{and} \qquad \frac{\mathrm{d}^2P}{\mathrm{d}v^2}=0 }[/math] the cubic virial equation can be solved to yield: [math]\displaystyle{ B = -v_c , }[/math] [math]\displaystyle{ C = \frac{v_c^2}{3}, }[/math] and [math]\displaystyle{ Z_c = \frac{P_c v_c}{RT_c} = \frac 1 3. }[/math] [math]\displaystyle{ Z_c }[/math] is therefore 0.333, comparing to 0.375 from Van der Waals equation of state.

Between the critical point and the triple point is the saturation region of fluids. In this region, the gaseous phase coexists with the liquid phase under saturation pressure [math]\displaystyle{ P_\text{sat} }[/math], and the saturation temperature [math]\displaystyle{ T_\text{sat} }[/math]. Under the saturation pressure, the liquid phase has a molar volume of [math]\displaystyle{ v_\text{l} }[/math], and the gaseous phase has a molar volume of [math]\displaystyle{ v_\text{g} }[/math]. The corresponding molar densities are [math]\displaystyle{ \rho_\text{l} }[/math] and [math]\displaystyle{ \rho_\text{g} }[/math]. These are the saturation properties needed to compute second and third virial coefficients.

A valid equation of state must produce an isotherm which crosses the horizontal line of [math]\displaystyle{ P_\text{sat} }[/math] at [math]\displaystyle{ v_\text{l} }[/math] and [math]\displaystyle{ v_\text{g} }[/math], on [math]\displaystyle{ T_\text{sat} }[/math]. Under [math]\displaystyle{ P_\text{sat} }[/math] and [math]\displaystyle{ T_\text{sat} }[/math], gas is in equilibrium with liquid. This means that the PρT isotherm has three roots at [math]\displaystyle{ P_\text{sat} }[/math]. The cubic virial equation of state at [math]\displaystyle{ T_\text{sat} }[/math] is: [math]\displaystyle{ P_\text{sat} = RT_\text{sat} \left(1 + B\rho + C\rho^2\right) \rho }[/math] It can be rearranged as: [math]\displaystyle{ 1 - \frac{RT_\text{sat}}{P_\text{sat}} \left(1 + B\rho + C\rho^2\right) \rho = 0 }[/math] The factor [math]\displaystyle{ RT_\text{sat}/P_\text{sat} }[/math] is the volume of saturated gas according to the ideal gas law, and can be given a unique name [math]\displaystyle{ v^\text{id} }[/math]: [math]\displaystyle{ v^\text{id} = \frac{RT_\text{sat}}{P_\text{sat}} }[/math] In the saturation region, the cubic equation has three roots, and can be written alternatively as: [math]\displaystyle{ \left(1 - v_\text{l} \rho \right) \left(1 - v_\text{m} \rho \right) \left(1 - v_\text{g} \rho \right) = 0 }[/math] which can be expanded as: [math]\displaystyle{ 1 - \left(v_\text{l} + v_\text{g} + v_m\right)\rho + \left(v_\text{l} v_\text{g} + v_\text{g} v_\text{m} + v_\text{m} v_\text{l}\right)\rho^2 - v_\text{l} v_\text{g} v_\text{m} \rho^3 = 0 }[/math] [math]\displaystyle{ v_\text{m} }[/math] is a volume of an unstable state between [math]\displaystyle{ v_\text{l} }[/math] and [math]\displaystyle{ v_\text{g} }[/math]. The cubic equations are identical. Therefore, from the linear terms in these equations, [math]\displaystyle{ v_m }[/math] can be solved: [math]\displaystyle{ v_\text{m} = v^\text{id} - v_\text{l} - v_\text{g} }[/math] From the quadratic terms, B can be solved: [math]\displaystyle{ B = -\frac{\left(v_\text{l} v_\text{g} + v_\text{g} v_\text{m} + v_\text{m} v_\text{l}\right)}{v^\text{id}} }[/math] And from the cubic terms, C can be solved: [math]\displaystyle{ C = \frac{v_\text{l} v_\text{g} v_\text{m}}{v^\text{id}} }[/math] Since [math]\displaystyle{ v_\text{l} }[/math], [math]\displaystyle{ v_\text{g} }[/math] and [math]\displaystyle{ P_\text{sat} }[/math] have been tabulated for many fluids with [math]\displaystyle{ T_\text{sat} }[/math] as a parameter, B and C can be computed in the saturation region of these fluids. The results are generally in agreement with those computed from Benedict-Webb-Rubin and Starling equations of state. However, the accuracy in B and C are critically dependent on the measurements of [math]\displaystyle{ v_\text{g} }[/math] and [math]\displaystyle{ P_\text{sat} }[/math], which are very difficult to measure accurately at low temperatures. The measurement errors thus introduced into B and C should be considered when one compares the values thus derived with those derived from second order regression analysis of PρT isotherms.

Gas-liquid-solid equilibrium

The cubic virial equation of state accurately represents the gas-liquid equilibrium of most substance from the critical point down to the triple point, where solid phase starts to appear. It can be extended to account for the gas-liquid-solid equilibrium: [math]\displaystyle{ P = \frac{RT}{V} \left(1+\frac{B}{V}+\frac{C}{V^2}+\frac{U}{V^n}+\frac{W}{V^{2n}}\right) }[/math] In this equation, the first term [math]\displaystyle{ \frac{RT}{V} }[/math] represents the pressure generated by kinetic energy of molecules. The second term [math]\displaystyle{ \frac{RTB}{V^2} }[/math] represents long-range bimolecular attraction. The third term [math]\displaystyle{ \frac{RTC}{V^3} }[/math] represents short-range tri-molecular repulsion. The second term pulls the PVT isotherm down as volume is reduced, while the third term pushes the isotherm up. When the temperature is below the critical point, the PVT isotherm thus has an S-shaped bent which allows a liquid phase to coexist with the prevalent gaseous phase.

Now, if we had a [math]\displaystyle{ \frac{RTU}{V^n} }[/math] term to pull the PVT isotherm down in the liquid phase, and a [math]\displaystyle{ \frac{RTW}{V^{2n}} }[/math] terms to push it back up, a solid phase could be created, as these two terms producing another S-shaped bend between liquid and solid. It was demonstrated that such an S-shaped bend could be synthesized using a [math]\displaystyle{ \delta }[/math]-function like Lorentzian function over a van der Waals equation of state. Such an equation of state is difficult to manipulate mathematically. In contrast, a virial equation is easier to handle.

Argon is used to evaluate realistically this extended virial equation for gas-liquid-solid equilibrium. Data will be analyzed in reduced forms. All PVT variables are scaled by their respective critical values. It is expect from the principle of corresponding states that the results would apply to other well behaved fluids. The relevant data of argon are summarized in the following table:

Property of Argon	Value	Reduced value
Critical point volume (dm3/mol)	0.07459	1
Critical point temperature (K)	150.687	1
Critical point pressure (MPa)	4.863	1
Critical compressibility	0.291	0.291
Triple point vapor molar volume (dm3/mol)	9.853	132.1
Triple point liquid molar volume (dm3/mol)	0.0282	0.378
Triple point solid molar volume (dm3/mol)	0.246	0.330
Triple point temperature (K)	83.8058	0.553
Triple point pressure (MPa)	0.06889	0.0142

When the variables P, V, and T are replaced by their reduced equivalents, [math]\displaystyle{ p = \frac{P}{P_c} }[/math], [math]\displaystyle{ v = \frac{V}{V_c} }[/math], and [math]\displaystyle{ t = \frac{T}{T_c} }[/math], the virial equation takes the following form: [math]\displaystyle{ p = \frac{t}{vZ_c}\left(1-\frac{b}{v}+\frac{c}{v^2}-\left(\frac{v_u}{v}\right)^n+\left(\frac{v_w}{v}\right)^{2n}\right) }[/math] where [math]\displaystyle{ Z_c = \frac{P_c V_c}{RT_c} }[/math], [math]\displaystyle{ b = \frac{B}{V_c} }[/math], [math]\displaystyle{ c = \frac{C}{V_c^2} }[/math], [math]\displaystyle{ U = v_u^n }[/math], and [math]\displaystyle{ W=v_w^{2n} }[/math]. We will be concerned mostly with condition at the triple point of argon, where b = 3.424 and c = 1.152 from an earlier study.

[math]\displaystyle{ V_w }[/math] must be slightly larger than the volume of solid argon, 0.33, and [math]\displaystyle{ v_u }[/math] must be between the volumes of liquid and solid argon. Initially, [math]\displaystyle{ v_w }[/math] is set to the volume of solid, to produce the last sharply rising edge of the isotherm where solid phase appears at very low volume. The exponential n must be then determined so that the valley in the n-2n potential must fit between the volumes of solid (0.33) and liquid (0.378). After the exponential n is determined, the value of [math]\displaystyle{ v_u }[/math] can be adjusted to satisfy the Gibbs Rule, which requires that the Gibbs free energy of liquid phase and that of solid phase must be equal under the triple point temperature and pressure.

To produce a solid phase in argon, the exponential value n must be very large, larger than 20; otherwise, the PVT isotherm would not bend to an S-shape between liquid and solid. The best estimation is that n = 30, [math]\displaystyle{ t_t=0.553 }[/math], [math]\displaystyle{ v_w=0.335 }[/math], and [math]\displaystyle{ v_u=0.3443 }[/math]. The isotherm is shown in the right figure, in which three virial terms are plotted separately for clarity: [math]\displaystyle{ p_1 = \frac{t_t}{vZ_c}\left(1-\frac{b}{v}+\frac{c}{v^2}\right) }[/math] [math]\displaystyle{ p_2 = \frac{t_t}{vZ_c}\left(\frac{v_u}{v}\right)^n }[/math] [math]\displaystyle{ p_3 = \frac{t_t}{vZ_c}\left(\frac{v_w}{v}\right)^{2n} }[/math] [math]\displaystyle{ p = p_1-p_2+p_3 }[/math]

PVT isotherm of argon with n = 30.

In this figure, [math]\displaystyle{ p_1 }[/math] represents the sum of the first three virial terms, of the cubic virial equation, and it shows the behavior of argon in its gaseous and liquid phases. [math]\displaystyle{ p_2 }[/math] represents the potential contributed from the [math]\displaystyle{ \left(\frac{v_u}{v}\right)^n }[/math] term, and [math]\displaystyle{ p_3 }[/math] represents the contributions from the [math]\displaystyle{ \left(\frac{v_w}{v}\right)^{2n} }[/math] term. When n < 30, [math]\displaystyle{ p_2 }[/math] would interfere with [math]\displaystyle{ p_1 }[/math] and lower significantly the volume of liquid.

It is surprising that n has to be greater than 20 to produce a solid phase. It is well known that the Lenard-Jones 6-12 potential can be used to compute the second virial coefficient from quantum mechanical principles. Quantum mechanical reasoning relates the second virial coefficient to bimolecular attraction, and the third coefficients to tri-molecular repulsion, etc. In the liquid phase of argon, one atom is surrounded by 12 nearest neighbors and up to 32 next-to-nearest neighbors. In the solid phase, all atoms are locked in place, and the number of interacting neighbors is infinite. Therefore, n = 30, or even greater, is reasonable.

As a result, we have a virial equation of state to describe quantitatively the gas-liquid-solid equilibrium for argon, and all fluids which observe the principle of corresponding states, at its triple point: [math]\displaystyle{ p = \frac{t_t}{vZ_c}\left(1-\frac{3.424}{v}+\frac{1.152}{v^2}-\left(\frac{0.3443}{v}\right)^{30}+\left(\frac{0.3350}{v}\right)^{60}\right) }[/math]

Virial equations of state

From Benedict-Webb-Rubin and Starling equations of state, it was determined that the best virial equation of state should take the form of [math]\displaystyle{ Z=1+b\rho_r+c\rho_r^2+f\rho_r^5. }[/math] The second and third virial coefficients in this equation can be computed from experimental PρT data using a linear regression. After removing the first three virial terms, the residue in compressibility factor Z could then be used to obtain the sixth virial coefficient.

See also

• Virial theorem
• Statistical mechanics
• Equation of state