Molecular dynamics simulation of multi-constituent gas diffusion properties in Venus atmosphere

The purpose of studying the Venus’s atmosphere is to model and simulate Venus’s environment. One of the key parameters of the Venus’s atmosphere is diffusion coefficient. Experimental measurements of diffusion coefficients are particularly difficult under Venus’s environmental conditions. Molecular dynamics have become an important method for studying the properties and dynamics of microscopic systems. In this paper, the equilibrium molecular dynamics (EMD) simulations are used to calculate the interdiffusion coefficients of carbon dioxide (CO 2 ) and nitrogen (N 2 ) at room temperature and pressure in combination with Darken's equation. And the results are compared with experimental values and empirical equations to verify the rationality of the calculation method and the accuracy of the results. The interdiffusion coefficients of trace gases on the surface of Venus for the CO 2 system in different states and the CO 2 -N 2 interdiffusion coefficients with altitude in the Venus environment are given. The results show that the diffusion coefficients of the gases on the surface of Venus are two orders of magnitude smaller than those in the Earth's atmosphere and molecular dynamics simulations can well predict the diffusion characteristics of the Venus’s atmosphere and support the simulation of the Venus’s surface environment and the Venus’s atmosphere model.


Introduction
Venus is the closest planet to the Earth, not only in distance but also in size and mass 1 . As the sister planet of the Earth, the detection of Venus is of great significance for the study of the origin and future evolution of the Earth. However, the environmental conditions on the Venus's surface are much harsher than Earth 2 , with high temperature, high pressure and acidic gases being the main characteristics of the Venus surface environment 1 . The study of the Venus's atmosphere is of great importance for the simulation of the Venus's environment and the modeling of the Venus's atmosphere. Venus landing probe are one of the most important ways to understand the surface environment of Venus. According to the Venus International Reference Atmosphere 3 (VIRA) and the deep atmospheric temperature obtained by the probe landing, the Venus's atmosphere has a maximum temperature of 470°C and a maximum pressure of 9 MPa, and the gas environment is composed of a molar fraction of 96.5% CO2 ,3.5% N2 and other trace gases such as SO2.
The atmosphere of Venus has been studied by many researchers. Kremic et al. 4 carried out a simulation study of the Venus surface gas environment with nine gases: CO2, N2, SO2, H2O, OCS, CO, H2S, HCl and HF. Sébastien et al. 5,6 carried out a verification experiment of the gradient of concentration in the atmosphere of Venus, and the diffusion of CO2 and N2 to a uniform distribution of concentration when the container was filled sequentially took about 100 h. Ando et al. 7 carried out the thermal structure of the Venusian atmosphere from the sub-cloud. According to the abovementioned studies on the Venus atmosphere, the gas diffusion coefficient is one of the key parameters in the study and environment simulation of the Venus atmosphere, which not only enables the Venus atmosphere environment to be modeled, but also guides the Venus atmosphere environment simulation test procedure, shortens the overall test cycle and improves the reliability of the system. Methods for calculating diffusion coefficients can be divided into empirical equations and experimental methods. Empirical equations for binary diffusion coefficients, which are generally a function of temperature, pressure, and the molar mass of the components, are typical: The Hirrschfelder-Bird-Spotz equation, Chen-Othmer equation, Fuller-Schettler and Giddings (FSG) equation 8 .
The Hirrschfelder-Bird-Spotz equation: Among them , , for absolute temperature, gas pressure in atm and molar mass. Ω is the collision integral, which depends on temperature and the interaction energy of molecular collisions, and . σ is the collision diameter between the gas molecules A and B. Most such values have been obtained by viscosity measurement.
Due to the complexity of the HBS equation, Chen and Othmer use the critical temperature with critical volume to provide an approximation, the Chen and Othmer equation: Fuller-Schettler and Giddings obtained an equation in which atomic and structural volume increments and other parameters were obtained by a least-squares fit to over 300 measurements. The Fuller-Schettler and Giddings equation: ∑ is determined by summing the atomic contributions shown in Table 1. The method of empirical equation varies considerably in accuracy for different binary systems due to limited experimental data and different fitting intervals 8 . However, for the empirical correlation method, the diffusion coefficient is exponentially related to the temperature when the pressure and the composition of the binary gas system are determined and does not reflect the physical changes that occur in certain states, such as the critical point state. This is the reason why more research is conducted into unconventional states and molecular dynamics simulations are used for calculations. In addition to empirical equation, diffusion coefficients can also be obtained by experiment. However, due to the high ambient temperature and pressure on the surface of Venus, the CO2-N2 gas mixture is in a supercritical state 6 , and it is difficult to measure the diffusion coefficient using physical tests at this temperature and pressure.
The molecular dynamics method of calculating diffusion coefficients is well generalized. Many scholars have calculated the diffusion coefficients of gases under supercritical conditions. Ge et al. 9 calculated the diffusion coefficients of krypton and argon under supercritical conditions and analyzed the difference between Fick and M-S diffusion coefficients.
Molecular dynamics simulation has always been uniquely advantageous for the validation of diffusion coefficient calculations, but the accuracy of molecular dynamics calculations is affected by the fact that there are many different molecular models of gases, such as water molecules, with different application conditions and scenarios. Keffer et al. 10 summarized four different methods of calculating the interdiffusion coefficient. In this paper, the equilibrium molecular dynamics (EMD) combined with the Darken's equation is used to calculate the diffusion coefficients of the Venus's atmosphere, because of its simplicity and accuracy 10 .

Comparison with experiment data at room temperature and pressure
The mean squared displacement (MSD) of N2 at different temperatures at one atmosphere is shown in Figure 1 as a function of time. It can be observed from the diagram that the slope of the MSD gradually increases as the temperature increases from 313.7K to 365.1K. The measurement of the diffusion coefficient experimental values based on a gas chromatographic method for measuring the interdiffusion coefficient in a gas atmosphere with a CO2 concentration of more than 99.8%, at different temperatures and under 1 atm. According to Darken's equation, at a CO2 gas concentration of more than 99.8%, the diffusion coefficient approximates the product of the N2 self-diffusion coefficient and the activity correlation coefficient in this gas environment. To reduce the large number of calculations required by introducing a large number of CO2 molecules in order to guarantee the number of N2 molecules, and small variation of self-diffusion coefficient between 10% and 1% N2 concentration makes it possible to use the self-diffusion coefficient at 10% N2 concentration instead of the N2 selfdiffusion coefficient of the above environment. The results are shown in Table 2 and Figure 2. The diffusion coefficient based on molecular dynamics calculations deviates from the experimental value by around 5%, and the correlation equation gives a diffusion coefficient that deviates from the experimental value by more than 10%.  Both the MD simulations and the empirical equation calculate diffusion coefficients were increasing with temperature. A comparison between the MD simulations and the empirical equation shows that the empirical correlation formula has a relatively high calculation error, with the Chen-Othmer equation having a lower calculation error of 10%, which is better than the MD simulations (which do not consider the molecular structure as well as the Coulomb potential).
By comparing the two methods of MD calculations, the error in calculating the diffusion coefficient is around 5% when the molecular structure as well as the Coulomb potential is taken into account, whereas it is 18.5% when the above conditions are not taken into account. This showing that the model that considers the molecular structure and the Coulomb potential in calculating the diffusion coefficient has the advantage of being more accurate.
When taking a molecular dynamic approach to calculating diffusion coefficients, considering more elaborate models (including bond-angle interactions) may give more accurate results, but this undoubtedly increases the cost of the calculation, so the method and model need to be chosen considering the conditions of the problem involved and the requirement for accuracy in calculating the diffusion coefficients.

Calculation results at room temperature and high pressure
After preliminary verification of the reasonableness of the calculation procedure, the mutual diffusion coefficients of trace gases and carbon dioxide gases are calculated for each trace gas at room temperature and high pressure (298 K, 3 MPa) and at the surface of Venus (773 K, 9 MPa). And the other parameters are set up in the same way as when molecular structure and Coulombic potential are taken into account (OCS molecules are not taken into account).
In the simulation of the Venus gas environment, it is necessary to pre-mix the gas components at room temperature and pressure in order to accurately measure whether the concentration of each gas component meets the requirements. The diffusion coefficient of CO2-N2 at room temperature and pressure (3 MPa, 298 K) is calculated using Darken's equation in combination with the MD method, as shown in Figure 3. It is clear from the calculations that the diffusion coefficient at room temperature and pressure is two orders of magnitude smaller than that at room temperature and high pressure, and that both the self-diffusion coefficient and the CO2-N2 interdiffusion obtained under these conditions increase with increasing N2 molar fraction. The necessary condition for a small variation in the coefficient of mutual diffusion to be satisfied is that the two gases in a binary gas system each have approximately the same coefficient of self-diffusion when they are each trace gases. However, it cannot be proven that these conditions are met in all binary gas systems under all conditions.
Calculate the diffusion coefficient of each gas at room temperature and high pressure as shown in Table 3. Table 3 Gas diffusion coefficient at room temperature and high pressure (298 K, 3 MPa) ( Figure 4 the profile of diffusion coefficients The predicted CO2 and N2 interdiffusion coefficients are calculated by the Venus Atmosphere Model, as shown in Figure 4. The diffusion coefficients exponential increase with the altitude of the Venus's atmospheric environment. The variation of the diffusion coefficients with height refines the model of the Venus's atmospheric model.

Discussion
Using equilibrium molecular dynamics, the Einstein's relation was used to calculate the binary gas diffusion coefficient of CO2-N2 at room temperature and pressure in conjunction with the Darken equation, considering the Coulomb potential as well as the molecular structure, and compared with the experimental and empirical value. The self-diffusion coefficients of other trace gases on the surface of Venus under CO2 were simulated in both two states. The simulations show: 1) The results of diffusion coefficient calculations based on molecular dynamics are in good agreement with the experimental values and considering the Coulomb potential as well as the molecular structure, the results are more accurate with an error of around 5% with experimental values.
2) The self-diffusion coefficients of CO2-N2 at room temperature and high pressure with different concentration ratios were calculated statistically, and the variation of the coefficients of mutual diffusion with nitrogen concentration was small in this state, the difference between the maximum and minimum values is in the order of 5%, so that the effect of concentration on the diffusion coefficient can be ignored for simplified analysis.
3) The self-diffusion coefficients of the gas components on the surface of Venus under room temperature and high pressure (298K, 3MPa), high temperature and high pressure (773K, 9MPa) with carbon dioxide were calculated, which provides a basis for analyzing the mass transfer characteristics of the gases on the surface of Venus and diffusion of gases released from Venus surface materials.
4) The interdiffusion coefficients of carbon dioxide and nitrogen in the Venus's bottom atmosphere shows that interdiffusion coefficient increases with the height, the maximum value is ten times the minimum value. The calculation improves the model of the Venus's bottom atmosphere and provides a basis for the analysis of the Venus's atmospheric environment.

Methods
The simulations in this paper are based on the DREIDING force-field model. 11 is a superposition of valence (or bonded) interactions. is non-bonding energy (depending only on the interatomic distance), which is determined by the van der Waals potential ,Coulomb potential , hydrogen bonding potential ℎ . Hydrogen bonding has been ignored to simplify the model in this work.
Treating gas molecules as rigid balls, the force field model only considers .
where is the depth of the potential energy traps, is the distance between atoms at zero interaction potential energy, is the interatomic distance, is the amount of electrical charge, is the energy conversion parameter, 0 is the dielectric constant. The atomic structure and potential energy parameters for the main gas components of the Venus surface environment are listed in Table 5 and Table 6. The potential energy model for all triatomic molecules in Table 5 is the L-J12-6 potential energy model, as shown in eq (5).  Table 6 shows the two-atom molecular gas model, LJ-12-6 for N2 13 and HCl 16 and LJ-9-6 for CO gas molecules from. 12 The gas model of HF is the Morse model 17 . Subscripts m in N2, HF, HCl represent imaginary atoms, which represent the location of charge aggregation and are calculated taking only the Coulombic force into account. For the atomic interaction potentials required in the calculations, the joint Lorentz-Berthelot rule is adopted. 18 Diffusion characterizes the migration of components in space, and for gases in general, Fick's law is often used to describe the diffusion process. When Fick's law describes a binary gas system, the diffusion coefficient is called the interdiffusion coefficient and the driving force for the diffusion process is generally considered to be a concentration gradient. When only a single component gas is considered and there are no diffusion driving forces, it is called self-diffusion. Self-diffusion characterizes the ability of a substance to migrate by mass at a given temperature and pressure, and when the definition of self-diffusion is extended to multi-gas mixtures, it can be used in combination with Darken's equation to describe the interdiffusion coefficient of a gas mixture.

Statistical equation
Self-diffusion coefficients can be calculated in equilibrium molecular dynamics using either the Einstein relation or the Green-Kubo relation, which is equivalent to the Green-Kubo relation and is the basis for calculating self-diffusion coefficients in equilibrium molecular dynamics simulations. In this paper the Einstein relation is used to calculate the self-diffusion coefficient.
In the equation ⃗ is the position vector of atom i, d represents the analogue dimension, Represents the number of particles.

The Darken equation
Darken's equation was proposed to describe the diffusion of binary alloy systems but has since been widely applied to other binary systems. The expression is.
In the equation, is the activity of the component, is the molar fraction of the components. The activity correlation coefficient can be calculated using the definition of activity in combination with the van der Waals equation of state with the following equation Error! Reference source not found. .

Simulation Setting
Calculation of the diffusion coefficient of the CO2-N2 binary system at various temperatures and at atmospheric pressure (1 atm), taking into account the molecular structure and Coulombic potential; total number of molecules is 4096, total number of atoms is 12288. The temperature and pressure were maintained constant in the isothermal-isobaric ensemble using the Nose-Hoover thermostat and the cut-off radius is 12Å. The time coupling constant is 100 fs and the time step is 2 fs. The first 200 thousand steps are used for system equilibrium and the last 500 thousand steps are used for calculation of the fit.
For comparison, we performed a control calculation without considering the molecular structure and the Coulomb potential, where the LJ parameters used are shown in Table 7. The diffusion coefficient of the CO2-N2 binary system was calculated under the same conditions and the total number of atoms was 27000.