The multicomponent self-consistent Ornstein—Zernike application for CO2, N2, O2 shock Hugoniots simulation

A multicomponent equation of state with wide range of applicability is required to simulate shock waves in CxNyOz mixtures. This problem demands fine molecular interaction model due to competition between repulsion and attraction forces during shock compression process. A self-consistent Ornstein-Zernike application (SCOZA) based on distribution function integral equation theory can be used for it. The hypernetted-chain/soft core mean spherical approximation (HMSA) for SCOZA has been successfully applied to dense fluid systems with ambidextrous interactions. However, it was not designed to simulate mixtures, such as shock products of CxNyOz system. The convenient way to simulate multicomponent systems is the van der Waals one-fluid model (vdWlf). It has been shown, that vdWlf is not good enough for CO2 shock products at pressures higher, than 50 GPa. The multicomponent HMSA closure application based on partial version of the virial theorem has been offered in this paper. It is verified by molecular Monte-Carlo simulation at pressures up to 160 GPa with accuracy about 1–2%.


Introduction
Distribution function integral equation theory (DFIET) is one of the classical and well-known approaches in study of a dense fluid behavior. Amount of studies, based on multicomponent version of DFIET increases drastically during last decade. This growing of interest resides that, on the one hand, dense fluid mixtures is under the focus of different branches of science. On the other, DFIET provides fine simulation results with good computation efficiency, especially when compared with the Monte-Carlo (MC) methods.
A condition of self-consistency is required to simulate processes with wide range of thermodynamic parameters.
In this way, hypernetted-chain/soft core mean spherical approximation (HMSA) closure [1] with EXP-6 pair potentials was successfully applied for dense and overcritical fluids with ambidextrous molecular interactions. Unfortunately, this closure was designed for a single component. The van der Waals one-fluid model (vdW1f) model was proposed by Ree [2] to overcome this drawback.
It has been shown [3], that vdW1f model is not fine enough for mixtures of components, which pair potential parameters differ drastically, such as CO 2 shock products. This discrepancy becomes unacceptable at pressures higher, than 50 GPa [4]. The main aim of this paper is to present the application of HMSA closure for multicomponent systems, based on partial version of virial theorem for mixtures of shock products.

Theory
The full multicomponent equation of state (EOS), offered in this paper, is based on the Ornstein-Zernike (OZ) integral equation: where g ij represents radial distribution function of component i relatively to component j, h ij -total correlation function, c ij -pair correlation function, m-amount of components, ρ knumerical density of component k. An additional closure relation for pair and total correlation functions is needed [1] to compute radial distribution functions of system. An ordinary form of EOS as a relation of three thermodynamic parameters can be obtained via virial equation for previously received distribution function: where ρ-numerical density, An information about temperature dependence of distribution functions of the system is not represented in (1)-(2), so, it should be involved in a closure relation. Unfortunately, there are no theoretically proved closure relations for a general case of ambidextrous interactions. It can be constructed by interpolation between particular solutions in accordance with self-consistency conditions, given by theory. The HMSA [1] closure relation (4) is based on two asymptotic solutions analytically obtained for exponential repulsion at low radii and polynomial attraction at high radii: where φ R ij -repulsive potential, φ A ij -attractive potential, f ij -switching function. The switching functions interpolates solutions with unknown parameters α ij : to be determined from self-consistency conditions for each m 2 − m pair interactions of components.
The first m conditions of self-consistency can be obtained from the virial theorem. There are two ways to compute an isothermal compressibility: through the definition relation [5] and through density derivatives of the virial equation (3): Another one is a Maxwell equation: where E-internal energy, P -volume. Due to lack of self-consistency conditions, parameters were determined in this paper by numerical minimization of discrepancy between the definition and the finite difference virial derivative relations for isothermal compressibility: provided linear interpolation for unknown cross component parameters

Algorythm
To solve the equation system (1)- (4), a change of variables should be done: The HMSA (4) equation in this term takes a form: It should be noted, that HMSA produce a system of m algebraic equations in a multicomponent case. OZ equation, however, takes a form of matrix equation: A calculation of the convolution integral in (14) can be avoided by Fourier transform: Functions T ij could be retrieved from (15) as: A solution of algebraic equation system (13) and matrix equation (16) could be found by performing direct iterations in the case of dense fluids, far enough from a critical point: where upper index l denotes number of iteration. Up to 10 000 iterations for (17)- (18) and 30 iterations for (8) were performed in this work to obtain accuracy about 0.1% for isothermal compressibility. Unfortunately, the convergence of this method is not satisfactory for near-and supercritical fluids with compressibility factor βP ρ > 16. More sophisticated methods are required [6] for this case.

Simulation and results
All simulations in this paper were performed using EXP-6 pair potentials: where ǫ ij -potential well depth, α ij -stiffness of exponential repulsion, r max ij -radius of potential maximum, r min ij -radius of potential minimum (figure 1). Internal cutoff radius r max ij is not a free parameter. It should be determined from ǫ ij , α ij and r min ij . External cutoff radius was 20.48 times higher than r min ij in this paper. The potential parameters (table 1) was found in [7,8] by solving the inverse problem of shock waves and isothermal compression simulations through method of characteristic function extremum (MCFE) with Kang-Lee-Ree-Ree thermodynamic perturbation theory (KLRR-T) and vdW1f model. Chemical composition in view of dissociation and recombination was found by the same way via IVTANTHERMO [9] thermodynamic potentials available at http://www.chem.msu.su/rus/handbook/ivtan/. An independent realization of MCFE is required to validate the multicomponent HMSA closure application (MHMSA) for shock products simulation; so, it is not a goal of this paper. But perspectives of it can be assessed by joint verification of MHMSA, HMSA+vdW1f and KLRR-T+vdW1f via molecular Monte-Carlo (MC) simulation with the same composition. The MC simulation was performed by MCCCS Towhee software [10][11][12] avaliable at http://towhee.sourceforge.net/. Each run for system of 1000 molecules included 500 000 configurations.
The simulation of shock wave experiments [13][14][15]. was performed at first. During it, N 2 at T = 77 K, ρ = 0.808 g/cm 3 , E = −2.842 kcal/mole was shocked up to 90 GPa. No significant discrepancy was found for any theory (table 2) due to close values of pair interaction potential parameters for the dissociation products N 2 and N. Maximal deviation of pressure was about 1% for MHMSA and 4% for KLRR-T+vdW1f relatively to MC results. MCFE with KLRR-T+vdW1f provides accuracy about 7% relatively to experimental data. The same result (table 3) was obtained for O 2 experiment [13]. Initial state was at T = 77 K, ρ = 1.202 g/cm 3 , E = −1.413 kcal/mole. Maximal deviation of pressure was about 1% for MHMSA and 7% for KLRR-T+vdW1f relatively to MC results. MCFE with KLRR-T+vdW1f provides accuracy about 8% relatively to experimental data. One may suggest more interesting results for liquid CO 2 because of rich set of pair potential parameters of products. The initial state was at T = 218 K, ρ = 1.173 g/cm 3 , E = −98.486 kcal/mole [16]. Two simulations were performed. The chemical composition of the first (table 4) one was represented as binary mixture of CO 2 and O. All theoretical and experimental data are in a good agreement up to 50 GPa (figure 2). The HMSA+vdW1f and KLRR-T+vdW1f schemes coincide with each other but underestimate the pressure relatively to MC at higher pressures. High accuracy was found for MHMSA model with the same closure relation as HMSA. Therefore, the reason of previous lack is in vdW1f model. The accuracy of MCFE with KLRR-T+vdW1f can be analyzed via comparison of density-temperature dependency (figure 3) with first-principle models [17][18][19].      The structure of multicomponent fluid mixture can be analyzed by visualization of radial distribution functions (figure 4). One may see discrepancies between amplitudes as well as radii of maxima. Therefore, the disadvantage of vdW1f one fluid model is not surprising.

Conclusions
Multicomponent equations of state with wide range of applicability have been verified by MC for C x N y O z shock products mixture simulation. The drawback of vdW1f model has been confirmed at pressures, higher than 50 GPa. The new application of multicomponent HMSA closure for self-consistent OZ application has been offered. The multicomponent EOS for CO 2 with accuracy about 1-2% at pressures up to 160 GPa was developed.