Thermodynamics of Phase Transitions for Pure Substances and Mixtures

Jan Helm

Department of Electrical Engineering, Technical University Berlin, Berlin, Germany.

DOI: 10.4236/jmp.2026.175027                

Abstract

This paper presents in a concise way the mathematical models and experimental data of phase transitions of pure substances and solutions with new ansatz and calculation methods. In chap. 2 and chap. 3 we deal with the theoretical basics of solutions and mixtures. In chap. 4, the theory of equations-of-state is formulated, and in chap. 5 the calculation results for four pure substances (benzene, ethanol, argon, carbon dioxide) are presented and commented. In chap. 6, the equation-of-state and mixing rules for mixtures are formulated, and the calculation results for an example solution (benzene-ethanol) are discussed and compared with measured values. The paper introduces two novel methods: 1) Exact algebraic solution for phase diagrams based on Peng-Robinson and Mie-Grueneisen equation-of-state; 2) A new ansatz for mixture phase diagrams based on the weighted sum of partial pressures.

Keywords

Gibbs Free Energy, Helmholtz Free Energy, Enthalpy, Equation-of-State, Phase Transition, Van-Der-Waals Equation-of-State, Peng-Robinson Equation-of-State, Mie-Grueneisen Equation-of-State, Maxwell Rule, Saturation Curve, Binary Mixture

1. Introduction

Equations-of-state (eos) are an area, which is of considerable theoretical, but also technological importance.

Eos for fluid-gas phase are relatively well-understood, starting with the famous vdWaals equation.

There are several cubic extensions of vdWaals eos, which are in satisfactory agreement with measurements.

The best is arguably Peng-Robinson eos with 3 parameters and a precision of about 10%.

There exist more precise quartic eos (Shah’s eos), and polynomial-exponential eos.

Beside this, there are individual numerical approximations for selected substances, e.g. water.

Eos for solids are few, most prominent is Mie-Grueneisen eos.

For saturation curves however, the situation is quite different.

For fluid-gas saturation curves, there is only one closed solution, namely Lekner’s solution based on the vdWaals eos. For phase diagrams (solid-fluid-gas), there are only numerical parametric approximations for selected substances, based on direct measurement.

The Maxwell-Gibbs equation, which yields the exact solution for the saturation curve, is known in theory as the Maxwell-rule, but has not been solved for the Peng-Robinson and the Mie-Grueneisen eos.

In this paper, the Maxwell-Gibbs equation is solved for both fluid-gas and the solid phase in the form p(E_th), based on Peng-Robinson and Mie-Grueneisen eos, and is used to calculate the triple point and the phase diagram

(E_th = k_B T thermal energy, is used throughout this paper in place of temperature T).

As for binary mixtures, there exist only purely phenomenological ansatzes, there is no general theoretical treatment based on molecular data of the components.

In this paper we formulate a theoretical basis for binary solutions, based on the weighted sum of partial eos pressures, and including the 1-2-interaction of the components (i.e. non-ideal and irregular solutions).

2. Basics of Classical and Quantum Statistics

Statistical mechanics describes the thermodynamic behavior of large systems ([1]-[3]).

Key features of a thermodynamic ensemble are its partition function (probability distribution) and its macroscopic function (extremal variable of the ensemble: e.g. entropy = max, free energy = min). They are functions of the thermodynamic state variables: temperature, volume pressure, number of particles, chemical potential [4].

Fundamental variable is: partition function Z

$Z={\displaystyle \sum _i\exp \left(-\frac{E_i}{k_{B}T}\right)}$ classical Boltzmann statistics (1a)

$Z={\displaystyle \sum _i\frac{1}{\exp \left(\frac{E_i}{k_{B}T}\right)+1}}$ quantum fermion Fermi-Dirac statistics (1b)

$Z={\displaystyle \sum _i\frac{1}{\exp \left(\frac{E_i}{k_{B}T}\right)-1}}$ quantum boson Bose-Einstein statistics (1c)

Important thermodynamic variables mean energy U, Helmholtz free energy F, Gibbs energy G, entropy S, heat capacity C_v.

$U=\langle E\rangle =-\frac{\partial \log Z}{\partial \beta }=k_B{T}^2\frac{\partial \log Z}{\partial T}$ mean energy (2a)

$F=\langle E\rangle -TS=-k_{B}T\log Z$ Helmholtz free energy (2b)

$G=F+PV=\langle E\rangle -TS+PV=-\frac{1}{\beta }\frac{\partial \left(V\log Z\right)}{\partial V}$ Gibbs energy, (2c)

with beta parameter $\beta =\frac{1}{k_{B}T}$

$S=\frac{\partial }{\partial T}\left(k_{B}T\log Z\right)$ entropy (2d)

$C_v=\frac{\partial \langle E\rangle }{\partial T}=\frac{1}{k_B{T}^2}\langle {\left(\Delta E\right)}^2\rangle =k_{B}T\left(2\frac{\partial \log Z}{\partial T}+T\frac{{\partial }^2\log Z}{\partial T^2}\right)$ heat capacity (2e)

for

$Z\approx W\exp \left(-\frac{\langle E\rangle }{k_{B}T}\right)$ , $S=\frac{\partial }{\partial T}\left(k_{B}T\log Z\right)\approx \frac{\partial }{\partial T}\left(-\langle E\rangle +k_{B}T\log W\right)=k_B\log W$ (2f)

which gives the famous Boltzmann approximate formula

$S\approx k_B\log W$ (2g)

2.1. The Three Most Important Ensembles

The microcanonical ensemble describes [5] a system with fixed energy and fixed number of particles.

The canonical ensemble describes a system of fixed number of particles that is in thermal equilibrium with a heat bath of a fixed temperature.

The grand canonical ensemble describes a system with non-fixed particle numbers that is in thermal and chemical equilibrium with a thermodynamic reservoir with fixed temperature and fixed number of particles.

Thermodynamic classical and quantum ensembles with their partition function and the macroscopic function [6] are shown in Table 1.

Table 1. Thermodynamic classical and quantum ensembles.

	Thermodynamic classical ensembles
	Microcanonical	Canonical	Grand canonical
Fixed variables	N, V, E	N, V, T	μ, V, T
Microscopic features	Number of microstates W, with width ω $W={\displaystyle \sum _{k}f\left(\frac{H_k-E}{\omega }\right)}$ $f\left(x\right)=\exp \left(-\pi x^2\right)$	Canonical partition function $Z={\displaystyle \sum _k\exp \left(-\frac{E_k}{k_{B}T}\right)}$ $\begin{array}{l}Tr\left(\exp \left(-\beta H\right)\right)={\displaystyle \sum _k\exp \left(-\beta E_k\right)}\\ =Z\left(\beta \right)\end{array}$	Grand partition function $Z={\displaystyle \sum _k\exp \left(-\frac{E_k-\mu N_k}{k_{B}T}\right)}$ $Z\left(\beta ,{\mu }_1,{\mu }_2,\cdots \right)=Tr\left(\exp \left(\beta \left({\displaystyle \sum _k{\mu }_k{N}_k}-H\right)\right)\right)$
Probability and density matrix	$\widehat{\rho }=\frac{1}{W}{\displaystyle \sum _{k}f\left(\frac{H_k-E}{\omega }\right)}|{\psi }_k\rangle \langle {\psi }_k|$	$\rho =\frac{{\text{e}}^{-\beta H}}{Tr{\text{e}}^{-\beta H}}$ $p\left(E_m\right)=\frac{\exp \left(-\beta E_m\right)}{{\displaystyle \sum _k\exp \left(-\beta E_k\right)}}$	$\rho =\frac{{\text{e}}^{-\beta \left(H-{\displaystyle \sum _i{\mu }_i{N}_i}\right)}}{Tr{\text{e}}^{-\beta \left(H-{\displaystyle \sum _i{\mu }_i{N}_i}\right)}}$ $P\left(E_m\right)=\frac{{\text{e}}^{-\beta \left(E_m-{\displaystyle \sum _i{\mu }_i{N}_i}\right)}}{{\displaystyle \sum _n{\text{e}}^{-\beta \left(E_n-{\displaystyle \sum _i{\mu }_i{N}_i}\right)}}}$
Minimal principle	Boltzmann entropy $S=k_B\log W=\max$ v.Neumann entropy $S=-k_{B}Tr\left(\rho \ln \rho \right)=\max$	Helmholtz free energy $F=-k_{B}T\log Z=\min$	Grand potential $\Omega =-k_{B}T\log Z=\min$

The macroscopic function obeys the minimum principle of statistics:

• the macroscopic function attains an extremum in equilibrium

• microcanonical: entropy is maximal

• canonical: free energy is minimal

• grand canonical: grand potential is minimal

The partition function describes the statistical properties of a system in thermodynamic equilibrium^.

Partition functions are functions of the thermodynamic state variables.

For the canonical discrete ensemble the partition function reads

classical $Z={\displaystyle \sum _k{g}_k\exp \left(-\beta E_k\right)}$ , $\beta =\frac{1}{k_{B}T}$ (3a)

quantum $Z=Tr\left(\exp \left(-\beta \widehat{H}\right)\right)$ or $Z={\displaystyle \sum _k\frac{g_k}{\exp \left(\frac{E_k}{k_{B}T}\right)\pm 1}}$ (fermion+, boson−) (3b)

probability of state s reads $p_s=\frac{1}{Z}\exp \left(-\beta E_s\right)$ (3c)

The three probability distributions for the average number in a state ε, $N\left(\epsilon \right)$ with degeneracy g(ε) and chemical potential μ are derived from the maximization of number of states $W\left({\nu }_i\right)$ with occupation numbers ${\nu }_i$ of N particles for states with degeneracy g_r and energy levels ${\eta }_r$ $r=1,2,\cdots ,N^{\prime }$ under energy condition

${\displaystyle \sum _r{\nu }_r{\eta }_r=E}$ and particle number condition ${\displaystyle \sum _r{\nu }_r=N}$ .

• Fermions $S=n+1/2$ (Pauli principle)

Fermi-Dirac statistics $W\left({\nu }_i\right)=\underset{r}{\Pi }\frac{g_r!}{{\nu }_r!\left(g_r-{\nu }_r\right)!}$ (4a)

follows from Pauli-principle ${\nu }_r=0,1$ because the wave function is antisymmetric in particles $\Psi \left(\cdots {\nu }_i\cdots {\nu }_k\cdots \right)=-\Psi \left(\cdots {\nu }_k\cdots {\nu }_i\cdots \right)$ :

$N\left(\epsilon \right)=\frac{g\left(\epsilon \right)}{\exp \left(\frac{\epsilon -\mu }{k_{B}T}\right)+1}$ (4b)

• Bosons $S=n$

Bose-Einstein statistics $W\left({\nu }_i\right)=\underset{r}{\Pi }\frac{\left(g_r+{\nu }_r-1\right)!}{{\nu }_r!\left(g_r-1\right)!}$ (4c)

follows from ${\nu }_r=0,1,2,\cdots$ because the wave function is symmetric in particles $\Psi \left(\cdots {\nu }_i\cdots {\nu }_k\cdots \right)=\Psi \left(\cdots {\nu }_k\cdots {\nu }_i\cdots \right)$ :

$N(\epsilon )=\frac{g(\epsilon )}{\exp \left(\frac{\epsilon -\mu }{k_{B}T}\right)-1}$ (4d)

• Boltzmann-Maxwell statistics

$N\left(\epsilon \right)=\exp \left(-\frac{\epsilon }{k_{B}T}\right)$ (4e)

is derived from both FD and BM statistics for high temperature $\exp \left(\frac{\epsilon }{k_{B}T}\right)\gg 1$ .

The second principle of thermodynamics states that $\frac{\text{d}S}{\text{d}t}\ge 0$ entropy is non-decreasing in time t, or (Lorentz-invariant) $\frac{\text{d}S}{\text{d}\tau }\ge 0$ entropy is non-decreasing in proper time τ.

2.2. Statistical Mechanics Basics

[5]

Boltzmann theorem: The probability p_α of the system being found in the microstate α is proportional to

$p_{\alpha }=\exp \left(-\frac{E_{\alpha }}{k_{B}T}\right)$ (5a)

Gibbs probability is $p_{\alpha }=\frac{1}{Z}\exp \left(-\frac{E_{\alpha }}{k_{B}T}\right)$ , (5b)

partition function $Z={\displaystyle \sum _{\alpha }\exp \left(-\frac{E_{\alpha }}{k_{B}T}\right)}$ , where the beta-parameter is $\beta =1/k_{B}T$ , the thermal average (X) of any property X of the system is then

$\langle X\rangle ={\displaystyle \sum _{\alpha }{X}_{\alpha }{p}_{\alpha }}=\frac{1}{Z}{\displaystyle \sum _{\alpha }{X}_{\alpha }\exp \left(-\frac{E_{\alpha }}{k_{B}T}\right)}$ (5c)

2.3. Thermodynamic Quantities

Important thermodynamic quantities are [7]

• mean energy $U\equiv \langle E\rangle ={\displaystyle \sum _{\alpha }{X}_{\alpha }{p}_{\alpha }}=\frac{1}{Z}{\displaystyle \sum _{\alpha }{E}_{\alpha }\exp \left(-\beta E_{\alpha }\right)}=-\left(\frac{\partial \log Z}{\partial \beta }\right)$ (6a)

• heat capacity $C_v\equiv {\left(\frac{\partial U}{\partial T}\right)}_V=k_B{\beta }^2{\left(\frac{{\partial }^2\log Z}{\partial {\beta }^2}\right)}_V$ , $C_v\text{d}T=T\text{d}S$ , $C_v=-\beta {\left(\frac{\partial S}{\partial \beta }\right)}_V$ (6b)

where entropy reads $S=-k_B\beta {\left(\frac{\partial \log Z}{\partial \beta }\right)}_V+k_B\log Z$ (6c)

from this follows approximately $Z\approx \Omega \exp \left(-\beta U\right)$ and the famous Boltzmann formula

$S\approx k_B\log \Omega$ (6d)

where $\Omega =\Omega \left(U,U+\delta U\right)$ is the number of states,

and with $\delta U=\text{stdev}\left(U\right)$ , $\delta U=T{\left(k_B{C}_v\right)}^{1/2}$ (6e)

• Helmholtz free energy $F\left(T,V,N\right)\equiv U-TS=-\frac{1}{\beta }\log Z$ (6f)

$\text{d}F=-S\text{d}T-p\text{d}V+\mu \text{d}N$

• pressure p_i conjugate to volume V_i is $p\equiv -\left(\frac{\partial F}{\partial V}\right)=\frac{1}{\beta }\frac{\partial \log Z}{\partial V}$ (6g)

• we introduce the average particle distance λ, and the specific volume $v=\frac{V}{N}$ , where $v={\lambda }^3$

• equation-of state (eos) $p\left(\lambda ,\beta \right)=-\frac{\partial F}{\partial V}=\frac{1}{\beta }\frac{\partial \log Z}{\partial V}=\frac{1}{\beta }\frac{1}{3{\lambda }^2}\frac{1}{Z}\frac{\partial Z\left(\lambda ,\beta \right)}{\partial \lambda }$ , (6h)

where the ideal gas law is $p=\frac{1}{v\beta }=\frac{1}{{\lambda }^3\beta }$ (6i)

• Gibbs free energy $G\left(T,p,N\right)\equiv F+pV=\mu N=-\frac{1}{\beta }\frac{\partial }{\partial V}\left(V\log Z\right)$ (6j)

$\text{d}G=-S\text{d}T+V\text{d}p+\mu \text{d}N$ , $V=\frac{\partial G}{\partial p}$

• Entropy $S=-k_B\beta {\left(\frac{\partial \log Z}{\partial \beta }\right)}_V+k_B\log Z$ (6k)

$S=-k_B\beta {\left(\frac{\partial \log Z}{\partial \beta }\right)}_V+k_B\log Z$

$S={-\frac{\partial G}{\partial T}|}_{p,N}={-\frac{\partial G}{\partial T}|}_{V,N}$

• Enthalpy $H\equiv U+pV$

$\begin{array}{c}H=U+pV=G-F+U\\ =-\frac{1}{\beta }\frac{\partial }{\partial V}\left(V\log Z\right)+\frac{1}{\beta }\log Z-\left(\frac{\partial \log Z}{\partial \beta }\right)\\ =-\left(\frac{1}{\beta }\frac{\partial \log Z}{\partial V}+\left(\frac{\partial \log Z}{\partial \beta }\right)\right)\end{array}$ (6l)

$\text{d}H=\text{d}U+\text{d}\left(pV\right)$ , $\text{d}H=T\text{d}S+V\text{d}p+\mu \text{d}N=\text{d}G+\text{d}\left(TS\right)$

• compressibility $\kappa \left(\lambda ,\beta \right)=-\frac{1}{V}{\left(\frac{{\partial }^{2}G}{\partial p^2}\right)}_T$ , (6m)

from this follows heat capacity $C_V\left(\lambda ,\beta \right)=-T{\left(\frac{{\partial }^{2}F}{\partial T^2}\right)}_V$ (6n)

• chemical potential of a species ${\mu }_i=\left(\frac{\partial U}{\partial N_i}\right)$ (6o)

• activity of a species (specific = per particle) $a_i=\exp \left(\frac{{\mu }_i-{\mu }_{0,i}}{k_{B}T}\right)$

or equivalently ${\mu }_i={\mu }_{0,i}+k_{B}T\log a_i$ (6p)

3. Basics of Solutions

3.1. Basic Equations

Partition function is $Z=Z\left(V,T,{\mu }_{\alpha }\right)$

with moles n_α of chemical species α, x_α = n_α/n for the mole fraction of component α,

$\text{d}Z=\left(\partial Z/\partial V\right)\text{d}V+\left(\partial Z/\partial T\right)\text{d}T+{\displaystyle \sum _{\alpha }\left(\partial Z/\partial {\mu }_{\alpha }\right)\text{d}{\mu }_{\alpha }}$ (7)

from this follows

$\left(\partial Z/\partial p\right)\text{d}p+\left(\partial Z/\partial T\right)\text{d}T={\displaystyle \sum _{\alpha }{n}_{\alpha }\text{d}{\mu }_{\alpha }}$ (8)

$p\text{d}V-S\text{d}T={\displaystyle \sum _{\alpha }{n}_{\alpha }\text{d}{\mu }_{\alpha }}$ Gibbs-Duhem equation [8]-[10] (9)

The chemical potential of component α is ${\mu }_{\alpha }\equiv \frac{\partial G}{\partial n_{\alpha }}\left(p,T,n_{\alpha }\right)$ (10a)

molar chemical potential of ideal gas is $\mu \left(T,p\right)={\mu }_0\left(T\right)+RT\log \left(p/p_0\right)$ (10b)

where p is the pressure, p₀ is the reference pressure (1 bar), and μ₀ is the reference chemical potential.

The ideal gas eos is for n mols $V=\partial G/\partial p=nRT/p$ (11)

free energy becomes $F=G-pV=n\left({\mu }_0\left(T\right)-RT+RT\log \left(nRT/V{p}_0\right)\right)$ (12)

For perfect gas mixtures

$F_{mix}={\displaystyle \sum _{\alpha }{n}_{\alpha }\left({\mu }_{0,\alpha }\left(T\right)-RT+RT\log \left(n_{\alpha }RT/V{p}_0\right)\right)}$ (12a)

$G_{mix}={\displaystyle \sum _{\alpha }{n}_{\alpha }\left({\mu }_{0,\alpha }\left(T\right)+RT\log \left(p_{\alpha }/p_0\right)\right)}$ (13)

$H_{mix}=\frac{\partial \left(\mu /T\right)}{\partial \left(1/T\right)}=\frac{\partial \left(\mu \beta \right)}{\partial \beta }$ Gibbs-Helmholtz equation [8] (14)

for a perfect mixture

$p=p{\displaystyle \sum _{\alpha }{x}_{\alpha }}={\displaystyle \sum _{\alpha }{p}_{\alpha }}$ Dalton’s law (15)

where partial volume is $V_{\alpha }=\partial {\mu }_{\alpha }/\partial p=RT\partial \log p_{\alpha }/\partial p=RT/p$ (16)

Gibbs free energy of the mixture is $\left(p_{\alpha },{p^{\prime }}_{\alpha }\right)$ is the partial pressure before and after mixing)

${\Delta }_{mix}G=RT{\displaystyle \sum _{\alpha }{n}_{\alpha }\log \left(p_{\alpha }/{p^{\prime }}_{\alpha }\right)}$ (17a)

and the entropy

${\Delta }_{mix}S=-\frac{\partial }{\partial T}{\Delta }_{mix}G=-R{\displaystyle \sum _{\alpha }{n}_{\alpha }\log \left(p_{\alpha }/{p^{\prime }}_{\alpha }\right)}$ (17b)

A perfect solution is defined as having the same Gibbs energy of mixing as the perfect gas mixture

${\Delta }_{mix}S=-R{\displaystyle \sum _{\alpha }{x}_{\alpha }\log \left(x_{\alpha }\right)}$ (17c)

${\Delta }_{mix}\mu =RT{\displaystyle \sum _{\alpha }{x}_{\alpha }\log \left(x_{\alpha }\right)}$ (17d)

equivalently we have the partial pressure of perfect solution

${\mu }_{\alpha }\left(T,p\right)={\mu }_{1,\alpha }\left(T,p\right)+RT\log \left(x_{\alpha }\right)$ (18a)

where ${\mu }_{1,\alpha }\left(T,p\right)$ is the partial chemical potential of pure component.

From μ(sol) = μ(vap) follows

${\mu }_{1,\alpha }+RT\log \left(x_{\alpha }\right)={\mu }_{0,\alpha }+RT\log \left(p_{\alpha }/p_0\right)$ (18b)

so

${\mu }_{1,\alpha }={\mu }_{0,\alpha }+RT\log \left(p_{\alpha }/p_0{x}_{\alpha }\right)={\mu }_{0,\alpha }+RT\log \left(K_{H,\alpha }/p_0\right)$ (18c)

where

$p_{\alpha }=K_{H,\alpha }{x}_{\alpha }$ , $K_{H,\alpha }=\exp \left(\frac{{\mu }_{1,\alpha }-{\mu }_{0,\alpha }}{RT}\right)$ (18d)

is Henry’s constant, independent from $x_{\alpha }$ .

For diluted binary solutions we have for the solvent $x_1\approx 1$ , $x_2\ll 1$ , $p_1=K_{H,1}$ , and $p_1\approx p_1^{*}{x}_1$ Raoult’s law (19)

where $p_1^{*}=p_1^{*}\left(T,V\right)$ is the pressure of pure solvent, is exact for perfect solutions.

Values of Henry’s constant are given for various substances in Table 2.

Table 2. [8] Henry’s constant K_H for dissolved gases at T = 25˚C.

Gas	KH (water) in GPa	KH (benzene) in GPa
CH4 methane	4.185	0.0569
C2H2 acetylene	0.135
C2H4 ethylene	1.155
C2H6 ethane	3.06
air	7.295
N2 nitrogen	8.765	0.239
O2 oxygen	4.438
H2 hydrogen	7.16	0.367
He helium	12.66
CO carbon monoxide	5.79	0.163
CO2 carbon dioxide	1.67	0.0144
H2S hydrogen sulfide	0.055

3.2. Excess Energy, Freezing-Boiling Point Shift

Excess values (compared to ideal solution) are given by [8] [11]

Excess Gibbs energy

$G_E={\displaystyle \sum _{\alpha }{x}_{\alpha }\left({\mu }_{\alpha }-{\mu }_{\alpha }^{*}\right)}$ (20a)

Excess volume

$V_E=V_m-{\displaystyle \sum _{\alpha }{x}_{\alpha }{V}_{\alpha }}$ (20b)

where V_m is the volume of the mixture.

As an example, excess volume V_E for a mixture of water-ethanol [8] is given in Figure 1.

Figure 1. Excess volume V_E for a mixture of water-ethanol.

Freezing point, boiling point, osmotic pressure in binary mixtures

Freezing point depression is $T_f=T^{*}-{\Delta }_{f}T$ ,

we have

${\mu }_1\left(s\right)={\mu }_1\left(l\right)+RT\ln \left(x_1\right)$ (21a)

change in Gibbs energy ${\Delta }_{f}G\equiv {\mu }_1\left(l\right)-{\mu }_1\left(s\right)$ , $-\frac{{\Delta }_{f}G}{RT}=\log \left(x_1\right)$ (21b)

change in enthalpy ${\Delta }_{f}H=R{T}^2\frac{\partial \ln \left(x_1\right)}{\partial T}$ (21c)

from which follows ${\Delta }_{f}T=K_f{x}_2=\frac{R{T}^{*}{}^2}{{\Delta }_{f}H}$ (21d)

Correspondingly, boiling point elevation is $T_b=T^{*}+{\Delta }_{b}T$ ,

change in enthalpy is ${\Delta }_{b}T=K_b{x}_2=\frac{R{T}^{*}{}^2}{{\Delta }_{b}H}$ (21e)

Values of freezing-boiling constants are given for various substances in Table 3.

Table 3. Cryoscopic constants K_f and ebullioscopic constants K_b for some compounds [8].

Compound	Kf (K kg/mol)	Kb (K kg/mol)
acetic acid	3.9	3.07
benzene	5.12	2.53
CS2 carbon disulfide	3.8	2.37
CCl4 carbon tetrachloride	30	4.95
naphthalene	6.94	5.8
phenol	7.27	3.04
water	1.86	0.51
camphor	40

Osmotic (molar) pressure of a binary mixture ${\Delta }_{O}p$ is

${\Delta }_{O}p=c_{O}RT$ , $c_O=\frac{c_2}{x_1{V}_1^{*}+x_2{V}_2^{*}}$ Hoff’s law (22)

Thermodynamic variables for binary mixtures

Given a fluid with coordination number z, interaction energies w₁₁, w₂₂, w₁₂, molecule number N = N₁ + N₂, we have [8] [9] [11]

internal energy change $\Delta U=w_{12}-\left(w_{11}+w_{22}\right)/2=0$ (23a)

entropy change $\Delta S=k_B\log \left(W\right)=-N{k}_B\left(x_1\log x_1+x_2\log x_2\right)$ (23b)

free energy $F=U-TS$ (23c)

the configurational partition function $Z=\exp \left(-\beta F\right)$ reads

$Z=\frac{N!}{N_1!N_2!}\exp \left(-\frac{\beta z}{2}\left(N_1{w}_{11}+\frac{N_2{w}_{22}}{2}\right)\right)$ (23d)

and the free energy becomes

$F=-k_{B}T\log Z_m=z\left(\frac{N_1{w}_{11}}{2}+\frac{N_2{w}_{22}}{2}\right)+k_{B}T\left(N_1\log N_1+N_2\log N_2\right)$ (23e)

with partial free energy

$F_{\alpha }=-k_{B}T\log Z_{\alpha }=\frac{z{N}_{\alpha }{w}_{\alpha \alpha }}{2}$ (23f)

we obtain mixing free energy $\Delta F=F-\left(x_1{F}_1+x_1{F}_1\right)=N{k}_{B}T\left(x_1\log x_1+x_2\log x_2\right)$ (23g)

and the mixing entropy $\Delta S=-\frac{\partial \Delta F}{\partial T}=-N{k}_B\left(x_1\log x_1+x_2\log x_2\right)$ (23h)

the chemical potential ${\mu }_{\alpha }=\frac{\partial F}{\partial N_{\alpha }}=\frac{z{w}_{\alpha \alpha }}{2}+k_{B}T\log x_{\alpha }$ , ${\mu }_{\alpha }-{\mu }_{0,\alpha }=k_{B}T\log x_{\alpha }$ (23i)

The activity λ is given by ${\mu }_{\alpha }=k_{B}T\log {\lambda }_{\alpha }$ , and we obtain $\frac{{\lambda }_{\alpha }}{{\lambda }_{0,\alpha }}=x_{\alpha }$ (23j)

3.3. Ideal and Regular Molecular Solutions

Different types of solutions are given by Table 4.

Table 4. Types of solutions.

ideal solution	$\Delta H=0$	$\Delta S=\Delta S_{ideal}$
athermal solution	$\Delta H=0$	$\Delta S\ne \Delta S_{ideal}$
regular solution	$\Delta H\ne 0$ , $\Delta H=A{x}_1{x}_2$	$\Delta S=\Delta S_{ideal}$
irregular solution	$\Delta H\ne 0$	$\Delta S\ne \Delta S_{ideal}$

Interchange energy is $w\equiv w_{12}-\frac{\left(w_{11}+w_{22}\right)}{2}$ (24a)

N = nN_A, where N_A is the Avogadro number, n number of moles change in interaction inner energy is after mixing

$\Delta U=zN\left(x_1{x}_2{w}_{12}+\frac{x_1{}^2{w}_{11}}{2}+\frac{x_2{}^2{w}_{22}}{2}-\frac{x_1{w}_{11}}{2}-\frac{x_2{w}_{22}}{2}\right)$ (24b)

per mole inner energy $\Delta U_m=\frac{\Delta U}{n}=z{N}_A{x}_1{x}_{2}w$ (24c)

for ideal solution

$U_m=x_1{u}_1+x_2{u}_2=\frac{z\left(x_1{w}_{11}+x_2{w}_{22}\right)}{2}$ , $\Delta H_m=\Delta U_m$ , $G_{ideal}={\displaystyle \sum _{\alpha }\left(n_{\alpha }{\mu }_{0,\alpha }+RT{n}_{\alpha }\log x_{\alpha }\right)}$ (24d)

for regular solution

${\Delta }_{mix}{S}_m=-R{\displaystyle \sum _{\alpha }{x}_{\alpha }\log \left(x_{\alpha }\right)}$ , $S_E={\Delta }_{mix}{S}_m={\Delta }_{mix}{S}_{m,ideal}$ (24e)

${\Delta }_{mix}{G}_m=N_A\left(zw{x}_1{x}_2+k_{B}T\left(x_1\log x_1+x_2\log x_2\right)\right)$ (24f)

excess Gibbs energy $G_E={\Delta }_{mix}G-{\Delta }_{mix}{G}_{ideal}=N_{A}zw{n}_1{n}_2/\left(n_1+n_2\right)$ (24g)

Activity coefficient of regular solutions γ_α

From $G={\displaystyle \sum _{\alpha }{n}_{\alpha }{\mu }_{\alpha }}$ , $n_{\alpha }=x_{\alpha }N/N_A$ , ${\mu }_{\alpha }={\mu }_{0,\alpha }+RT\log \left({\gamma }_{\alpha }{x}_{\alpha }\right)$ (24h)

we obtain for activity coefficient γ_α

$G={\displaystyle \sum _{\alpha }\left(n_{\alpha }{\mu }_{0,\alpha }+n_{\alpha }RT\log \left({\gamma }_{\alpha }{x}_{\alpha }\right)\right)}$ (24i)

which gives excess Gibbs energy

$G_E={\Delta }_{mix}G-{\Delta }_{mix}{G}_{ideal}={\displaystyle \sum _{\alpha }{n}_{\alpha }RT\log {\gamma }_{\alpha }}=N_{A}zw{n}_1{n}_2/\left(n_1+n_2\right)$ (24j)

so

$\partial G_E/\partial n_{\alpha }=RT\log {\gamma }_{\alpha }$ , $\log {\gamma }_1=\beta zw{x}_2^2$ , $\log {\gamma }_2=\beta zw{x}_1^2$ (24k)

Total pressure becomes $p=p_1^{*}{x}_1\exp \left(\beta zw{x}_2^2\right)+p_2^{*}{x}_2\exp \left(\beta zw{x}_1^2\right)$ (24l)

where as usual $p_1^{*}$ and $p_2^{*}$ are the pressure of the pure component.

Phase separation and vapor pressure

The critical point for phase separation can be obtained from

$\frac{\partial {\mu }_1}{\partial x_1}=0$ , $\frac{{\partial }^2{\mu }_1}{\partial x_1^2}=0$ , $\frac{\partial {\mu }_2}{\partial x_2}=0$ , $\frac{{\partial }^2{\mu }_2}{\partial x_2^2}=0$ , (25a)

follows

$x_1=0.5$ , $x_2=0.5$ , ${\beta }_{cr}=2/zw$

$x_1=x_2=0.5$ , ${\left(\beta zw\right)}_{cr}=2$ , $T_{cr}=zw/2k$

activity coefficient ${\gamma }_{cr}^{}=1.649$ , activity $a_{cr}=\gamma x=0.824$ ,

critical pressure $p_{cr}=a_{cr}\left(p_1^{*}+p_2^{*}\right)=0.824\left(p_1^{*}+p_2^{*}\right)$ , (25b)

compared to $p=\left(p_1^{*}+p_2^{*}\right)$ for insoluble liquids (25c)

Regular solutions with correlation functions and volume fractions

• General energy calculation [8] [11]

The volume is

$V=\left(n_1{V}_1^{}+n_2{V}_2^{}\right)$ (26a)

with volume fractions ${\varphi }_i=x_i{V}_i/\left(x_1{V}_1+x_2{V}_2\right)$ (26b)

inner energy in dependence of correlation function is (radial-symmetrical potential u, radial distribution function g)

$U=\langle E_{kin}\rangle +\langle E_{pot}\rangle =\frac{3N{k}_{B}T}{2}+\frac{N^2}{2V}{\displaystyle \underset{0}{\overset{\infty }{\int }}u\left(r\right)g\left(r\right)4\pi r^2\text{d}r}$ (26c)

with average potential energy

$\begin{array}{c}\langle E_{pot}\rangle =V2\pi N_A^2(\frac{{\varphi }_1^2}{V_1^2}{\displaystyle \underset{0}{\overset{\infty }{\int }}{u}_{11}\left(r\right)g_{11}\left(r\right)4\pi r^2\text{d}r}+\frac{{\varphi }_2^2}{V_2^2}{\displaystyle \underset{0}{\overset{\infty }{\int }}{u}_{22}\left(r\right)g_{22}\left(r\right)4\pi r^2\text{d}r}\\ +\frac{2{\varphi }_1{\varphi }_2}{V_1{V}_2}{\displaystyle \underset{0}{\overset{\infty }{\int }}{u}_{12}\left(r\right)g_{12}\left(r\right)4\pi r^2\text{d}r})\end{array}$ (26d)

Subtracting the energy of the separate components

$U_1+U_2=V2\pi N_A^2\left(\frac{{\varphi }_1{\varphi }_2}{V_1{}^2}{\displaystyle \underset{0}{\overset{\infty }{\int }}{u}_{11}\left(r\right)g_{11}\left(r\right)4\pi r^2\text{d}r}+\frac{{\varphi }_1{\varphi }_2}{V_2{}^2}{\displaystyle \underset{0}{\overset{\infty }{\int }}{u}_{22}\left(r\right)g_{22}\left(r\right)4\pi r^2\text{d}r}\right)$ (26e)

gives mixing energy

$\begin{array}{c}{\Delta }_{mix}U=\langle E_{pot}\rangle -\left(U_1+U_2\right)\\ =V2\pi N_A^2(\frac{{\varphi }_1^2}{V_1^2}{\displaystyle \underset{0}{\overset{\infty }{\int }}{u}_{11}\left(r\right)}{g}_{11}\left(r\right)4\pi r^2\text{d}r+\frac{{\varphi }_2^2}{V_2^2}{\displaystyle \underset{0}{\overset{\infty }{\int }}{u}_{22}\left(r\right)}{g}_{22}\left(r\right)4\pi r^2\text{d}r\\ +{\varphi }_1{\varphi }_2(\frac{2}{V_1{V}_2}{\displaystyle \underset{0}{\overset{\infty }{\int }}{u}_{12}\left(r\right)}{g}_{12}\left(r\right)4\pi r^2\text{d}r-\frac{1}{V_1^2}{\displaystyle \underset{0}{\overset{\infty }{\int }}{u}_{11}\left(r\right)}{g}_{11}\left(r\right)4\pi r^2\text{d}r\\ -\frac{1}{V_2^2}{\displaystyle \underset{0}{\overset{\infty }{\int }}{u}_{22}\left(r\right)}{g}_{22}\left(r\right)4\pi r^2\text{d}r))\end{array}$ (26f)

We introduce scaling with energy scale ε and length scale σ,

$u_{ij}\left(r\right)={\epsilon }_{ij}{u}_0\left(r/{\sigma }_{ij}\right)$ , $g_{ij}\left(r\right)=g_0\left(r/{\sigma }_{ij}\right)$

and obtain

${\Delta }_{mix}U=V2\pi N_A^2{\varphi }_1{\varphi }_2\left({\displaystyle \underset{0}{\overset{\infty }{\int }}{u}_0\left(r\right)}{g}_0\left(r\right)4\pi r^2\text{d}r\right)\left(\frac{2{\epsilon }_{12}{\sigma }_{12}^3}{V_1{V}_2}-\frac{{\epsilon }_{11}{\sigma }_{11}^3}{V_1^2}-\frac{{\epsilon }_{22}{\sigma }_{22}^3}{V_2^2}\right)$ (26g)

Experimental approximation

Redlich-Kister expansion of excess Gibbs energy is a Taylor series in $\left(x_1-x_2\right)$ with coefficients (A, B, C) [8] [11]:

$G_E=x_1{x}_2\left(A+B\left(x_1-x_2\right)+C{\left(x_1-x_2\right)}^2+\cdots \right)$ (27a)

for ideal solution $A=0$ , $B=0$ , $C=0$ , $G_E=0$ , ${\gamma }_1=1$ , ${\gamma }_2=1$

$A=B=C=\cdots =0$ , $G^E=0$ , $RT\ln {\gamma }_1=RT\ln {\gamma }_2=0$ and ${\gamma }_1={\gamma }_2=1$ (27b)

for regular solution $B=0$ , $C=0$ , $G_E=A{x}_1{x}_2$ , $RT\log {\gamma }_1=A{x}_2^2$ , $RT\log {\gamma }_2=A{x}_1^2$ (27c)

for deviation from regular $C=0$ ,

$G_E=x_1{x}_2\left(A+B\left(x_1-x_2\right)\right)$ Margules equation (27d)

$RT\log {\gamma }_1=A{x}_2^2$ , $RT\log {\gamma }_2=A{x}_1^2$

In Table 5 below, are given values of Redlic-Kister parameters for various substances.

Table 5. Excess Gibbs function parameters for various solutions [8].

Components	T0 (K)	A/RT0	B/RT0	C/RT0
ethanol/methylcyclohexane	305	2.118	−0.239	0.375
methylcyclohexane/acetone	318	1.6907	−0.0001	0.1832
pyridine/acetone	303	0.1919	0.00050	0.0075
chloroform/furan	303	−0.1083	−0.0177	0.0071
pyridine/chloroform	303	−1.0271	0.2270	0.0930
chloroform/1-4-dioxan	303	−1.2006	−0.4131	0.0318

4. Equations-of-State

4.1. Cubic Fluid-Gas Equations-of-State

vdWaals eos

The well-known vdWaals equation-of-state (eos) for real gas reads in molar variables

$p=\frac{RT}{\left(V_m-b\right)}-\frac{a}{V_m^2}$ (28a)

and in specific (per particle) variables

$p=\frac{1}{\left(v-b_1\right)\beta }-\frac{a_1}{v^2}$ (28b)

with critical parameters

$\frac{1}{{\beta }_c}=k{T}_c=\frac{8a}{27b}$ , $v_c=3b$ , $p_c=\frac{a}{27b^2}$ (28c)

we obtain $a=27b^2{p}_c=\frac{27b}{8{\beta }_c}$ ,

so

$b=\frac{1}{8{\beta }_c{p}_c}$ , $a=\frac{27}{64{\beta }_c^2{p}_c}$ (28d)

with molecular parameters we obtain

$a=a_{W0}\epsilon {\sigma }^3$ (28e)

for Lennard-Jones potential $u_{LJ}\left(r,\sigma ,\epsilon \right)=4\epsilon \left({\left(\frac{\sigma }{r}\right)}^{12}-{\left(\frac{\sigma }{r}\right)}^6\right)$ ,

$a_{W0}=23.4$ [12] [13], (28f)

for dipole-dipole potential $u_{DD}\text{}\left(r,\sigma ,\theta \right)=4\epsilon \left({\Theta }_L\left(r,\sigma ,\Delta r\right)-\frac{{\Theta }_H\left(r,\sigma ,\Delta r\right)}{{\left(r/\sigma \right)}^3+\Delta r}\cos \left(\theta \right)\right)$ , $\Delta r=0.1\sigma$ ,

$a_{W0}=25.7$ [12] [13] (28g)

Parameter mixing rules

Parameter mixing rules determine the vdW parameters of the solution with component concentrations $x_i$ from the component parameters $a_i$ $b_i$ , in simplest form [14]

$b_1={\displaystyle \sum _i{x}_i}{b}_{1,i}$ , $a_1={\displaystyle \sum _{i,j}{x}_i}{x}_j{a}_{1,ij}$ , (29)

$a_{1,ij}=\sqrt{a_{1,i}{a}_{1,j}}$ geometric (GMA), $a_{1,ij}=\frac{2\sqrt{a_{1,i}{a}_{1,j}}+a_{1,i}+a_{1,j}}{4}$ expanded geometric (EGA), $a_{1,ij}=\left(a_{1,i}+a_{1,j}\right)/2$ simple arithmetic (SA).

Molecular mixing rules

Parameter mixing rules calculate the molecular parameters ${\epsilon }_{ij}$ (characteristic energy) and ${\sigma }_{ij}$ (effective hard-core radius) for the partial pressure $p_{ij}$ from the component parameters ${\epsilon }_i=a_i/{\sigma }_i^3$ and ${\sigma }_i={\left(3b_i/2\right)}^{1/3}$ , the resulting partial pressure is then [15]

$p_{ij}=\frac{1}{\left(v-2{\sigma }_{ij}^3/3\right)\beta }-\frac{{\epsilon }_{ij}{\sigma }_{ij}^3}{v^2}$ , and the total pressure is $p={\displaystyle \sum _{i,j}{x}_i{x}_j{p}_{ij}}$ (30a)

A widely used rule is Lorentz-Berthelot

${\sigma }_{ij}=\left({\sigma }_i+{\sigma }_j\right)/2$ , ${\epsilon }_{ij}=\sqrt{{\epsilon }_i{\epsilon }_j}$ (30b)

which generates the simple GMA parameter mixing rule above.

An improved rule is Halgren HHG rule [15]

${\sigma }_{ij}=\frac{{\sigma }_i^3+{\sigma }_j^3}{{\sigma }_i^2+{\sigma }_j^2}$ , ${\epsilon }_{ij}=\frac{4{\epsilon }_i{\epsilon }_j}{{\left(\sqrt{{\epsilon }_i}+\sqrt{{\epsilon }_j}\right)}^2}$ (30c)

A fit to experimental data gives the experimental Al-Matar rule [15]

${\sigma }_{ij}={\left(\frac{0.5640{\sigma }_i^6+0.9464{\sigma }_i^3{\sigma }_j^3+0.4896{\sigma }_j^6}{2}\right)}^{1/6}$ (30d)

${\epsilon }_{ij}=\left(\frac{0.0799{\epsilon }_i+1.9129\sqrt{{\epsilon }_i{\epsilon }_j}+0.0071{\epsilon }_j}{2}\right)\left(\frac{{\sigma }_i^3{\sigma }_j^3}{{\sigma }_{ij}^6}\right)$

PSRK model (Predictive Soave-Redlich-Kwong)

The PSRK model provides reliable predictions of VLE (vapor-liquid-equilibria) and gas solubilities [16].

Therefore, the PSRK model was implemented in the different process simulators and is well accepted as a predictive thermodynamic model for the synthesis and design of the different processes in the chemical, gas proc-essing, and petroleum industries. But also the group contribution equation of state PSRK shows a few weaknesses. Because the SRK (Soave-Redlich-Kwong) equation of state is used in PSRK, poor results are calculated for liquid densities of the pure compounds and the mixtures.

Better results are achieved with the improved Peng-Robinson eos.

Redlich-Kwong equation

The Redlich-Kwong equation is a real-gas equation and is formulated as (extended vdWaals) [17]

$p=\frac{RT}{\left(V_m-b\right)}-\frac{a}{\sqrt{T}{V}_m\left(V_m+b\right)}$ (31a)

where V_m = V/N_A is the molar volume, a, b are the generalized vdWaals constants, the constants a, b depend on critical values of the gas:

$a=0.42748\frac{R^2{T}_c^{2.5}}{p_c}$ , $b=0.08664\frac{R{T}_c}{p_c}$

The equation can be formulated in specific (per particle) volume v = V/N, $a_1=\frac{0.42748}{p_c{\beta }^{5/2}}$ , with specific generalized vdWaals parameters $b_1=\frac{b}{N_A}$ , $a_1=\frac{a}{N_A}\sqrt{k_B}$

$a_1=\frac{0.42748}{p_c{\beta }^{2.5}}$ , $b_1=0.08664\frac{1}{p_c{\beta }_c}$ (31b)

PSRK Mixing rule

$b={\displaystyle \sum _i{x}_i}{b}_i$ , $\frac{a}{bRT}={\displaystyle \sum _i{x}_i}\frac{a_{ii}}{b_{i}RT}+\frac{1}{A}\left(\frac{g^E}{RT}+{\displaystyle \sum _i{x}_i}\ln \frac{b}{b_i}\right)$ , $A=-0.64663$ (32a)

specific generalized variables

$b_1={\displaystyle \sum _i{x}_i}{b}_{1,i}$ , $\frac{a_1\beta }{b_1}={\displaystyle \sum _i{x}_i}\frac{a_{1,i}\beta }{b_{1,i}}+\frac{1}{A}\left(g^E\beta +{\displaystyle \sum _i{x}_i}\ln \frac{b_1}{b_{1,i}}\right)$ , (32b)

where $g_{}^E$ is the excess Gibbs energy.

Peng-Robinson equation

[18]

The Peng-Robinson equation is an improvement of the Redlich-Kwong equation in the form

$p=p\left(\beta ,v,p_c,{\beta }_c,\omega \right)=\frac{1}{\beta \left(v-b\right)}-\frac{\alpha a}{v\left(v+b\right)+b\left(v-b\right)}$ (33)

where

$a=\frac{0.45723}{p_c{\beta }_c^2}$ , $\alpha ={\left(1+\left(0.480+1.574\omega -0.176{\omega }^2\right)\left(1-\sqrt{\beta /{\beta }_c}\right)\right)}^2$ , $b=0.077796\frac{1}{p_c{\beta }_c}$ (33a)

ω is the acentric factor $\omega =-{\log }_{10}\left(\frac{p_{sat}\left(0.7T_c\right)}{p_c}\right)$ .

with critical parameters

${\beta }_c=\frac{b}{a}\frac{0.45723}{0.077796}=5.8773\frac{b}{a}$ , $p_c=\frac{0.077796}{b{\beta }_c}=\frac{0.077796}{5.8773}\frac{a}{b^2}=0.013236\frac{a}{b^2}$ (33b)

with vdWaals parameters

$\frac{1}{{\beta }_c}=k{T}_c=\frac{8a_W}{27b_W}$ , $v_c=3b_W$ , $p_c=\frac{a_W}{27b_W^2}$ (33c)

$b_W=\frac{2\pi {\sigma }^3}{3}=\frac{1}{8{\beta }_c{p}_c}=\frac{0.125}{0.077796}b=1.607b$ , $a_W=a_{W0}\epsilon {\sigma }^3=\frac{27}{64{\beta }_c^2{p}_c}=\frac{0.4219}{{\beta }_c^2{p}_c}$ ,

so

$a=\frac{0.45723}{p_c{\beta }_c{}^2}=1.083a_W$ , $b=0.077796\frac{1}{p_c{\beta }_c}=0.62223b_W$

Now we can reformulate the eos with

$\begin{array}{c}{\beta }_c=5.8773\frac{b_1}{a_1}=5.8773\frac{0.62223}{1.083}\frac{b_W}{a_W}=3.37676\frac{b_W}{a_W}\\ =\frac{3.37676}{a_{W0}\epsilon {\sigma }^3}\frac{2\pi {\sigma }^3}{3}=\frac{7.07227}{a_{W0}\epsilon },\end{array}$ (33d)

$\begin{array}{c}{p}_c=0.013236\frac{a}{b^2}=0.013236\frac{1.083a_W}{{0.62223}^2{b}_W^2}\\ =0.03702\frac{a_{W0}\epsilon {\sigma }^3}{{\sigma }^6{\left(2\pi /3\right)}^2}=0.00844a_{W0}\frac{\epsilon }{{\sigma }^3},\end{array}$

$p=p\left(\beta ,v,\omega ,\epsilon ,\sigma \right)=\frac{1}{\beta \left(v-b\right)}-\frac{\alpha a}{v\left(v+b\right)+b\left(v-b\right)}$

$a=1.083a_{W0}\epsilon {\sigma }^3$ , $b=0.62223b_W=0.62223\frac{2\pi {\sigma }^3}{3}=1.303{\sigma }^3$

$\alpha ={\left(1+\left(0.480+1.574\omega -0.176{\omega }^2\right)\left(1-\sqrt{0.141\frac{a}{b}\beta }\right)\right)}^2$

Values of ω [19]:

$\omega =-0.302$ vdWaals

$\omega =0.304$ acetone

$\omega =-0.644$ ethanol

$\omega =0$ argon

$\omega =0.353$ benzene

The material parameters here are the critical parameters $p_c,{\beta }_c$ , and the acentric factor ω.

Chen mixing rule [20]

$b_1={\displaystyle \sum _{i,j}{x}_i}{x}_j{b}_{1,ij}$ , $b_{1,ij}^{3/4}=\left(b_{1,i}^{3/4}+b_{1,j}^{3/4}\right)/2$ , $\frac{a_1}{b_1}={\displaystyle \sum _i{x}_i}\frac{a_{1,ii}}{b_{1,i}}+\frac{g_{res}^E}{A}$ , $A=-0.53087$ , (34)

where $g_{res}^E$ is the residual excess Gibbs energy.

Thermodynamic variables for Peng-Robinson

Thermodynamic variables for Peng-Robinson eos are [21]

pressure p

$p=p\left(\beta ,v,p_c,{\beta }_c,\omega \right)=\frac{1}{\beta \left(v-b\right)}-\frac{\alpha a}{v\left(v+b\right)+b\left(v-b\right)}$ (35a)

free energy F

$\begin{array}{c}F=-{\displaystyle \int p\left(v,T\right)}\text{d}v\\ =-\frac{1}{\beta }\log \left(v-b\right)-\alpha \left(\omega ,\beta \right)a\left(\mathrm{arctanh}\left(\frac{b+v}{\sqrt{2b}}\right)-\mathrm{arctanh}\left(\sqrt{2b}\right)\right)\end{array}$ (35b)

including a T-term, the complete expression is

$\begin{array}{c}F=-{\displaystyle \int p\left(v,T\right)}\text{d}v\\ =-\frac{1}{\beta }\left(1-\frac{3}{2}\log \beta +\log \left(v-b\right)\right)-\alpha a\left(\mathrm{arctanh}\left(\frac{b+v}{\sqrt{2b}}\right)-\mathrm{arctanh}\left(\sqrt{2b}\right)\right)\end{array}$ (35c)

partition function

$\begin{array}{c}Z=\exp \left(-\beta F\right)\\ ={\beta }^{3/2}\left(v-b\right)\exp \left(\alpha \left(\omega ,\beta \right)a\beta \left(\mathrm{arctanh}\left(\frac{b+v}{\sqrt{2b}}\right)-\mathrm{arctanh}\left(\sqrt{2b}\right)\right)\right)\end{array}$ (35d)

Gibbs energy G

$\begin{array}{c}G=F+pv\\ =-\frac{1}{\beta }\left(1-\frac{3}{2}\log \beta +\log \left(v-b\right)\right)\\ -\alpha a\left(\mathrm{arctanh}\left(\frac{b+v}{\sqrt{2b}}\right)-\mathrm{arctanh}\left(\sqrt{2b}\right)\right)+pv\end{array}$ (35e)

The material parameters here are the critical parameters $p_c,{\beta }_c$ , and the acentric factor ω.

The acentric factor ω is an independent third parameter alongside a, b.

The parameters are functions of the critical values: $\omega$ , $a=a\left(p_c,{\beta }_c\right)$ , $b=b\left(p_c,{\beta }_c\right)$ .

When formulated with reduced variables $T_r=\frac{T}{T_c}$ , $v_r=\frac{v}{v_c}$ , $p_r=\frac{p}{p_c}$ , the

Peng-Robinson equation in reduced variables depends only on material-specific ω, whereas a and b are material-independent

$p=p\left(\beta ,v,\omega \right)=\frac{1}{\beta \left(v-b\right)}-\frac{\alpha \left(\omega ,\beta \right)a}{v\left(v+b\right)+b\left(v-b\right)}$ (35f)

where

$a=0.45723$ , $b=0.077796$

$\alpha ={\left(1+\left(0.480+1.574\omega -0.176{\omega }^2\right)\left(1-\sqrt{\beta }\right)\right)}^2$

ω is the acentric factor $\omega =-{\log }_{10}\left(p_{sat}\left(1/0.7\right)\right)$ .

Below in Figure 2, the p(v) curves for different temperature values are shown for $\omega =0.2$ (carbon dioxide) [12] [13]

Figure 2. Peng-Robinson p(v) curves for different temperature values for $\omega =0.2$ (carbon dioxide).

Tait-Tammann equation of state

The Tait-Tammann equation [22] is adapted and intensely tested for water, it is an extension of the vdWaals equation in exponential form, for pressure p (bar) in dependence of $w=1/\rho$ (m³/kg) and temperature T (degree K).

$p=\left(p_0+B\right)\exp \left(-\frac{w-w_0}{C_0{w}_0}\right)-B$ , $p_0=p\left(w_0,T_0\right)$ normal pressure, $w=1/\rho$ (36a)

where $C_0$ is a material constant, for water $C_0=0.315$ .

$B=B_0+B_1\left(T-T_0\right)$ , with material constants $B_0,B_1$ , for water $B_0=2996\text{bar}$ , $B_1=7.555\text{bar}\cdot {\text{K}}^{-1}$

With the introduction of molecular mass m₀, it can be reformulated in standard variables (per particle) $p=P/N$ , $v=V/N$ , $\beta =1/k_{B}T$ in the form

$p=\left(p_0+b\right)\exp \left(-\frac{v-v_0}{C_0{v}_0}\right)-b$ , $p_0=p\left(v_0,{\beta }_0\right)$ , $b=B\frac{m_0}{1\text{kg}}$ , (36b)

Comparison of eos and mixing rules results

A comparison of eos’s for three mixing rules, namely, Geometric Mean Average (GMA), Expanded Geometric Average (EGA), and Simple Average (SA) for the ammonia-water binary system, based on absolute average deviation (AAD) in percent of the measured value is presented in [11].

Here, the following mixing rules are used for vdWaals parameter a

• geometric mean average (GMA)

$a_{ij}=\sqrt{a_i{a}_j}$ (37a)

• expanded geometric average (EGA)

$a_{ij}=\frac{2\sqrt{a_i{a}_j}+a_i+a_j}{4}$ (37b)

• the simple average (SA)

$a_{ij}=\frac{a_i+a_j}{2}$ (37c)

The below Table 6 gives an assessment of precision for different eos, the best relative error has the Peng-Robinson eos (about 10%).

Table 6. Absolute average deviation for GMA, EGA, and SA for ammonia – water binary system [14].

eos	mixing rule	% AAD
vdWaals	GMA	46.8
	EGA	46.7
	SA	59.7
Redlich-Kwong	GMA	20.6
	EGA	20.5
	SA	29.9
Peng-Robinson	GMA	9.9
	EGA	9.9
	SA	18.9

4.2. Extended Fluid-Gas Equations-of-State

Advanced Peng-Robinson (APR)

APR introduces an additional volume shift c in addition to the volume parameter b.

The modified (molar) Peng-Robinson eos has the form [23] [24]

$p=p\left(T,V_m\right)=\frac{RT}{\left(V_m-b\right)}-\frac{\alpha a}{V_m\left(V_m+b\right)+c\left(V_m-b\right)}$

where $c=\left(1-3\chi \right)\frac{R{T}_c}{p_c}$ , and V_m is the molar volume.

$\chi$ is the adjusted critical compressibility factor and for non-polar substances

$\chi =\left(0.329032-0.076788\omega +0.0211947{\omega }^2\right)$

APR modifies the α-factor in the ω-polynomial

$\alpha ={\left(1+m\left(\omega \right)\left(1-\sqrt{T_c/T}\right)\right)}^2$ ,

where

$m\left(\omega \right)=\left(0.452413+1.3098\omega -0.295937{\omega }^2\right)$

instead of original Peng-Robinson $m\left(\omega \right)=\left(0.480+1.574\omega -0.176{\omega }^2\right)$ .

Benedict-Webb-Rubin (BWR)

BWR is a polynomial-exponential eos, which in molar form reads [23] [25]

$p=p\left(T,V_m\right)=\frac{RT}{V_m}-\frac{B_{0}RT-A_0-C_0{T}^{-2}}{V_m^2}+\frac{\alpha a}{V_m^6}+c{T}^{-2}\frac{V_m{}^2+\gamma }{V_m^5}\exp \left(-\gamma /V_m^2\right)$

with 8 parameters $A_0,B_0,C_0,a,b,c,\alpha ,\gamma$ .

Because of the exponential term, it is no longer analytically solvable for V_m, but offers much higher precision.

Quartic Shah’s eos

The generalized quartic eos has the following molar form [26]-[28]

$p=p\left(T,V_m\right)=\frac{RT}{\left(V_m-k_0\beta \right)}+\frac{\beta k_{1}RT}{{\left(V_m-k_0\beta \right)}^2}-\frac{a{V}_m+k_0\beta c}{V_m\left(V_m+e\right)\left(V_m-k_0\beta \right)}$

with four parameters: β is the hard-sphere molar volume, $a\left(T,\omega \right)$ , $c\left(T,\omega \right)$ are T-dependent, $e\left(\omega \right)$ is T-independent (ω is the acentric factor).

This model introduces a more realistic repulsive term (first and second term above), making it applicable for polar substances.

Furthermore, the attractive α-term from Peng-Robinson is inverse-cubic (instead of inverse-quadratic) and more precise.

The eos is analytically solvable for v (Tartaglia’s solution), and yields four roots, from which one is physically not feasible (negative real part), the other three correspond to the three roots of Peng-Robinson, so the calculation of saturation curves and phase diagrams can be handled in the same way as for Peng-Robinson.

All considered, this quartic eos is considerably more precise than Peng-Robinson for pure substances, as well as for binary and ternary mixtures.

Comparison of cubic and quartic eos

The error of eos’s for pure substances is usually measured as relative deviation from measured value in partition function Z resp. enthalpy H, and in critical pressure p_c and volume V_c.

A typical relative error in enthalpy is ~30% for vdW, and ~10% for PR and APR ([25], chap. 6).

A typical relative error in critical variables is ~3% for PR and APR, and reduces to ~1% for quartic Shah-type eos [28].

4.3. Solid Equations-of-State

Mie-Grueneisen eos

The Mie-Grueneisen eos has the form [29]

$p-p_0=\frac{{\rho }_0{v}_s^2\left(\eta -1\right)\left(\eta -\frac{{\Gamma }_0}{2}\left(\eta -1\right)\right)}{{\left(\eta -s\left(\eta -1\right)\right)}^2}+{\Gamma }_{0}E$ , $\eta =\frac{\rho }{{\rho }_0}$ (38a)

where v_s is the bulk speed of sound, ρ₀ is the initial density (reference state), ρ is the current density, Γ₀ is Grueneisen’s gamma at the reference state, s = dv_s/dv_p is the Hugoniot coefficient, v_s is the shock wave velocity,v_p is the particle velocity, and E is the internal energy density.

The internal energy density e can be computed using

$E=\frac{1}{V_0}{\displaystyle \int C_v\text{d}T}\approx c_v\left(T-T_0\right)$

From the Dulong-Petit law follows $c_v=n\frac{n_{dof}}{2}{k}_B$ , $E\approx \frac{n_{dof}}{2}\frac{k_B\left(T-T_0\right)}{v}$ , where

n = particle density, $n_{dof}$ = number of degrees-of -freedom; for solids $n_{dof}\approx 6$ $c_v\approx 3n{k}_B$ , so $E=3k_B\left(T-T_0\right)=3\left(E_{th}-E_{th,0}\right)$ per particle.

Also, as speed of sound $v_s=\sqrt{\frac{Y}{\rho }}$ , Y = Young modulus, we obtain

$p-p_0=\frac{Y\left(\eta -1\right)\left(\eta -\frac{{\Gamma }_0}{2}\left(\eta -1\right)\right)}{\eta {\left(\eta -s\left(\eta -1\right)\right)}^2}+3{\Gamma }_0\frac{\left(E_{th}-E_{th,0}\right)}{v_0}$ , where the Young modulus is slightly dependent on particle density $\rho =1/v$ , $Y=y_1\rho -y_2$ per particle.

With these relations, the Mie-Grueneisen eos reads per particle

$p\left(\eta ,\beta \right)=p_0+\frac{\left(y_1\eta /v_0-y_2\right)\left(\eta -1\right)\left(\eta -\frac{{\Gamma }_0}{2}\left(\eta -1\right)\right)}{\eta {\left(\eta -s\left(\eta -1\right)\right)}^2}+3\frac{{\Gamma }_0}{v_0}\left(\frac{1}{\beta }-\frac{1}{{\beta }_0}\right)$ , where $\eta =\frac{v_0}{v}$ , (38b)

The equation is derived from the Mie ansatz $p-p_0={\Gamma }_0\left(E-E_0\right)$ .

In solids, as opposed to fluids/gasses, the molecules are located on a crystal lattice with a lattice constant a.

For the intermolecular potential $u=u\left(\overrightarrow{r},\sigma ,\epsilon \right)$ with repulsive (hardcore) radius $\sigma$ , characteristic (well depth) energy $\epsilon$ , and well-minimum radius $r_0$ , the solid’s lattice constant becomes $a=2r_0$ (fcc) resp. $a=2r_0/\sqrt{3}$ (bcc).

At the well-minimum, the potential has the form of a harmonic oscillator

$u_0\left(r\right)=\frac{1}{2}\frac{{\partial }^{2}u\left(r=r_0\right)}{\partial r^2}{\left(r-r_0\right)}^2$ ,

i.e. $Y\approx \frac{1}{r_0}\frac{{\partial }^{2}u\left(r=r_0\right)}{\partial r^2}$ from $\left(\frac{K}{A}=\frac{\partial u_0\left(r\right)}{r_0{}^2\partial r}\right)=\frac{1}{r_0}\frac{{\partial }^{2}u\left(r=r_0\right)}{\partial r^2}\frac{\left(r-r_0\right)}{r_0}$ (38c)

From the crystal lattice ansatz follows ${\Gamma }_0\approx 2$ and $s\approx 3/2$ , which is confirmed experimentally.

The Mie-Grueneisen eos is derived from the crystal lattice ansatz and from the Hugoniot equations for the conservation of mass, momentum, and energy [29] [30].

Below in Figure 3, the Mie-Grueneisen eos is shown [12] [13] for different temperatures for carbon dioxide, with $y_1=8.44$ , $y_2=17.4$ [19].

Figure 3. Mie-Grueneisen eos for different temperatures for carbon dioxide.

Thermodynamic variables for Mie-Grueneisen eos

Thermodynamic variables for Mie-Grueneisen eos are as follows [21]:

pressure

$p\left(v,\beta \right)=p_0+\frac{\left(y_1-y_{2}v\right)\left(v_0-v\right)\left(v_0-\frac{{\Gamma }_0}{2}\left(v_0-v\right)\right)}{v_0{\left(v_0-s\left(v_0-v\right)\right)}^2}+3\frac{{\Gamma }_0}{v_0}\left(\frac{1}{\beta }-\frac{1}{{\beta }_0}\right)$ (39a)

free energy F

$\begin{array}{c}F=-{\displaystyle \int p\left(v,\beta \right)}\text{d}v\\ =-\frac{1}{\beta }\left(1-\frac{3}{2}\log \beta \right)-v\left(p_0+3\frac{{\Gamma }_0}{v_0}\left(\frac{1}{\beta }-\frac{1}{{\beta }_0}\right)\right)-\frac{v_0\left(\left(s-1\right)v_0{y}_2-s{y}_1\right)\left(s-{\Gamma }_0\right)}{s^4\left(s-1\right)\left(s\left(v-v_0\right)+v_0\right)}\\ -\frac{{\Gamma }_0{y}_2{\left(s\left(v-v_0\right)+v_0\right)}^2}{2s^4{v}_0}+\frac{\left(s\left(v-v_0\right)+v_0\right)\left(s{\Gamma }_0\left(y_1+3v_0{y}_2\right)-s{v}_0{y}_2\left({\Gamma }_0+1\right)\right)}{s^4{v}_0}\\ -\frac{\left(s^2\left(v_0{y}_2-y_1\right)+2s{y}_1{\Gamma }_0+3v_0{y}_2{\Gamma }_0-2s{v}_0{y}_2\left({\Gamma }_0+1\right)\right)\log \left(s\left(v-v_0\right)+v_0\right)}{s^4}\end{array}$ (39b)

Gibbs energy G

$\begin{array}{c}G=F+pv\\ =-\frac{1}{\beta }\left(1-\frac{3}{2}\log \beta \right)-v\left(p_0+3\frac{{\Gamma }_0}{v_0}\left(\frac{1}{\beta }-\frac{1}{{\beta }_0}\right)\right)-\frac{v_0\left(\left(s-1\right)v_0{y}_2-s{y}_1\right)\left(s-{\Gamma }_0\right)}{s^4\left(s-1\right)\left(s\left(v-v_0\right)+v_0\right)}\\ -\frac{{\Gamma }_0{y}_2{\left(s\left(v-v_0\right)+v_0\right)}^2}{2s^4{v}_0}+\frac{\left(s\left(v-v_0\right)+v_0\right)\left(s{\Gamma }_0\left(y_1+3v_0{y}_2\right)-s{v}_0{y}_2\left({\Gamma }_0+1\right)\right)}{s^4{v}_0}\\ -\frac{\left(s^2\left(v_0{y}_2-y_1\right)+2s{y}_1{\Gamma }_0+3v_0{y}_2{\Gamma }_0-2s{v}_0{y}_2\left({\Gamma }_0+1\right)\right)\log \left(s\left(v-v_0\right)+v_0\right)}{s^4}\\ +v{p}_0+\frac{v\left(y_1-y_{2}v\right)\left(v_0-v\right)\left(v_0-\frac{{\Gamma }_0}{2}\left(v_0-v\right)\right)}{v_0{\left(v_0-s\left(v_0-v\right)\right)}^2}+3\frac{{\Gamma }_0}{v_0}v\left(\frac{1}{\beta }-\frac{1}{{\beta }_0}\right)\end{array}$ (39c)

$\begin{array}{c}G=F+pv\\ =-\frac{1}{\beta }\left(1-\frac{3}{2}\log \beta \right)-v\left(p_0+3\frac{{\Gamma }_0}{v_0}\left(\frac{1}{\beta }-\frac{1}{{\beta }_0}\right)\right)-\frac{v_0\left(\left(s-1\right)v_0{y}_2-s{y}_1\right)\left(s-{\Gamma }_0\right)}{s^3\left(s-1\right)\left(s\left(v-v_0\right)+v_0\right)}\\ -\frac{{\Gamma }_0{y}_2\left(s\left(v-v_0\right)+v_0\right)}{s^3}+\frac{y_1\left({\Gamma }_0-1\right)\log \left(sv\right)}{{\left(s-1\right)}^2}\\ +\frac{\left(s^3\left(y_1-v_0{y}_2\right)+2v_0{y}_2{\Gamma }_0+2s^2\left(-{\Gamma }_0{y}_1+v_0{y}_2\left({\Gamma }_0+1\right)\right)+s\left({\Gamma }_0{y}_1-v_0{y}_2\left(1+4{\Gamma }_0\right)\right)\right)\log \left(s\left(v-v_0\right)+v_0\right)}{s^3{\left(s-1\right)}^2}\\ +v{p}_0+\frac{\left(y_1-y_{2}v\right)\left(v_0-v\right)\left(v_0-\frac{{\Gamma }_0}{2}\left(v_0-v\right)\right)}{{\left(v_0-s\left(v_0-v\right)\right)}^2}+3\frac{{\Gamma }_0}{v_0}v\left(\frac{1}{\beta }-\frac{1}{{\beta }_0}\right)\end{array}$ (39d)

4.4. Fluid-Gas Transition

The vdWaals eos is [3]

$p=\frac{kT}{v-b}-\frac{a}{v^2}$ , $v=\frac{V}{N}$ (40a)

The critical temperature T_c in the liquid-gas transition results from

$\frac{\text{d}p}{\text{d}v}=0$ , $\frac{{\text{d}}^{2}p}{\text{d}{v}^2}=0$ , follows $p_c{\left(v-v_c\right)}^3=0$

$k{T}_c=\frac{8a}{27b}$ , $v_c=3b$ , $p_c=\frac{a}{27b^2}$

with reduced variables $T_r=\frac{T}{T_c}$ , $v_r=\frac{v}{v_c}$ , $p_r=\frac{p}{p_c}$ .

vdWaals equation becomes universal $p_r=\frac{8}{3}\frac{T_r}{v_r-1/3}-\frac{3}{v_r^2}$ (40b)

and also and the compressibility ratio is universal $\frac{p_c{v}_c}{k_B{T}_c}=\frac{3}{8}=0.375$ .

Below in Figure 4 are shown four isotherms of the universal vdWaals equation

in relative coordinates with the spinodal curve $\frac{\partial p}{\partial v}=0$ (black dash-dot curve) and the saturation curve (red dash-dot curve). The critical point lies at the turning point $\frac{{\partial }^{2}p}{\partial v^2}=0$ on the orange isotherm [3] [12] [13].

Figure 4. dWaals eos p(T,v).

The saturation curve (left wing = fluid, right wing = gas) left (low volume) wing ends at the triple point, its points are determined by Maxwell’s equal-area rule [3]

$G\left(v_1\left(p,T\right),T\right)=G\left(v_2\left(p,T\right),T\right)$ , $p=p\left(v_1,T\right)$ , $p=p\left(v_2,T\right)$

Saturation curve vdWaals

For vdWaals fluid-gas the saturation curve can be calculated in closed form from Maxwell’s equal-area rule (Maxwell-Gibbs equation) [31], using relative variables,

and the universal vdWaals form $p_{vdW}\left(v,\beta \right)=\frac{8}{3\beta }\frac{1}{v-1/3}-\frac{3}{v^2}$

$p=\frac{8}{3\beta }\frac{1}{v_f-1/3}-\frac{3}{v_f^2}$ , $p=\frac{8}{3\beta }\frac{1}{v_g-1/3}-\frac{3}{v_g^2}$

The Maxwell-Gibbs equation reads $G_f=G_g$ , $eqgr\left(v_f,v_b\right)=0$

$eqgr\left(v_f,v_g,\beta \right)=-\log \frac{v_g-1/3}{v_f-1/3}+\frac{9}{4}\beta \left(\frac{1}{v_f}-\frac{1}{v_g}\right)+\frac{v_g}{\left(v_g-1/3\right)}-\frac{v_f}{\left(v_f-1/3\right)}$ (41)

These are two equations for the variables p, T, v_f, v_g, from which v_f, v_g can be

eliminated, giving the saturation curve in the form $p_{sat}\left(E_{th}\right)$ , where $E_{th}=\frac{T}{T_c}$ is

the relative temperature, and p is the relative pressure relative to p_c.

The vdWaals equation is cubic in v, we insert the smallest and the largest of the three roots into v_f,v_g (second and third Cardano’s root) $v_f=v_{vdW,2}\left(p,\beta \right)$ , $v_g=v_{v_{vdW,3}}\left(p,\beta \right)$ , and obtain the condition function $eqgr\left(p,E_{th}\right)$ .

The condition function (real and imaginary part) has the form Figure 5 [12] [13], where in the real part (left) the edge of the grey area marks the zero-condition, and in the imaginary part (right) the “wall” is the boundary of the real-valued region in the upper left half.

Figure 5. Real and imaginary part of the Maxwell-Gibbs equation for vdWaals fluid-gas.

The saturation curve runs along the “wall” on the edge of the grey area, and ends at the end of the shaded part of the “wall” at about $T=0.5T_c$ , which is the triple point.

$p_{sat}\left(E_{th}\right)$ is calculated numerically [12] [13] in relative coordinates (Figure 6).

Figure 6. vdWaals saturation (fluid-gas) curve.

Analytic solution for saturation curve

Lekner [31] (Figure 7) found an analytic parametric solution for the vdWaals saturation curve, without use of the the cubic roots of the eos. The solution is analytic, but does not allow to determine the triple point.

Starting with the ansatz

$\rho =v^{*}/\left(v-v^{*}\right)$

we obtain the Maxwell-Gibbs equation in the form

$\log \left(\frac{{\rho }_f}{{\rho }_g}\right)=\frac{\left({\rho }_f-{\rho }_g\right)\left({\rho }_f+{\rho }_g+2\right)}{{\rho }_f+{\rho }_g+2{\rho }_f{\rho }_g}$ (42a)

which yields the solution

$f\left(y\right)=\frac{y\cosh y-\sinh y}{\sinh y\cosh y-y}$

${\rho }_f\left(y\right)=f\left(y\right)\exp \left(y\right)/2$ , ${\rho }_g\left(y\right)=f\left(y\right)\exp \left(-y\right)/2$

$v_f=v^{*}\left(1+1/{\rho }_f\right)$ , $v_g=v^{*}\left(1+1/{\rho }_g\right)$

$p=p^{*}\frac{v^{*}{}^2\left(v_f{v}_g-v^{*}\left(v_f+v_g\right)\right)}{v_f^2{v}_g^2}$ , $T=T^{*}\frac{v^{*}\left(v_f+v_g\right)\left(v_f-v^{*}\right)\left(v_g-v^{*}\right)}{v_f^2{v}_g^2}$ (42b)

As a result, one obtains the saturation curve in dependence of the parameter $0\le y\le \infty$ , critical point $y=0$ : $f\left(0\right)=1$ , ${\rho }_f\left(0\right)={\rho }_g\left(0\right)=1/2$ , $T\left(0\right)=8T^{*}/27$ , $p\left(0\right)=p^{*}/27$ , $v_f\left(0\right)=v_g\left(0\right)=3v^{*}$ , triple point: $y=\infty$ , $f\left(y\right)~2\left(y-1\right)\exp \left(-y\right)$ , ${\rho }_f\left(y\right)~\left(y-1\right)$ , ${\rho }_g\left(y\right)~2\left(y-1\right)\exp \left(-2y\right)$ , $v_f\left(\infty \right)=v^{*}$ , $v_g\left(\infty \right)\to \infty$ , $T\left(\infty \right)=0$ , $p\left(\infty \right)=0$ .

The saturation curve $\left(T\left(y\right),p\left(y\right)\right)$ parametric in y is shown below in the middle in relative coordinates [31], the outer curves are the two spinoidal curves (first and third cubic root of p).

Figure 7. Lekner’s solution for the vdWaals saturation curve.

Approximate saturation curve

For approximate solutions, the ansatz which $p_r=p\left(v_r,T_r,\varphi \right)$ has been suggested, where $\varphi$ is a substance-dependent dimensionless parameter, $\varphi =\frac{p_c{v}_c}{k_B{T}_c}$ ,

a better candidate is $\omega =-{\log }_{10}\left(p_r\left(T_r=0.7\right)\right)$ .

The approximate saturation curve is [22]

$\log \left(p_{rs}\right)=5.37\left(1-\frac{1}{T_r}\right)+\omega \left(\text{7}\text{.49}-11.18T_r{}^3+3.69T_r{}^6+17.93\log T_r\right)$ (43)

with $\omega =-{\log }_{10}\left(p_r\left(T_r=0.7\right)\right)$ .

The family of saturation curves, showing the vdW curve as a member (blue curve), is shown in Figure 8. The blue dots are calculated from Lekner’s solution. The orange dots are calculated from data in the ASME Steam Tables [32].

Figure 8. Saturation curves in dependence of ω, Lekner’s vdW saturation curve in blue.

Saturation curve Peng-Robinson

For Peng-Robinson fluid-gas the saturation curve can be calculated in closed form from Maxwell’s equal-area rule [21]

$a=0.45723$ , $b=0.077796$

$\alpha ={\left(1+\left(0.480+1.574\omega -0.176{\omega }^2\right)\left(1-\sqrt{\beta }\right)\right)}^2$

$p=p_{PR}\left(v_f,\beta \right)=\frac{1}{\beta \left(v_f-b\right)}-\frac{\alpha a}{v_f\left(v_f+b\right)+b\left(v_f-b\right)}$ (44a)

$p=p_{PR}\left(v_g,\beta \right)=\frac{1}{\beta \left(v_g-b\right)}-\frac{\alpha a}{v_g\left(v_g+b\right)+b\left(v_g-b\right)}$ (44b)

$G_f=G_g$ , $eqgr\left(v_f,v_b\right)=0$ (44c)

$\begin{array}{l}eqgr\left(v_f,v_b,\beta ,\omega \right)\\ =-\left(\log \frac{v_g-b}{v_f-b}\right)-\alpha \left(\omega ,\beta \right)\beta a\left(\mathrm{arctanh}\left(\frac{b+v_g}{\sqrt{2b}}\right)-\mathrm{arctanh}\left(\frac{b+v_f}{\sqrt{2b}}\right)\right)\\ +\left(\frac{v_g}{\left(v_g-b\right)}-\frac{v_f}{\left(v_f-b\right)}-a\beta \left(\frac{1}{v_g}-\frac{1}{v_f}\right)\right)\end{array}$ (44d)

These are two equations for the variables p, T, v_f, v_g, from which v_f,v_g can be

eliminated, giving the saturation curve in the form $p_{sat}\left(E_{th}\right)$ , where $E_{th}=\frac{T}{T_c}$ is

the relative temperature, and p is the relative pressure relative to p_c.

The Peng-Robinson equation is cubic in v, we insert the smallest and the largest of the three roots into v_f,v_g (second and third Cardano’s root): $v_f=v_{PR,2}\left(p,\beta \right)$ , $v_g=v_{PR,3}\left(p,\beta \right)$ , and obtain the condition function $eqgr\left(p,E_{th}\right)$ and the Maxwell-Gibbs equation $eqgr\left(p,E_{th}\right)=0$ .

Maxwell-Gibbs equation gives a solution $p_{sat}\left(E_{th}\right)$ , where $eqgr\left(p,E_{th}\right)$ is real ($\mathrm{Im}\left(eqgr\left(p,E_{th}\right)\right)=0$ ), because otherwise we have two equations (real part and imaginary part = zero) for two variables, which yields a point instead of a curve.

• Triple point

For sufficiently low value $E_{th}=E_{th,t}<E_{th,c}$ , $eqgr\left(p,E_{th}\right)$ becomes complex, and the saturation curve ends at triple point $E_{th,t}$ .

• Critical point

The upper end point of the saturation curve is the critical point $E_{th,c}$ , where

the pressure as a turning point $\frac{{\partial }^{2}p\left(E_{th},v\right)}{\partial v^2}=0$ .

Clausius-Clapeyron equation

For first order phase transition fluid-gas with s discontinuous we have with Gibbs free energy G and entropy S we have $G\left(V_f,T\right)=G\left(V_g,T\right)$ gives

$\frac{\text{d}p}{\text{d}T}=\frac{S_g-S_f}{V_g-V_f}$ , with latent heat $Q_{gf}=T\left(S_g-S_f\right)$ , we obtain $\frac{\text{d}p}{\text{d}T}=\frac{Q_{gf}}{T\left(V_g-V_f\right)}$ , and using $\left(V_g-V_f\right)\approx V_g$ , $p{V}_g=k_{B}T$ , gives $\frac{\text{d}p}{\text{d}T}=\frac{Q_{gf}p}{k_B{T}^2}$ ,

$p\left(T\right)=p_0\exp \left(-\frac{Q_{gf}}{k_{B}T}\right)$ , (45)

which is the Clausius-Clapeyron equation.

4.5. Solid-Fluid, Solid-Gas Transition

We calculate the phase diagram for solid-fluid-gaseous carbon dioxide using the Peng-Robinson eos for fluid-gas and the Mie-Grueneisen eos for solid-fluid transition. The substance in consideration is here carbon dioxide.

It has the phase diagram [19] Figure 9.

Figure 9. Phase diagram of carbon dioxide.

Saturation curve fluid-gas

We use the Peng-Robinson eos with $\omega =0.2$ for CO₂ [33]

$a=0.45723$ , $b=0.077796$

$\alpha ={\left(1+\left(0.480+1.574\omega -0.176{\omega }^2\right)\left(1-\sqrt{\beta }\right)\right)}^2$

with the Maxwell-Gibbs condition for the volume variables $v_f,v_g$ .

$p=p_{PR}\left(v_f,\beta \right)=\frac{1}{\beta \left(v_f-b\right)}-\frac{\alpha a}{v_f\left(v_f+b\right)+b\left(v_f-b\right)}$ , (46a)

solution $v_f=v_{PR,2}\left(p\right)$ second Cardano’s root

$p=p_{PR}\left(v_g,\beta \right)=\frac{1}{\beta \left(v_g-b\right)}-\frac{\alpha a}{v_g\left(v_g+b\right)+b\left(v_g-b\right)}$ , (46b)

solution $v_g=v_{PR,3}\left(p\right)$ third Cardano’s root

$G_f=G_g$ , $eqgr\left(v_f,v_g,\beta ,\omega \right)\equiv G_{PR}\left(v_g,\beta \right)-G_{PR}\left(v_f,\beta \right)=0$ , $E_{th}=1/\beta$ (46c)

$\begin{array}{l}eqgr\left(p,\beta ,\omega \right)\\ =-\left(\log \frac{v_g-b}{v_f-b}\right)-\alpha \left(\omega ,\beta \right)\beta a\left(\mathrm{arctanh}\left(\frac{b+v_g}{\sqrt{2b}}\right)-\mathrm{arctanh}\left(\frac{b+v_f}{\sqrt{2b}}\right)\right)\\ +\left(\frac{v_g}{\left(v_g-b\right)}-\frac{v_f}{\left(v_f-b\right)}-a\beta \left(\frac{1}{v_g}-\frac{1}{v_f}\right)\right)\end{array}$ (46d)

We obtain for the Maxwell-Gibbs equation $eqgr\left(\beta ,p\right)$ [34] in Figure 10, where in the real part (left) the edge of the grey area marks the zero-condition, and the imaginary part (right) is zero in the area under consideration.

Figure 10. Maxwell-Gibbs equation $eqgr\left(\beta ,p\right)$ for fluid-gas saturation curve.

$p_{sat}\left(E_{th}\right)$ is calculated numerically [12] [13] in relative coordinates in Figure 11.

Figure 11. Fluid-gas saturation curve Peng-Robinson eos.

Saturation curve solid-fluid (melting curve)

We use the Mie-Grueneisen eos for СO₂ for the volume variable v_s

$p_{MG}\left(v,\beta \right)=p_0+\frac{\left(y_1-y_{2}v\right)\left(v_0-v\right)\left(v_0-\frac{{\Gamma }_0}{2}\left(v_0-v\right)\right)}{v_0{\left(v_0-s\left(v_0-v\right)\right)}^2}+3{\Gamma }_0\left(\frac{1}{\beta }-\frac{1}{{\beta }_0}\right)$ (47a)

with the Young Modulus in relative coordinates $Y\left(v\right)=\left(y_1/v-y_2\right)$ , with $y_1=0.844$ , $y_2=1.74$ for CO₂ [19] and the Peng-Robinson eos

$p_{PR}\left(v,\beta \right)=\frac{1}{\beta \left(v-b\right)}-\frac{\alpha a}{v\left(v+b\right)+b\left(v-b\right)}$ (47b)

for the volume variable v_f,

$p=p_{MG}\left(v_s,\beta \right)$ , $p=p_{PR}\left(v_f,\beta \right)$ , (47c)

solution $v_f=v_{PR,3}\left(p,\beta \right)$ third Cardano’s root of Peng-Robinson eos, and solution $v_s=v_{MG,1}\left(p,\beta \right)$ first Cardano’s root of Mie-Grueneisen eos, and the Maxwell-Gibbs condition for the volume variables $v_s,v_f$

$G_f=G_s$ , $eqgmr\left(p,\beta \right)\equiv G_{PR}\left(v_f\left(p,\beta \right),\beta \right)-G_{MG}\left(v_s\left(p,\beta \right),\beta \right)=0$ , where $E_{th}=1/\beta$ (47d)

We obtain for the Maxwell-Gibbs equation $eqgmr\left(\beta ,p\right)=0$ [34] the real part (Figure 12(a)).

The melting curve lies at the steep left edge of the “boot”.

$p_{sat}\left(E_{th}\right)$ is calculated numerically [3] in relative coordinates in Figure 12(b).

As depicted in the measured phase diagram above, the solid-fluid curve rises steeply near the triple point to about the critical pressure p_c, and then bends off towards about 2p_c near the critical temperature T_c.

Saturation curve solid-gas (sublimation curve)

We use the Mie-Grueneisen eos for СO₂ for the volume variable v_s

$p_{MG}\left(v,\beta \right)=p_0+\frac{\left(y_1-y_{2}v\right)\left(v_0-v\right)\left(v_0-\frac{{\Gamma }_0}{2}\left(v_0-v\right)\right)}{v_0{\left(v_0-s\left(v_0-v\right)\right)}^2}+3{\Gamma }_0\left(\frac{1}{\beta }-\frac{1}{{\beta }_0}\right)$ (48a)

with the Young Modulus in relative coordinates $Y\left(v\right)=\left(y_1/v-y_2\right)$ , with $y_1=0.844$ , $y_2=1.74$ for CO₂ [19]

and the Peng-Robinson eos $p_{PR}\left(v,\beta \right)=\frac{1}{\beta \left(v-b\right)}-\frac{\alpha a}{v\left(v+b\right)+b\left(v-b\right)}$ (48b)

for the volume variable v_g, $p=p_{MG}\left(v_s,\beta \right)$ , $p=p_{PR}\left(v_g,\beta \right)$ , (48c)

solution $v_g=v_{PR,3}\left(p,\beta \right)$ third Cardano’s root of Peng-Robinson eos, and solution $v_s=v_{MG,3}\left(p,\beta \right)$ third Cardano’s root of Mie-Grueneisen eos, and the Maxwell-Gibbs condition for the volume variables $v_s,v_g$

$G_g=G_s$ , $eqgmg\left(p,\beta \right)\equiv G_{PR}\left(v_g\left(p,\beta \right),\beta \right)-G_{MG}\left(v_s\left(p,\beta \right),\beta \right)=0$ (48d)

(a) (b)

Figure 12. (a) Maxwell-Gibbs equation $eqgmr\left(\beta ,p\right)$ for solid-fluid saturation curve; (b) Solid-fluid saturation curve for Peng-Robinson & Mie-Grueneisen eos.

We obtain for the Maxwell-Gibbs equation $eqgmg\left(p,\beta \right)$ [12] [13] the real part in Figure 13.

Figure 13. Maxwell-Gibbs equation $eqgmg\left(p,\beta \right)$ for solid-gas sublimation curve.

The sublimation curve lies at the bottom of the “boot”.

$p_{sat}\left(E_{th}\right)$ is calculated numerically [12] [13] in relative coordinates in Figure 14.

Figure 14. Solid-gas sublimation curve for Peng-Robinson & Mie-Grueneisen eos.

5. Equations-of-State: Ansatz, Calculation and Results

5.1. Benzene

5.1.1. Material Data of Benzene

The important material data of benzene are as follows [19].

Benzene is an organic chemical compound with the molecular formula C₆H₆. The benzene molecule is composed of six carbon atoms joined in a planar hexagonal ring (Figure 15). The intermolecular potential is of Lennard-Jones type.

Figure 15. Benzene structure.

vdWaals parameters

$a_W=18.24{\text{L}}^2\cdot \text{bar}/{\text{mol}}^2=52.16\text{eV}\cdot {\text{A}}^3$

$b_W=0.1193\text{L}/\text{mol}=198{\text{A}}^3$

$p_0=a_W/b_W^2=1.33\text{meV}/{\text{A}}^3$

where $a=\frac{0.45723}{p_c{\beta }_c^2}=1.083a_W$ , $b=0.077796\frac{1}{p_c{\beta }_c}=0.62223b_W$ .

$\sigma \approx 3.70$ Å, $\epsilon /k_B\approx 400\text{K}$ .

ε σ parameters

$\sigma \left(\text{Benzene}\right)=3.7\dot{\text{A}}$, $\epsilon \left(\text{Benzene}\right)=k_{B}400\text{K}=\left(4/3\right)0.0259=0.0345\text{eV}$

crit. $T_c=562\text{K}$ (289˚C), $p_c=4.89\text{MPa}$ , ${\lambda }_c=7.47\text{A}$

${\beta }_c=1/0.0485{\text{eV}}^{-1}=20.6{\text{eV}}^{-1}$ , $p_c=30.5\mu \text{eV}\cdot {\text{A}}^{-3}=0.0229p_0$ , ${\lambda }_c=7.47\text{A}$

triple $T_t=278.5\text{K}$ , $p_t=4.83\text{kPa}$ , ${\lambda }_{tf}=5.31\text{A}$ , ${\lambda }_{ts}=5.33\text{A}$ , ${\lambda }_{tg}=75\text{A}$

${\beta }_t=1/0.0240{\text{eV}}^{-1}=41.6{\text{eV}}^{-1}$ , $p_t=30.1\times {10}^{-3}\mu \text{eV}\cdot {\text{A}}^{-3}$ , ${\lambda }_{tf}=5.31\text{A}$ ,

${\lambda }_{ts}=5.33\text{A}$ , ${\lambda }_{tg}=75\text{A}$

$Y=2\text{GPa}=12500\mu \text{eV}/{\text{A}}^3$

$\omega =0.212$

Phase diagram of benzene is shown in Figure 16.

Figure 16. Phase diagram of benzene.

5.1.2. Equation-of-State and Phase Diagram Benzene

The eos’s are shown in Figure 17(a), Figure 17(b) and Figure 18.

Figure 19(a), Figure 19(b), Figure 19(c) show the Maxwell-Gibbs equation and the diagram of the three saturation curves fluid-gas, solid-fluid, solid-gas.

Figure 20 shows the complete calculated phase diagram of benzene.

(a) (b)

Figure 17. (a) Peng-Robinson fluid-gas eos in relative coordinates; (b) Mie-Grueneisen solid eos.

Figure 18. The three branches (real part) of the volume function $v=v\left(p,\beta \right)$ Peng-Robinson eos.

(a)

(b)

(c)

Figure 19. (a) Maxwell-Gibbs eq. (real part) and diagram of the fluid-gas curve; (b) Maxwell-Gibbs eq. (real part) and diagram of the solid-fluid curve; (c) Maxwell-Gibbs eq. (real part) and diagram of the solid-gas curve.

Figure 20. Phase diagram of benzene.

5.2. Ethanol

5.2.1. Material Data of Ethanol

The important material data of ethanol are as follows [19].

Ethanol is an alcohol with the formula CH₃-CH₂-OH, and has a dipole intermolecular potential, mainly from H-H covalent bond (Figure 21).

Figure 21. Structure of ethanol.

vdWaals parameters

$a_W=12.56{\text{L}}^2\cdot \text{bar}/{\text{mol}}^2=35.9\text{eV}\cdot {\text{A}}^3$

$b_W=0.0871\text{L}/\text{mol}=144.6{\text{A}}^3$

where $a=\frac{0.45723}{p_c{\beta }_c^2}=1.083a_W$ , $b=0.077796\frac{1}{p_c{\beta }_c}=0.62223b_W$

$p_0=a_W/b_W^2=1.717\text{meV}/{\text{A}}^3$

ε σ parameters

${\epsilon }_h\left(\text{Eth}\right)=0.0542\text{eV}$ , $\sigma \left(\text{Eth}\right)=2.6\dot{\text{A}}$

crit. $T_c=513\text{K}$ , $p_c=6.25\text{MPa}$ , ${\lambda }_c=6.55\dot{\text{A}}$

crit. $\beta =1/0.0445{\text{eV}}^{-1}=22.5{\text{eV}}^{-1}$ , $p_c=39\mu \text{eV}\cdot {\text{A}}^{-3}=0.0227p_0$ , ${\lambda }_c=6.55\dot{\text{A}}$

triple $T_t=150\text{K}$ , $\beta =1/0.013{\text{eV}}^{-1}=77{\text{eV}}^{-1}$ , $p_t=4.3\times {10}^{-10}\text{MPa}=26.9\times {10}^{-10}\mu \text{eV}\cdot {\text{A}}^{-3}$ ,

${\lambda }_{tf}=4.6\dot{\text{A}}$, ${\lambda }_{ts}=4.4\dot{\text{A}}$

$Y={\rho }_0{v}_s^2=3.66\text{GPa}=22838\mu \text{eV}/{\text{A}}^3$

$\omega =0.62$

Phase diagram of ethanol is shown in Figure 22.

Figure 22. Phase diagram of ethanol.

5.2.2. Equation-of-State and Phase Diagram Ethanol

The eos’s are shown in Figure 23(a), Figure 23(b) and Figure 24.

(a) (b)

Figure 23. (a) Peng-Robinson fluid-gas eos in relative coordinates; (b) Mie-Grueneisen solid eos.

Figure 24. The three branches (real part) of the volume function $v=v\left(p,\beta \right)$ Peng-Robinson eos.

Figure 25(a), Figure 25(b), Figure 25(c) show the Maxwell-Gibbs equation and the diagram of the three saturation curves fluid-gas, solid-fluid, solid-gas.

Figure 26 shows the complete calculated phase diagram of ethanol.

(a)

(b)

(c)

Figure 25. (a) Maxwell-Gibbs eq. (real part) and the diagram of the fluid-gas curve; (b) Maxwell-Gibbs eq. (real part) and the diagram of the solid-fluid curve; (c) Maxwell-Gibbs eq. (real part) and the diagram of the solid-gas curve.

(a) (b)

Figure 26. (a) Phase diagram of ethanol; (b) Phase diagram of ethanol and benzene, relative to ethanol.

5.3. Argon

5.3.1. Material Data of Argon

The important material data of argon are as follows [19].

Argon belongs to the noble gasses, is practically chemically inactive, and its intermolecular potential is of Lennard-Jones type, it is even a typical Lennard-Jones substance.

vdWaals parameters

$a_W=1.355{\text{L}}^2\cdot \text{bar}/{\text{mol}}^2=3.87\text{eV}\cdot {\text{A}}^3$

$b_W=0.03201\text{L}/\text{mol}=53{\text{A}}^3$

$p_0=a_W/b_W^2=1.378\text{meV}/{\text{A}}^3$

where $a=\frac{0.45723}{p_c{\beta }_c^2}=1.083a_W$ , $b=0.077796\frac{1}{p_c{\beta }_c}=0.62223b_W$

ε σ parameters

$\sigma \left(\text{Ar}\right)=3.4\dot{\text{A}}$, $\epsilon \left(\text{Ar}\right)=0.0123\text{eV}$

crit. $T=150.8\text{K}$ , $p_c=4.83\text{MPa}=0.0218p_0$ , ${\lambda }_c=5.0\text{A}$

${\beta }_c=1/0.0130{\text{eV}}^{-1}=76.8{\text{eV}}^{-1}$ , $p_c=30.1\mu \text{eV}\cdot {\text{A}}^{-3}$

triple $T_t=83.7\text{K}$ , $p_t=0.068\text{MPa}$ , $v_{tf}=47.2{\text{A}}^3$

${\beta }_t=1/0.0072{\text{eV}}^{-1}=138.4{\text{eV}}^{-1}$ , $p_t=0.42\mu \text{eV}\cdot {\text{A}}^{-3}$ , ${\lambda }_{tf}=3.6\text{A}$ ,

${\lambda }_{tg}=25.8\text{A}$

$\omega =0.001$

$Y=1.6\text{GPa}=9.98\text{meV}/{\text{A}}^3$

Phase diagram of ethanol is shown in Figure 27.

Figure 27. Phase diagram of argon.

5.3.2. Equation-of-State and Phase Diagram Argon

Figure 28(a), Figure 28(b), Figure 28(c) show the Maxwell-Gibbs equation and the diagram of the three saturation curves fluid-gas, solid-fluid, solid-gas.

Figure 29 shows the complete calculated phase diagram of argon.

(a)

(b)

(c)

Figure 28. (a) Maxwell-Gibbs eq. (real part) and diagram of the fluid-gas curve; (b) Maxwell-Gibbs eq. (real part) and diagram of the solid-fluid curve; (c) Maxwell-Gibbs eq. (real part) and diagram of the solid-gas curve.

Figure 29. Phase diagram of argon.

5.4. Carbon Dioxide

5.4.1. Material Data of Carbon Dioxide

The important material data of carbon dioxide are as follows [19].

Carbon dioxide is a chemical compound with the chemical formula CO₂, made up of molecules that each have one carbon atom covalently double bonded to two oxygen atoms [33].

Carbon dioxide has the intermolecular potential of Lennard-Jones type, from O-O covalent binding (Figure 30).

Figure 30. Structure of carbon dioxide.

vdWaals parameters

$a_W=3.6{\text{L}}^2\cdot \text{bar}/{\text{mol}}^2=10.3\text{eV}\cdot {\text{A}}^3$

$b_W=0.0427\text{L}/\text{mol}=70.9{\text{A}}^3$

$p_0=a_W/b_W^2=2.05\text{meV}/{\text{A}}^3$

where $a=\frac{0.45723}{p_c{\beta }_c^2}=1.083a_W$ , $b=0.077796\frac{1}{p_c{\beta }_c}=0.62223b_W$

ε σ parameters

$\sigma \left({\text{CO}}_2\right)=3.03\dot{\text{A}}$, $\epsilon \left({\text{CO}}_2\right)=k_{B}125\text{K}=\left(1.25/3\right)0.0259=0.0108\text{eV}$

crit. $T_c=304.1\text{K}$ , $p_c=7.38\text{MPa}$ , $v_c=93.9{\text{cm}}^3/\text{mol}=156{\text{A}}^3$ , ${\rho }_c=0.0106\text{mol}/{\text{cm}}^3=0.00642{\text{A}}^{-3}$ , ${\lambda }_c=5.38\text{A}$

crit. ${\beta }_c=1/0.0262\text{eV}=38.1{\text{eV}}^{-1}$ , ${\lambda }_c=5.38\text{A}$ , $p_c=46\mu \text{eV}\cdot {\text{A}}^{-3}=0.0224p_0$

triple $T_{tr}=216.6\text{K}$ , $p_{tr}=0.52\text{MPa}$ , ${\rho }_{str}=1.562\text{g}/{\text{cm}}^3=0.0353\text{mol}/{\text{cm}}^3$ , ${\rho }_{ftr}=1.178\text{g}/{\text{cm}}^3=0.0266\text{mol}/{\text{cm}}^3$ , ${\rho }_{gtr}=0.00198\text{g}/{\text{cm}}^3=0.0000447\text{mol}/{\text{cm}}^3$

${\beta }_t=1/0.0186\text{eV}=53.5{\text{eV}}^{-1}$ , $p_c=3.2\mu \text{eV}\cdot {\text{A}}^{-3}$ , ${\lambda }_{tf}=4.0\text{A}$ , ${\lambda }_{tg}=33.6\text{A}$ ,

${\lambda }_{ts}=3.6\text{A}$ , $v_{tf}=64{\text{A}}^3$

$\omega =0.152$

$Y=y_1/\left(v/v_c\right)-y_2$

$y_1=8.439$ , $y_2=17.38$

$Y_{max}\approx 790\text{MPa}=4.93\text{meV}/{\text{A}}^3$

Phase diagram of carbon dioxide is shown in Figure 9.

5.4.2. Equation-of-State and Phase Diagram Carbon Dioxide

Figure 31(a), Figure 31(b), Figure 31(c) show the Maxwell-Gibbs equation and the diagram of the three saturation curves fluid-gas, solid-fluid, solid-gas.

Figure 32 shows the complete calculated phase diagram of carbon dioxide.

(a)

(b)

(c)

Figure 31. (a) Maxwell-Gibbs eq. (real part) and diagram of the fluid-gas curve; (b) Maxwell-Gibbs eq. (real part) and diagram of the solid-fluid curve; (c) Maxwell-Gibbs eq. (real part) and diagram of the solid-gas curve.

Figure 32. Phase diagram of carbon dioxide.

6. Solutions: Ansatz, Calculation and Results

6.1. Equation-of-State for Mixtures and Solutions

General ansatz

A simple and reliable ansatz for the eos of mixtures/solutions is the summing-up of partial pressures arising from pairwise interaction of the components

$p={\displaystyle \sum _{i,j}{x}_i{x}_j{p}_{ij}}$

where $x_i$ is the relative concentration of component i, and $p_{ij}=p\left(v,T,{\sigma }_{ij},{\epsilon }_{ij},\cdots \right)$ is the eos of the (i, j) components. This ansatz is theoretically well-founded and yields results in agreement with measurement within the precision margin of the Peng-Robinson eos (10%).

Pressure is in the molecular description momentum transfer, and addition of momenta is linear and independent of source, therefore this ansatz is generally valid in thermodynamics.

Furthermore, this ansatz has the following advantage for binary solutions. In this approach, the non-ideal character of a binary solution (the interaction between the two substances 1 and 2) is represented solely by the substance 12 (in the chosen example benzene-ethanol-50%), its parameters must be inserted in the model, on equal footing with the parameters of substance1 and substance 2. There is no need to introduce additional interaction relations in the model, in order to describe the non-ideal behavior of the solution.

Using Peng-Robinson eos with 3 parameters $\left(p_c,T_c,\omega \right)$ , a binary solution is completely described by 3 × 3 = 9 parameters.

The molecular parameters $\left({\sigma }_{ij},{\epsilon }_{ij},\cdots \right)$ are calculated from individual parameters $\left({\sigma }_i,{\epsilon }_i,\cdots \right)$ , $\left({\sigma }_j,{\epsilon }_j,\cdots \right)$ using molecular mixing rules (which introduces an additional error from mixing rules), or are measured directly (which is preferable of course).

Eos and mixing rules

For $\left(\sigma ,\epsilon \right)$ we use Lorentz-Berthelot

${\sigma }_{ij}=\left({\sigma }_i+{\sigma }_j\right)/2$ , ${\epsilon }_{ij}=\sqrt{{\epsilon }_i{\epsilon }_j}$ , (49a)

for volume v we apply the σ-rule and we obtain $v_{ij}={\left(\frac{{\left(v_i\right)}^{1/3}+{\left(v_j\right)}^{1/3}}{2}\right)}^3$ (49b)

• vdWaals eos

We formulate the dimensionless eos in the form in molecular parameters

$p_{ij}=\frac{1}{\left(v-2\pi {\sigma }_{ij}^3/3\right)\beta }-\frac{a_{W0}{\epsilon }_{ij}{\sigma }_{ij}^3}{v^2}$ (50a)

and obtain the total eos $p={\displaystyle \sum _{i,j}{x}_i{x}_j{p}_{ij}}$ (50b)

• Peng-Robinson eos parameter mixing

The dimensionless eos has the form in molecular parameters

$p=p\left(\beta ,v,\omega ,a,b\right)=\frac{1}{\beta \left(v-b\right)}-\frac{\alpha a}{v\left(v+b\right)+b\left(v-b\right)}$

where $a=\frac{0.45723}{p_c{\beta }_c{}^2}=1.083a_W$ , $b=0.077796\frac{1}{p_c{\beta }_c}=0.62223b_W$

$\alpha ={\left(1+\left(0.480+1.574\omega -0.176{\omega }^2\right)\left(1-\sqrt{0.141\frac{a}{b}\beta }\right)\right)}^2$

and${\beta }_c=\frac{b}{a}\frac{0.45723}{0.077796}=5.8773\frac{b}{a}$ , $p_c=\frac{0.077796}{b{\beta }_c}=\frac{0.077796}{5.8773}\frac{a}{b^2}=0.013236\frac{a}{b^2}$ , $v_c=3b_W=b3/0.622223=4.8214b$

and $b=0.62223\ast 2\pi {\sigma }^3/3=1.3032{\sigma }^3$ , $a=1.083a_{W0}\epsilon {\sigma }^3$ ,

For ω, we use its definition in the form $\omega =-{\log }_{10}\left(p_r\right)$ , $p_r=\frac{p_{sat}\left(0.7T_c\right)}{p_c}$ .

The energy ratio in p_r is constant $\frac{T}{T_c}=0.7$ , so p_r transforms like 1/v,

$p_{r,ij}={\left(\frac{2}{{\left(1/p_{r,i}\right)}^{1/3}+{\left(1/p_{r,j}\right)}^{1/3}}\right)}^3$ , (51a)

$\begin{array}{c}{\omega }_{ij}=-{\log }_{10}\left(p_{r,ij}\right)=3\left({\log }_{10}\left({\left(1/p_{r,i}\right)}^{1/3}+{\left(1/p_{r,j}\right)}^{1/3}\right)-{\log }_{10}2\right)\\ =3\left({\log }_{10}\left(\frac{{10}^{{\omega }_i/3}+{10}^{{\omega }_j/3}}{2}\right)\right)\end{array}$ (51b)

We apply the v-rule to b: $b_{ij}={\left(\frac{{\left(b_i\right)}^{1/3}+{\left(b_j\right)}^{1/3}}{2}\right)}^3$ (51c)

As for a, we see that $\tilde{a}=\frac{a}{b}$ is an energy, so we can apply the ε-rule

${\tilde{a}}_{ij}=\sqrt{{\tilde{a}}_i{\tilde{a}}_j}$ , $a_{ij}=\sqrt{{\tilde{a}}_i{\tilde{a}}_j}{b}_{ij}=\sqrt{{\tilde{a}}_i{\tilde{a}}_j}\frac{{\left({\left(b_i\right)}^{1/3}+{\left(b_j\right)}^{1/3}\right)}^3}{8}=\sqrt{\frac{a_i}{b_i}\frac{a_j}{b_j}}\frac{{\left({\left(b_i\right)}^{1/3}+{\left(b_j\right)}^{1/3}\right)}^3}{8}$ (51d)

• Mie-Grueneisen eos parameter mixing

$p\left(\eta ,\beta \right)=p_0+\frac{Y\left(\eta -1\right)\left(\eta -\frac{{\Gamma }_0}{2}\left(\eta -1\right)\right)}{\eta {\left(\eta -s\left(\eta -1\right)\right)}^2}+3\frac{{\Gamma }_0}{v_0}\left(\frac{1}{\beta }-\frac{1}{{\beta }_0}\right)$ per particle

$\frac{p\left(\eta ,\beta \right)}{p_0}=1+\frac{\frac{Y}{p_0}\left(\eta -1\right)\left(\eta -\frac{{\Gamma }_0}{2}\left(\eta -1\right)\right)}{\eta {\left(\eta -s\left(\eta -1\right)\right)}^2}-\frac{3{\Gamma }_0}{p_0{v}_0{\beta }_0}\left(1-\frac{{\beta }_0}{\beta }\right)$ , $\eta =\frac{v_0}{v}$

The reference point $\left(p_0,{\beta }_0,v_0\right)$ becomes the triple point $\left(p_t,{\beta }_t,v_t\right)$ in the phase diagram, so it must be calculated for the mixture/solution.

We apply the σ mixing rule for length to v

$v_{t,ij}={\left({\lambda }_{t,ij}\right)}^3={\left(\frac{{\lambda }_{t,i}+{\lambda }_{t,j}}{2}\right)}^3={\left(\frac{{\left(v_{t,i}\right)}^{1/3}+{\left(v_{t,j}\right)}^{1/3}}{2}\right)}^3$ (52a)

We apply the ε mixing rule for energy to β, T, pv and Yv

${\beta }_{t,ij}=\sqrt{{\beta }_{t,i}{\beta }_{t,j}}$ , $T_{t,ij}=\sqrt{T_{t,i}{T}_{t,j}}$ (52b)

$p_{t,ij}{v}_{t,ij}=\sqrt{p_{t,i}{v}_{t,i}{p}_{t,j}{v}_{t,j}}$ , so $p_{t,ij}=\frac{\sqrt{p_{t,i}{v}_{t,i}{p}_{t,j}{v}_{t,j}}}{v_{t,ij}}=\frac{\sqrt{p_{t,i}{v}_{t,i}{p}_{t,j}{v}_{t,j}}}{{\left(\frac{{\left(v_{t,i}\right)}^{1/3}+{\left(v_{t,j}\right)}^{1/3}}{2}\right)}^3}$ (52c)

$Y_{ij}=\frac{\sqrt{Y_i{Y}_j{v}_{t,i}{v}_{t,j}}}{v_{t,ij}}$ (52d)

6.2. Example Solutions

• benzene-ethanol solution [35]

Ethanol dissolves in benzene because of a balance of intermolecular interactions and favorable mixing entropy that allows the two liquids—polar protic ethanol and nonpolar aromatic benzene—to form a stable homogeneous solution at common concentrations.

• Dispersion forces dominate benzene

Benzene is nonpolar; its primary attractive forces are London (induced-dipole) dispersion forces between aromatic rings.

• Ethanol is amphiphilic

Ethanol (CH₃CH₂OH) has a small nonpolar ethyl group (${\text{CH}}_3{\text{CH}}_2^{-}$ ) and a polar hydroxyl group (–OH). The ethyl portion is compatible with benzene’s nonpolar environment; the hydroxyl can form weaker interactions with benzene than with water, but does not prohibit mixing.

• Dispersion interactions between ethanol’s alkyl portion and benzene are favorable

The ethyl group can engage in London dispersion with benzene, lowering the energy cost of breaking ethanol-ethanol and benzene-benzene contacts.

• Hydrogen bonding is not an absolute barrier

While ethanol forms hydrogen bonds with itself, those bonds are not so strong that they prevent disruption. In a benzene solution, some ethanol molecules retain hydrogen bonding clusters; many simply lose some H-bonds while gaining stabilizing dispersion contacts with benzene. The net energetic change can be small or slightly favorable.

• Entropy of mixing helps

Mixing increases configurational entropy, which favors solution formation even when enthalpic changes are modestly unfavorable. For small organic molecules like ethanol and benzene, the entropy term often offsets small positive enthalpy changes.

• Concentration dependence and limits

Ethanol is fully miscible with many organic solvents including benzene over a wide range, but its solubility in purely nonpolar solvents decreases as the solvent becomes less able to accommodate the hydroxyl group. Benzene can solvate moderate amounts of ethanol because ethanol’s nonpolar tail is sufficient to integrate into the benzene network; at very high ethanol fractions the system tends toward ethanol-like behavior with retained hydrogen bonding.

Measured values for the benzene-ethanol solution [16] are shown in Figure 33 for ${\gamma }^{\infty }$ activity coefficient at infinite dilution, $y_i$ component concentration in vapor, T_m melting temperature at normal pressure, H_E excess enthalpy.

• acetone-water

Acetone-water mixtures are miscible in all proportions and do not form a binary azeotrope. The boiling point depends on the mixture’s composition, ranging between the boiling points of pure acetone (56˚C - 56.5˚C) and pure water (100˚C).

Figure 33. Measured values for the benzene-ethanol solution.

Measured values for the acetone-water solution [4] are shown in Figure 34 for $y_i$ component concentration in vapor, ${\alpha }_i$ component separation factor, $\gamma$ activity coefficient, T_b boiling temperature at normal pressure. The double curves on the right refer to heated liquid, resp. cooled vapor, which show a hysteresis effect. The continuous line is the calculated value using the NRTL model.

Figure 34. Measured values for the acetone-water solution.

6.3. Saturation Curves of Solutions Benzene-Ethanol

$\left(E_{th},p,v\right)$ are given in relative coordinates (critical values $\left(E_{th,c},p_c,v_c\right)$ ) rel. benzene-ethanol-50 (i.e. 50%).

The concentration parameter x₀ is the relative benzene concentration $x_0=\left(1,0.75,0.5,0.4,0.25,0.1,0\right)$

For the eos we use here the ansatz $p={\displaystyle \sum _{i,j}{x}_i{x}_j{p}_{ij}}$ (53a)

This ansatz is theoretically well-founded and yields results in agreement with measurement within the precision margin of the Peng-Robinson eos (10%).

Specifically, for the Peng-Robinson eos we have

$\begin{array}{l}{p}_{PR}\left(E_{th},v,x_0\right)\\ =x_0^2{p}_{PR}\left(E_{th},v;E_{thc1},p_{c1},{\omega }_1\right)+{\left(1-x_0\right)}^2{p}_{PR}\left(E_{th},v;E_{thc2},p_{c2},{\omega }_2\right)\\ +2x_0\left(1-x_0\right)p_{PR}\left(E_{th},v;E_{thc12},p_{c12},{\omega }_{12}\right)\end{array}$ (53b)

where indices 1 = benzene, 2 = ethanol, 12 = benzene-ethanol-50.

Peng-Robinson eos is calculated in the individual relative coordinates, scaled from benzene-ethanol-50.

The parameters $E_{thc12},p_{c12},{\omega }_{12}$ of benzene-ethanol-50 are calculated according to the mixing rules (see above).

Specifically, for the Mie-Grueneisen eos we have

$\begin{array}{l}{p}_{MG}\left(E_{th},v,x_0\right)\\ =x_0^2{p}_{MG}\left(E_{th},v;Y_1,E_{th0},v_0,p_0\right)+{\left(1-x_0\right)}^2{p}_{MG}\left(E_{th},v;Y_2,E_{th0},v_0,p_0\right)\\ +2x_0\left(1-x_0\right)p_{MG}\left(E_{th},v;Y_{12},E_{th0},v_0,p_0\right)\end{array}$ (53c)

The parameter $Y_{12}$ of benzene-ethanol-50 is calculated according to the mixing rules.

The parameters $E_{th0},v_0,p_0$ are set to the triple point $E_{tht},v_t,p_t$ for the given x₀, calculated from the fluid-gas curve for this x₀.

The calculated Maxwell-Gibbs eq. (real part) and the diagram of the fluid-gas curve for different relative benzene concentration values of the benzene-ethanol solution are given in the following Figures 35(a)-(g).

The calculated Maxwell-Gibbs eq. (real part) and the diagram of the solid-fluid curve for different relative benzene concentration values of the benzene-ethanol solution are given in the following Figures 36(a)-(g).

(a)

(b)

(c)

(d)

(e)

(f)

(g)

Figure 35. (a) Fluid-gas curve of benzene-ethanol solution for x₀ = 1 (benzene); (b) Fluid-gas curve of benzene-ethanol solution for x₀ = 0.75; (c) Fluid-gas curve of benzene-ethanol solution for x₀ = 0.50; (d) Fluid-gas curve of benzene-ethanol solution for x₀ = 0.40; (e) Fluid-gas curve of benzene-ethanol solution for x₀ = 0.25; (f) Fluid-gas curve of benzene-ethanol solution for x₀ = 0.10; (g) Fluid-gas curve of benzene-ethanol solution for x₀ = 0 (ethanol).

(a)

(b)

(c)

(d)

(e)

(f)

(g)

Figure 36. (a) Solid-fluid curve of benzene-ethanol solution for x₀ = 1 (benzene); (b) Solid-fluid curve of benzene-ethanol solution for x₀ = 0.75; (c) Solid-fluid curve of benzene-ethanol solution for x₀ = 0.50; (d) Solid-fluid curve of benzene-ethanol solution for x₀ = 0.40; (e) Solid-fluid curve of benzene-ethanol solution for x₀ = 0.25; (f) Solid-fluid curve of benzene-ethanol solution for x₀ = 0.10; (g) Solid-fluid curve of benzene-ethanol solution for x₀ = 0 (ethanol).

The calculated Maxwell-Gibbs eq. (real part) and the diagram of the solid-gas curve for different relative benzene concentration values of the benzene-ethanol solution are given in the following Figures 37(a)-(g).

6.4. Phase Diagrams, Enthalpy, Characteristic Points of Solutions Benzene-Ethanol

Phase diagrams

The calculated phase diagrams for different relative benzene concentration values of the benzene-ethanol solution are shown in Figures 38(a)-(g).

(a)

(b)

(c)

(d)

(e)

(f)

(g)

Figure 37. (a) Solid-gas curve of benzene-ethanol solution for x₀ = 1 (benzene); (b) Solid-gas curve of benzene-ethanol solution for x₀ = 0.75; (c) Solid-gas curve of benzene-ethanol solution for x₀ = 0.50; (d) Solid-gas curve of benzene-ethanol solution for x₀ = 0.40; (e) Solid-gas curve of benzene-ethanol solution for x₀ = 0.25; (f) Solid-gas curve of benzene-ethanol solution for x₀ = 0.10; (g) Solid-gas curve of benzene-ethanol solution for x₀ = 0 (ethanol).

The combined phase diagram of all solutions is shown below Figure 39 in coordinates relative to benzene-ethanol-50%.

Each fluid-gas saturation curve descends steeply to the triple point of the corresponding solid-fluid curve, the triple point is the intersection of the solid-fluid and the solid-gas curve at the low pressure p_t < 0.01.

(a) (b)

(c) (d)

(e) (f)

(g)

Figure 38. (a) Phase diagram of benzene-ethanol solution for x₀ = 1 (benzene); (b) Phase diagram of benzene-ethanol solution for x₀ = 0.75; (c) Phase diagram of benzene-ethanol solution for x₀ = 0.50; (d) Phase diagram of benzene-ethanol solution for x₀ = 0.40; (e) Phase diagram of benzene-ethanol solution for x₀ = 0.25; (f) Phase diagram of benzene-ethanol solution for x₀ = 0.10; (g) Phase diagram of benzene-ethanol solution for x₀ = 0 (ethanol).

Figure 39. The combined phase diagram benzene-ethanol solution for different relative benzene concentrations.

Enthalpy

The excess enthalpy of the fluid-gas transition is the difference between the two enthalpies

$H_E\left(E_{th},p\right)=G\left(E_{th},v_g\left(E_{th},p\right),p\right)-G\left(E_{th},v_f\left(E_{th},p\right),p\right)$ , (54)

measured at normal pressure $H_{E,0}\left(E_{th}\right)=H_E\left(E_{th},p_0\right)$ , $p_0=1\text{bar}=0.1\text{MPa}$ .

Excess enthalpy $H_{E,0}\left(E_{th}\right)$ in kJ/mol vs. relative benzene concentration x₀ is shown below in Figure 40.

Figure 40. Calculated excess enthalpy $H_{E,0}\left(E_{th}\right)$ in kJ/mol vs. relative benzene concentration for the benzene-ethanol solution.

Measured excess enthalpy in J/mol is [16] roughly in agreement (Figure 41).

Figure 41. Measured excess enthalpy in J/mol for the benzene-ethanol vs. relative ethanol concentration.

The vaporization enthalpy is the enthalpy of the gas phase at normal pressure

$H_v\left(E_{th},p_0\right)=G\left(E_{th},v_g\left(E_{th},p_0\right),p_0\right)$ (55)

Vaporization enthalpy in kJ/mol is shown below in Figure 42.

Vaporization enthalpy depends only weakly on the temperature.

Our calculated values agree roughly with the calculation in [36], but deviate from the measured values for benzene at x₀ = 1 in [36] (Figure 43).

Figure 42. Calculated vaporization enthalpy in kJ/mol for the benzene-ethanol solution.

Figure 43. Measured and calculated [36] vaporization enthalpy in kJ/mol for the benzene-ethanol solution at T = 35˚C.

Triple point temperature

The calculated triple point temperature $E_{th,t}$ in K is shown in Figure 44.

Figure 44. Calculated triple point temperature in K for the benzene-ethanol solution.

The dependence on concentration is approximately linear.

In order to determine the melting point, we take the point on the steep descend part of the fluid-gas curve at normal pressure (Figure 45).

Figure 45. Calculated melting point for the benzene-ethanol solution vs. relative benzene concentration.

The resulting curve shows an accumulation near the benzene melting point at x₀ = 1.

Pure benzene melts (at normal pressure) at 5.53˚C (278.6 K), while pure ethanol melts at −114.14˚C (159 K).

The measured values of melting point [16] also show an accumulation effect near the melting point of benzene, and agree roughly with the calculated values (Figure 46).

Figure 46. Measured melting point for the benzene-ethanol solution vs. relative ethanol concentration.

7. Conclusions

In this paper, there are two important results:

• Exact solution of the Maxwell-Gibbs equation and calculation of phase diagrams

The Maxwell-Gibbs equation is the equality of enthalpy for both phases along the saturation curve, the other condition is continuity of pressure. The two equations yield the Maxwell-Gibbs condition $eqgr\left(p,E_{th}\right)=0$ after inserting an appropriate branch v_i(E_th,p) from the three algebraic solution (i = 1, 2, 3) of the cubic eos (vdWaals, Peng-Robinson and Mie-Grueneisen eos are all cubic in volume v). The condition $eqgr\left(p,E_{th}\right)=0$ is an algebraic-transcendent equation for (E_th, p),which yields the saturation curve in the form p(E_th), including the triple point, which is the algebraic branching point, where the three saturation curves meet.

This method is used to calculate the complete phase diagram for four selected substances (benzene, ethanol, argon, carbon dioxide), and the results are compared with measured data. The results agree with measured data within the accuracy of the Peng-Robinson eos (about 10%).

• Exact solution for the eos of binary solutions and calculation of their phase diagrams

We formulate an exact theoretical basis for binary solutions, based on the weighted sum of partial eos pressures, and including the 1-2-interaction of the components (i.e. non-ideal and irregular solutions).

Using this ansatz, we calculate the general eos in dependence on relative concentration x₀ of the first component. Furthermore, we calculate the eos for seven concentrations for the solution benzene-ethanol and compare the results with measurements. Again, the agreement is satisfactory and the deviation is within the accuracy of the Peng-Robinson eos (about 10%).

To achieve this, we introduce two novel methods.

• Exact algebraic solution for phase diagrams based on Peng-Robinson and Mie-Grueneisen equation-of-state

• A theoretically exact ansatz for mixture phase diagrams based on the weighted sum of partial pressures

• Limitations of the ansatz

The calculation ansatz for phase diagrams is based on the Maxwell-Gibbs condition for saturated curves, is theoretically exact.

The results for pure substances are limited in validity, however, by the precision of the underlying eos: the Peng-Robinson eos and the Mie-Grueneisen eos. In particular, Peng-Robinson is not well adapted to polar and ionic substances.

The precision will be increased by factor ~3, when the quartic Shah’s eos is used instead of Peng-Robinson, furthermore this eos is also applicable for polar fluids.

The calculation ansatz for solutions, pressure as the sum of weighted partial pressures, is also theoretically exact. For binary solutions, one needs reliable measurement data for parameters (v_c, T_c, ω, Y) for substance1, substance2, and substance12 = 50%-mixture. The data for the first two are usually available, but not for the latter. In consequence, the data for substance12 are calculated using mixing rules (Lorentz-Berthelot).

So, the precision is reduced by the error of the substance 12 parameters using mixing rules.

Conflicts of Interest

The author declares no conflicts of interest regarding the publication of this paper.

References

[1]	Jansen, H.J.F. (2008) Statistical Mechanics. Lecture Notes, Oregon State University.
[2]	Tuckerman, M.E. (2013) Statistical Mechanics. Lecture Notes, New York University.
[3]	Tong, D. (2012) Statistical Physics. Lecture Notes, University of Cambridge.
[4]	Rarey, J. and Gmehling, J. (2009) Factual Data Banks and Their Application to the Synthesis and Design of Chemical Processes and the Development and Test of Physical Property Estimation Methods. Pure and Applied Chemistry, 81.
[5]	Binney, J.J. (2002) The Theory of Critical Phenomena. Oxford University Press.
[6]	Helm, J. (2021) Physics Fundamentals. http://www.researchgate.net/publication/330737990
[7]	Goldenfeld, N. (2018) Lectures on Phase Transitions and the Renormalization Groups. CRC Press.[CrossRef]
[8]	deWith, G. (2013) Liquid-State Physical Chemistry. Wiley.
[9]	Subramanian, M. (2025) Solution Thermodynamics. Lecture Notes, Nadar College of Engineering.
[10]	Blandamer, M.J. and Reis, J.C. (2025) A Notebook for Topics in Thermodynamics of Solutions and Liquid Mixtures. Libre Texts Chemistry.
[11]	Thomsen, K. (2009) Electrolyte Solutions Thermodynamics. Lecture Notes, Denmark Technical University.
[12]	Helm, J. (2025) Statistical Mechanics of States of Matter. Mathematica code ThermoPhase.nb. http://www.researchgate.net/publication/393004134
[13]	Helm, J. (2026) Statistical Mechanics of States of Matter Version 2. Mathematica Code ThermoPhaseN.nb. https://www.researchgate.net/publication/393004134
[14]	Babalola, F.U., Akanji, I.O. and Oyegoke, T. (2021) Comparative Analysis of the Performance of Mixing Rules for Density Prediction of Simple Chemical Mixtures. Journal of Engineering Sciences, 8, F25-F31.[CrossRef]
[15]	Al‐Matar, A.K. and Rockstraw, D.A. (2004) A Generating Equation for Mixing Rules and Two New Mixing Rules for Interatomic Potential Energy Parameters. Journal of Computational Chemistry, 25, 660-668.[CrossRef] [PubMed]
[16]	Gmehling, J. (2003) Potential of Group Contribution Methods for the Prediction of Phase Equilibria and Excess Properties of Complex Mixtures. Pure and Applied Chemistry, 75, 875-888.[CrossRef]
[17]	Murdock, J.W. (1993) Fundamental Fluid Mechanics for the Practicing Engineer. CRC Press.
[18]	Tsonopoulos, C. and Heidman, J.L. (1985) From Redlich-Kwong to the Present. Fluid Phase Equilibria, 24, 1-23.[CrossRef]
[19]	Helm, J. (2025) Thermodynamics of Materials, Compendium of Data. http://www.researchgate.net/%20publication/400831845
[20]	Chen, J., Fischer, K. and Gmehling, J. (2002) Modification of PSRK Mixing Rules and Results for Vapor-Liquid Equilibria, Enthalpy of Mixing and Activity Coefficients at Infinite Dilution. Fluid Phase Equilibria, 200, 411-429.[CrossRef]
[21]	Helm, J. (2025) Statistical Mechanics of Phase Transitions. Journal of Modern Physics, 16, 1491-1557.[CrossRef]
[22]	Li, Y. (1967) Equation of State of Water and Sea Water. Journal of Geophysical Research, 72, 2665-2678.[CrossRef]
[23]	Loffredo, M. (2023) Optimization of Peng-Robinson and Redlich-Kwong-Soave Equations of State. Master’s Thesis, Politecnico di Torino.
[24]	Patel, N.C. and Teja, A.S. (1982) A New Cubic Equation of State for Fluids and Fluid Mixtures. Chemical Engineering Science, 37, 463-473.[CrossRef]
[25]	Poling, B.E., et al. (2001) The Properties of Gases and Liquids. 5th Edition, McGraw-Hill.
[26]	Lin, Y. (1994) Application of a Generalized Quartic Equation of State. Master’s Thesis, University of Tennessee.
[27]	Shah, V.M., Bienkowski, P.R. and Cochran, H.D. (1994) Generalized Quartic Equation of State for Pure Nonpolar Fluids. AIChE Journal, 40, 152-159.[CrossRef]
[28]	Li, D., Cao, J. and Yun, Z. (2011) A New Quartic Equation of State Based on a General Form and Its Application to Pure Fluids. Industrial & Engineering Chemistry Research, 50, 13576-13584.[CrossRef]
[29]	Burshtein, A.I. (2008) Introduction to Thermodynamics and Kinetic Theory of Matter. Wiley.
[30]	Lemons, D.S. and Lund, C.M. (1999) Thermodynamics of High Temperature, Mie-Gruneisen Solids. American Journal of Physics, 67, 1105-1108.[CrossRef]
[31]	Lekner, J. (1982) Parametric Solution of the Van Der Waals Liquid-Vapor Coexistence Curve. American Journal of Physics, 50, 161-163.[CrossRef]
[32]	Goodstein, D.L. (1985) States of Matter. Dover.
[33]	(2025) NIST Webbook. http://webbook.nist.gov/chemistry
[34]	Helm, J. (2025) Thermodynamics of Solutions. Mathematica code ThermoSolution.nb. http://www.researchgate.net/publication/400833974
[35]	Denisov, V., Denisov, V., Esina, Z., Esina, Z., Korchuganova, M. and Korchuganova, M. (2017) Phase Equilibrium in Systems Based on Aliphatic Bebzene Hydrocarbons. Science Evolution, 2, 33-39.[CrossRef]
[36]	Nguyen, M. (1989) Determination of Excess Enthalpy of Binary Mixtures. Master’s Thesis, Texas Tech University.