**3. EoS models in the thermodynamic perturbation theory context**

Unlike macroscopic EoS models such as those described in the previous section, major efforts have been oriented toward the development of microscopic approaches based on statistical mechanics where the individual behavior of each particle in a fluid substance is considered. The repulsive and attractive forces are handled separately and combined to provide a description of the thermodynamic properties of fluids through methods based on statistical physics.

The basic idea is to study the microscopic behavior of a set of molecules by considering many instances of the set, each one corresponding to one possible state. This ensemble is described through the statistical properties averaged over all possible states. The basic components for this task is the pair potential function *u r*ð Þ and the pair correlation function *g r*ð Þ respectively, both functions of the distance *r* away from the center of some molecule. By defining them one can generate expressions to compute the system free energy and eventually all thermodynamic properties of interest [3].

Arriving to the energy expression while starting from *u r*ð Þ and *g r*ð Þ is a very complex task from the mathematical treatment point of view. Complex expressions of the two functions might correspond to more accurate description of the molecules dynamics but they also lead to intractable mathematical expressions. For this task perturbation theory has greatly enhanced the derivation of EoS models by firstly utilizing known closed form solutions for simple reference functions. Subsequently, the small changes between the accurate *u r*ð Þ and *g r*ð Þ functions and the reference ones are treated in a very elegant way by means of perturbation theory thus leading to approximate closed form solutions for complex pair functions [12].

### **3.1 The correlation function formalism to derive EoS models**

Consider a thermodynamically large system comprising of a fixed number of molecules, at fixed temperature and volume, which is allowed to exchange heat with the environment. Subsequently, consider a collection of many such probable systems forming what is known as a *canonical ensemble*. The aggregate thermodynamic properties of such systems can be described as functions of the statistic properties of the ensemble. For this task the *canonical partition function* is defined by

$$Q = \sum \exp\left(-\beta \mathcal{E}\_i\right),\tag{16}$$

where *Ei* corresponds to the energy of each possible microstate, *β* ¼ 1*=kBT* is the thermodynamic beta and *kB* is the Boltzmann constant. Note that *Q* is dimensionless and as it will be shown later it relates macroscopic thermodynamic properties of the system to the energy of the microscopic systems forming the ensemble.

For a system comprising of N identical molecules the partition function *QN*ð Þ *V*, *T* at given volume and temperature is given by [3].

$$Q\_N(V, T) = \frac{Z\_N(V, T)}{N!\Lambda^3},\tag{17}$$

where

$$Z\_N(V, T) = \int\_N \exp\left(-\beta U\_N\left(\mathbf{r}^N\right)\right) d\mathbf{r}^N, \quad \Lambda = \sqrt{\frac{h^2}{2\pi mk\_B T}},\tag{18}$$

and *ZN*ð Þ *V*, *T* is known as the configuration integral. It is easy to show that if the system potential energy *UN* is assumed to be zero then the configuration integral *ZN* simplifies to the system volume and the application of the related partition function leads simply to the ideal gas law. On the other hand, when *UN* 6¼ 0, it is often represented by a sum of pair-wise potentials, i.e.

$$U(\mathbf{r}^N) = \sum\_{i} \sum\_{j>i} u\left(r\_{\vec{\eta}}\right). \tag{19}$$

*Perturbation Theory and Phase Behavior Calculations Using Equation of State Models DOI: http://dx.doi.org/10.5772/intechopen.93736*

Various pair potential models *u r*ð Þ have been presented with the hard-sphere, the square well and the Lennard-Jones being the most pronounced ones [12]. The hard sphere model assumes that the particles are perfect spheres of diameter *σ*, the potential at distances less than the sphere diameter is equal to infinite (hence the "hard" sphere) and zero beyond that. Therefore, *u r*ð Þ and the corresponding Boltzmann factor are given by

$$u(r) = \begin{cases} \Leftrightarrow & r < \sigma \\ 0 & r > \sigma \end{cases}, \quad \exp\left(-\beta u(r)\right) = \begin{cases} 0 & r < \sigma \\ 1 & r > \sigma \end{cases}.\tag{20}$$

The square-well model [13] further allows for a negative value at some distance beyond the hard sphere diameter:

$$u(r) = \begin{cases} \infty & r < \sigma \\ -\varepsilon & \sigma < r < \chi\sigma \\ 0 & r > \chi\sigma \end{cases} \tag{21}$$

The Lennard-Jones model [14] offers the advantage of being defined by a continuous function of the distance *r*:

$$u(r) = 4e\left[\left(l/r\right)^{12} - \left(l/r\right)^{6}\right].\tag{22}$$

In the equations above σ is the sphere diameter, parameter *γ* is used to scale the well width, *l* is the length parameter and *ε* is the energy parameter. A detailed description on how to use the hard-sphere model pair potential function to develop an EoS model is given in Section 3.3.

#### **3.2 Derivation of fluid properties for specific pair functions**

By selecting the pair potential model *u r*ð Þ and incorporating it the configuration integral *ZN* and eventually to the canonical partition function expression, internal energy can be obtained by noting that

$$E = k\_B T^2 \frac{\partial}{\partial T} Q\_N(V, N). \tag{23}$$

By utilizing the pair-wise potential energy model of Eq. (19) it can be shown that

$$E = \frac{3}{2} Nk\_B T + 2\pi\rho N \int\_0^\infty u(r)g(r)r^2 dr. \tag{24}$$

Therefore, internal energy can be obtained as a function of the particle properties *u r*ð Þ and *g r*ð Þ. Clearly, the first term corresponds to the kinetic energy of the particles, that is the ideal gas contribution of the system.

Using similar arguments, pressure can be obtained by as the volume derivative of the configuration integral *ZN*, that is

$$p = k\_B T \frac{\partial}{\partial V} Z\_N(V, N) = k\_B T \frac{N}{V} - \frac{2\pi}{3} \rho^2 \int\_0^\infty r^3 \frac{\partial u(r)}{\partial r} g(r) dr. \tag{25}$$

Again, the first term corresponds to the ideal gas pressure term. The system Helmholtz energy is defined by

$$H = E - TS = -k\_B T \ln Q\_N(V, T). \tag{26}$$

The chemical potential corresponds to the energy required to add one more particle in the collection it is given by

$$\mu = H(V, T, N) - H(V, T, N - 1) = \frac{\partial H}{\partial \mathbf{N}}\Big|\_{V, T},\tag{27}$$

and eventually

$$
\mu = k\_B T \ln \rho \Lambda^3 + 4\pi \rho \int\_0^1 \int\_0^\infty r^2 u(r) g(r, \lambda) dr d\lambda. \tag{28}
$$

Given the expression above, entropy can be obtained by

$$\mathcal{S} = \frac{E - H}{T}.\tag{29}$$

#### **3.3 The hard-sphere model**

The generic fluid properties expressions derived in the previous section are now applied to the hard sphere model for the pair potential energy. By noting that the derivative of the Boltzmann factor of the hard-sphere model is simply the Dirac delta function [15] and replacing it to the generic properties' expressions of the previous paragraph, it follows for pressure that

$$p = p^{lG} + p^{EX} = \rho k\_B T + \rho k\_B T \frac{4\eta - 2\eta^2}{\left(1 - \eta\right)^3} = \rho k\_B T \frac{\mathbf{1} + \eta + \eta^2 + \eta^3}{\left(1 - \eta\right)^3},\tag{30}$$

where the pressure is now split into the ideal gas and the excess part and the packing function *η* which corresponds to the ratio of the particles volume over the total one is given by

$$\eta = \frac{1}{V} N \frac{4\pi}{3} \left(\frac{\sigma}{2}\right)^3 = \frac{\pi}{6} \rho \sigma^3. \tag{31}$$

The expressions for the other properties of interest are obtained similarly and they are given by

$$\begin{split} H &= H^{\text{IG}} + H^{\text{EX}} = Nk\_B T \left( \ln \rho \Lambda^3 - 1 \right) + Nk\_B T \frac{4\eta - 3\eta^2}{\left(1 - \eta\right)^2} \\ S &= S^{\text{IG}} + S^{\text{EX}} = -Nk\_B \left( \ln \rho \Lambda^3 - \frac{5}{2} \right) - Nk\_B \frac{4\eta - 3\eta^2}{\left(1 - \eta\right)^2} \\ \mu &= \mu^{\text{IG}} + \mu^{\text{EX}} = k\_B T \left( \ln \rho \Lambda^3 - 1 \right) + k\_B T \frac{1 + 5\eta - 6\eta^2 + 2\eta^3}{\left(1 - \eta\right)^3} .\end{split} \tag{32}$$

#### **3.4 Thermodynamic perturbation theory**

Although the hard-sphere pair potential model allows for an explicit calculation of thermodynamic properties of interest, its results are not that accurate mostly due to the inherent simplicity of the hard-sphere model itself. Nevertheless, many

*Perturbation Theory and Phase Behavior Calculations Using Equation of State Models DOI: http://dx.doi.org/10.5772/intechopen.93736*

researchers have pointed out that the comparison between the experimental structure factor and the one obtained computationally from the hard-sphere model indicates that the two curves are quite close to each other. To get a better match a more complex pair potential model could be sought which, however, would inevitably lead to mathematically intractable expressions of the properties. Alternatively, perturbation methods can be applied to the original simple hard-sphere model to add thermodynamic complexity under controlled extra computational burden.

The idea, firstly presented by Zwanzig [16], is to divide the total potential energy into two terms, *U*<sup>0</sup> and *Up* respectively, where the first term corresponds to a reference system and the second one corresponds to the perturbation, which needs to be significantly smaller than the reference one for the perturbation method to be applied successfully. The total energy is then given by

$$U = U\_0 + \lambda U\_p.\tag{33}$$

The perturbation parameter *λ* allows for various mixtures of *U*<sup>0</sup> and *Up* whereas the original fluid energy is obtained for *λ* ¼ 1. By replacing that expression to the configuration integral we obtain

$$Z\_N(V, T) = Z\_N^{(0)}(V, T) \langle \exp\left(-\beta \lambda U\_p\right) \rangle\_0,\tag{34}$$

where the h i*:* operator denotes the statistical average of the reference system. Replacing Eq. (34) to the expression for the Helmholtz energy we obtain

$$-\beta H = \ln \frac{Z\_N(V, T)}{N! \Lambda^{3N}} + \ln \left< \exp \left( -\beta \lambda U\_1 \right) \right>\_0 = -\beta H\_0 - \beta H\_p,\tag{35}$$

where the first term corresponds to a multiple of the Helmholtz energy of the reference system and the second term accounts for the energy of the perturbation. By combining the Taylor expansion forms of the exponential term and of the logarithmic term we obtain

$$-\beta \mathcal{H}\_p = \ln \left< \exp \left( -\beta \lambda U\_p \right) \right>\_0 = -\left( \lambda \beta \right) \mathbf{c}\_1 + \left( \lambda \beta \right)^2 \mathbf{c}\_2 - \left( \lambda \beta \right)^3 \mathbf{c}\_3 + \dots,\tag{36}$$

where

$$\begin{aligned} c\_1 &= \left< U\_p \right>\_0 \\ c\_2 &= \frac{1}{2!} \left( \left< U\_p^2 \right>\_0 - \left< U\_p \right>\_0^2 \right) \\ c\_3 &= \frac{1}{3!} \left( \left< U\_p^3 \right>\_0 - 3 \left< U\_p \right>\_0 \left< U\_p^2 \right>\_0 + 2 \left< U\_p \right>\_0^3 \right). \end{aligned} \tag{37}$$

The treatment above has allowed the energy to be described by simpler expressions of the perturbation energy term based on the *c*1,*c*2,*c*<sup>3</sup> parameters. Therefore, to get the full energy expression one needs to choose the *Up* model, compute the values of the *c*1,*c*2,*c*<sup>3</sup> parameters and replace then in Eq. (36) while setting *λ* ¼ 1. All thermodynamic properties of interest can then be computed as functions of the energy as shown in Section 3.2.

The beauty of the perturbation theory is that although the calculation of the *c*1,*c*2,*c*<sup>3</sup> parameters, which have been introduced by the application of perturbation theory and the assumption of a simple reference system, is not an easy task it still is significantly easier than replacing a complex pair potential and pair correlation function and running mathematical operations in Eq. (18).

As an example application of perturbation theory in statistical mechanics based thermodynamics, consider the use of the model obtained by perturbation theory to generate the Van der Waals EoS, this time from a statistical mechanics point of view instead from the classic macroscopic one. Firstly, let us state the following assumptions:


By introducing those assumptions to Eq. (37) the calculation of coefficient *c*<sup>1</sup> simplifies to

$$\mathcal{L}\_1 = \frac{\rho^2}{2} \int d\mathbf{r} \int u\_p(\mathbf{r}) \mathbf{g}\_o(\mathbf{r}) d\mathbf{r}.\tag{38}$$

To proceed we further need to introduce the following assumptions:


and we end up with

$$\sigma\_1 = -a\rho N, \quad a = -2\pi \int\_{\sigma}^{\infty} u\_p(r)r^2 dr. \tag{39}$$

Finally, by utilizing first order approximation only (up to *c*1), replacing *c*<sup>1</sup> in the free energy Eq. (36) and differentiating over volume to obtain pressure (i.e. *<sup>p</sup>* ¼ �*∂H=∂V*<sup>j</sup> *<sup>T</sup>*) the well known Van der Waals equation is obtained

$$p = \frac{Nk\_B T}{V - Nb} - a \frac{N^2}{V^2}.\tag{40}$$

Interestingly, the perturbed energy term *up* has not been defined explicitly but it has been incorporated into the EoS *a* parameter. From a perturbation theory point of view, the accuracy of the Van der Waals equation of state can be improved by further considering the *c*<sup>2</sup> term, the calculation of which, however, is quite more complicated.
