Molecular dynamics comparative study of methane–nitrogen and methane–nitrogen–ethane systems

1 Introduction

Molecular dynamics (MD) is a talented technique which consists in modelling a simulation of macroscopic systems involving a few molecules. Nowadays, molecular dynamics tends to become an alternative to experiments in order to provide transport properties (Morriss et al., 1991; Padilla and Toxvaerd, 1992). These processes happen in many cases extremely variables in the time and the space. Because these phenomenons are irreversible, they occur in material systems at the equilibrium state. Favour to thermal agitation, they can transport a molecular size like as matter, energy and quantity of movement.

The liquefied natural gas (LNG) is a subject of numerous studies. The LNG properties have frequently been investigated because of their industrial importance. LNG consists predominantly of methane (95%). It will contain higher fraction of ethane (3%) and some propane. We quote the works of Murad et al. (1979) which concern the study of the liquid methane by simulation of molecular dynamics. We also find the works of simulation by MD of Habenschuss et al. (1981) whom studied the site–site interaction of the liquid methane.

Recent studies of thermodynamic, structural and of transport properties by MD simulation have just been finalized the modelling and statistical simulation methods of quantum chemistry. These computations involve the study of the rigid methane molecule in the spherical interaction (Tchouar et al., 1998), and in the site–site interaction (Belkacem et al., 2005). The Gibbs ensemble is a simulation technique developed specifically for the study of phase equilibrium prediction of thermodynamic and structural properties of mixtures by different numerical methods as Monte-Carlo (MC) (McDonald and Singer, 1972; Sesé, 1992; Kincaid and Scheraga, 1982) and molecular dynamics (MD) (Murad et al., 1979; Nicolas et al., 1979). Different approaches derived from Monte-Carlo simulation involve the Quadratic Feynmann–Hibbs (QFH) (Sesé, 1993), the Gaussian Feynmann–Hibbs (GFH) effective potentials, the Path-Integral Monte-Carlo (PIMC) (Berne and Thirumalai, 1986), the Path-Integral Brownian Dynamics (PIBD) (Singer and Smith, 1988), and the other Path-Integral methods (PIBC) (MS–PI) (Powles and Abascal, 1983). In these calculations, different properties are obtained for the methane interactions and supposed that the system is spherical (Tchouar et al., 2003).

In our earliest work, we have studied the different properties of CH₄–N₂ system by MD at the SP₁, SP₄, and SP₈ points (Mesli and Mahboub, 2010). Here, we apply the same method for complex systems. So, we are interested in binary (CH₄–N₂) and ternary (CH₄–N₂–C₂H₆) mixtures from the SP₁ to SP₉ points. We detail in this article the various techniques of simulation which make it possible to study the phenomenon of transport. We will present here recent results on the thermo-diffusion of binary and ternary mixtures of hydrocarbons. To perform this search and to approach the reality, we have taking into account of all the site–site interactions between molecules.

This study is realized on the mixtures of methane-nitrogen and methane-nitrogen-ethane systems. In this case, the study is very complex because it concerns to build a model of 2 sites for the nitrogen, 5 sites for the methane and 8 sites for the ethane. So, we have numerically simulated the site–site interactions of binary and ternary mixtures to visualise the evolution of our systems in the time. Moreover, our study is carried simulations for both canonical NVT and isothermal-isobaric NPT systems in order to evaluate thermodynamic, structural and transport properties of LNG constituents and to optimise the liquefaction process.

In this work, we have applied the new Buckingham potential model to binary (CH₄–N₂) and ternary (CH₄–N₂–C₂H₆) systems. We have preferred the Buckingham model because all site–site interactions of molecules (C⋯C, C⋯H, H⋯H, N⋯H, N⋯N and C⋯N) in the mixture are kept of bond environments comparatively of classical approach which considered only the spherical interactions. The new results are compared with those of literature, especially with Williams work relative to lattice energy of hydrocarbons.

In Section 2 of this paper, we give the numerical model for the Buckingham potential. We have determined Buckingham potential parameters by adjusting the two curves (LJ and Buckingham) until the total minimization of potential energies have been occurred. Radial distribution functions are described in Section 3 to determine the structural properties of LNG constituents. In the same section, the properties of transport are calculated from temporal correlation functions. In the next section, we have studied the evolution of the various thermodynamic properties of binary (CH₄–N₂CH₄–N₂) and ternary (CH₄–N₂–C₂H₆) systems, for two units NVT and NPT from the SP₁ to SP₉ points. We give our comparative study in Section 4 and we present the conclusion in Section 5.

2 Numerical model

In literature, much potential are proposed (Sesé, 1991; Catlow and Harker, 1975; Nagy et al., 1995). The first potential having a theoretical significance is the inverse potential (Maghri and Jalili, 2004). Hoover has used the nine, six, and four inverse potentials to study the transition phase fluid-solid (Evens and Morriss, 1984). After, the searchers have employed the square well potential. Consider it insufficiency; scientists have preferred Lennard-Jones potential for the spherical interactions (Lennard–Jones, 1924) and Buckingham potential for the site–site interactions (Mirsky et al., 1978; Buckingham, 1938). There are two other potentials: Exponential-6 (Evens and Hoover, 1986), and Kihara (1963) which are recommended in the spherical interaction studies (Shadman et al., 2009).

From these tools, we can simulate the macroscopic phenomenon in order to understand the physical properties. The simulation must establish a band between the macroscopic properties and intermolecular forces. So, the interaction potential model used must approach the real physical system. In the present work, we have studied methane-nitrogen and methane-nitrogen-ethane systems by MD using the new model based on the Buckingham potential. Then, we have compared the results between binary and ternary systems.

In this work, we have studied the methane-nitrogen and the methane-nitrogen-ethane mixtures which consist of 256 molecules in freely rotating state. In these calculations, all site–site interactions (C⋯C, C⋯H, H⋯H, N⋯H, N⋯N and C⋯N) between any two separated molecules are considered through the site–site model of Buckingham (Buck) model.

We have preferred the potential of Buckingham (Buck) because it: (i) takes of all bond interaction environments, (ii) has been frequently used to describe the non-bonded energy interaction between united atoms, and (iii) used for the study of the complex systems (Eq. (1)):

We have calculated the new A, B, and C parameters by our method which we describe hence. Williams has only determined the parameters A, B, and C for C⋯C, H⋯H, and C⋯H interactions (Narten and Levy, 1972) from crystal lattice energies of hydrocarbons (Williams, 1967). The distance between the C and H sites is 1.026 Å in Williams’s model. This value is less than the distance between the C and H nucleus (the usual bond length is 1.094 Å).

First, we have taken the Lennard-Jones potential parameters: σ (interaction diameter) and ɛ (depth of well) from literature (Lennard-Jones, 1924), for estimate the C⋯C, H⋯H, and C⋯H interactions (Eq. (2)):

Secondly, we have used from literature the parameters σ and ɛ, for N⋯N interaction (Allen, 1987). To determine the new σ and ɛ parameters for C⋯N and N⋯H interactions, we have used the standard Lorentz–Berthelot combining rules to calculate the interactions between unlike molecules:

Finally, to estimate the new Buckingham A, B, and C values, we have adjusted by an iterative algorithm with a mean square method, the two curves (LJ and Buckingham) until the total minimization of potential energies has been occurred. In Table 1, we give all the site–site potential parameters in real units.

The thermodynamic functions (temperature and density) are calculated as follow: T∗ = T·K_B/ɛ, ρ∗ = σ³, ρ = ∑ρ_iX_i, and ∑X_i = 1, for the nine state points of the methane-nitrogen and methane-nitrogen-ethane mixtures (Table 2).

• K_B: Boltzmann constant, ρ_i: density for constituent i and X_i: molar fraction for constituent i.

• CH₄–N₂: ${\text{X}}_{{\text{CH}}_4}$  = 0.600 and ${\text{X}}_{{\text{N}}_2}$  = 0.400.

• CH₄–N₂CH₄–N₂–C₂H₆: ${\text{X}}_{{\text{CH}}_4}$  = 0.500, ${\text{X}}_{{\text{N}}_2}$  = 0.200 and ${\text{X}}_{{\text{C}}_2{\text{H}}_6}$  = 0.300.

• MCRB: ${\text{X}}_{{\text{CH}}_4}$  = 0.42 and ${\text{X}}_{{\text{N}}_2}$  = 0.03.

• MCRT: ${\text{X}}_{{\text{CH}}_4}$  = 0.55, ${\text{X}}_{{\text{N}}_2}$  = 0.03 and ${\text{X}}_{{\text{C}}_2{\text{H}}_6}$  = 0.4

The cutoff distances for the Lennard–Jones intermolecular potentials were always 2.5σ_ij. The overall runtime is 100 ps with an equilibration period of 50 ps and data production period reaching up to 50 ps to ensure a good accuracy in calculated ensemble averages of most properties. Calculations are performed for both the microcanonical NVT ensemble and the isothermal-isobaric NPT ensemble. In the previous, the Nosé thermostat process is applied, whereas the isothermal-isobaric NPT incorporate the Parinnello–Raman integration scheme. The equations of motion (both for the translation and rotation) are solved using the fifth order Gear predictor and periodic boundary conditions around the central cubic box and the minimum image truncation were included in the calculations, long-range corrections are also applied. The integration scheme uses a constant time step algorithm Δt = 5 fs.

When the equilibrium is reached, several properties are calculated: pressure P, energy of configuration U, enthalpy H, and pair correlation function (Allen, 1987). The self-diffusion coefficient D was obtained from the slop of the mean square displacement as function of time, at long time (Hansen and McDonald, 1986). The viscosity coefficient η was evaluated from Green–Kubo integral of the computed stress correlation function (Hansen and McDonald, 1986).

3 Results and discussions

We have realised our simulation from the SP₁ to SP₉ points taken from the phase diagram (Tchouar et al., 1998) using the Buckingham potential model. The mixtures of methane–nitrogen and methane–nitrogen–ethane are taken to be rigid. All the site points are covered by molecules to interact through the Buckingham site–site potential. Then, we have compared the obtained properties for the binary and ternary systems.

Using an approach by numerical simulation of molecular dynamics, we have studied the evolution of the various structural, transport and thermodynamic properties of the CH₄–N₂ and CH₄–N₂–C₂H₆ mixtures in both NVT and NPT systems using the cited potential model.

3.1

3.1 Structural properties

The structural properties g(r) are calculated to describe the average structure of the fluid. The site–site correlation function g(r) is proportional to the probability density of finding the β site of some molecule at a distance r_αβ from the α site of some different molecule. It is given in terms of the angular pair correlation function (Allen, 1987; Streett and Gubbins, 1977; Lowden and Chasndler, 1974).

The equation of g(r) can be easily obtained from equation cited in reference (Lowden and Chasndler, 1974)

(4)

$$g(r_{\alpha \beta })=<g(r_{\alpha \beta }+r_{c\alpha 1}-r_{c\beta 2}\text{;}w1w2){>}_{w1w2}$$ where r is the vector from site α in molecule 1 to site β in molecule 2, r_αβ is the vector from the centre of molecule 1 to site α, and <⋯>_w1w2 is an outweighed average over the molecular orientations. The resulting site–site correlation functions g(r) are shown in Fig. 1.

Close

Figure 1

Radial distribution functions of site–site interactionsite–site interaction in NVT system at SP₁ point. (a) CH₄–N₂. (b) CH₄–N₂–C₂H₆.

Export to PPT

All these figures give us significant structural information for both the NPT and NVT systems. The radial distribution function forms inform us that all the interactions occur in liquid state. For the two figures (Fig. 1a and b), the NPT system characterizes much better as for the binary as the ternary system than the NVT one. So, we have obtained similar results with those of Murad et al. (1979). In addition, the curves of the ternary system present many sketches which are absent in the binary system. Ternary mixture has multi-components with numerous interactions.

In Fig, 2, we have given the most important orientations responsible for the peaks. So, we have only drawn the configurations which reproduced the picks of radial distribution functions in Fig. 1 obtained by MD. These configurations were drawn using the ACD/ChemSketch.

Close

Figure 2

Main pair interactions of two molecules: (A) CH₄–CH₄, (B) C₂H₆–CH₄, and (C) N₂–CH₄.

Export to PPT

For the H⋯H interaction, we have considered two molecules of methane or ethane in position such two atoms of hydrogen and β are in contact. The distance r₁ between the α hydrogen in molecule1 and the atom γ hydrogen in molecule 2 (with α and are in contact) reaches its maximum. This situation corresponds to the peak in g(r)_HH. The highest separation is happening when the three atoms: α,β and γ are collinear (configurations 1-A and 1-B in Fig. 2). This separation is self-sufficient from the orientation of molecule1 around its α hydrogen atom.

The configurations responsible of the peaks in Fig. 1 are shown in the representations 2-A, 2-B and 3-A, 3-B for the H⋯C and respectively C⋯C interactions which are represented by the distances r₂ and r₃ (Fig. 2).

From the configurations 1-C, 2-C and 3-C, we observe that the distances r₄ and r₅ are independents of molecule orientation around its site axis. These are relatives to the C⋯N and H⋯N interactions.

From these results, we confirm that the discontinuities and the peaks in the site–site correlation functions are the consequence of the intermolecular potential Buckingham model. So, these characteristics are obvious in the ternary system because the Buckingham potential takes of all bond interaction environments.

The average number neighbours n(R) for a molecule is obtained from the following equation:

(5)

$$n(R)=4\pi \rho {\int }_0^R{r}^{2}g(r)\mathit{dr}$$ where the integration limit R is taken as the position of successive minima in g(r).

When the value of r becomes close to the diameter of collision, g(r) increases quickly until the maximum r = r_max. When the separation r continues to increase g(r) decreases which imply that at long distance the influence of the central molecule disappeared. The advantages of these curves allow us to calculate the position of the various maxima and minima and the number of close neighbour. This last one is obtained by integration of the radial distribution function.

The obtained value of this number is compared to the experience measured by X-rayon diffraction which giving a first peak equals to 12 and the second peak equals to 55 (Habenschuss et al., 1981). Sesé (1991) has obtained him results from quantum simulations. The first maximum is R_max1= 4.05 Å for the three state points. The second maximum is R_max2= 5.75 Å for SP₁ and R_max2 = 5.85 Å for SP₄ and SP₈. The first minimum for the three state points is R_min1 = 7.75 Å. or the binary and ternary systems, the number of nearest neighbours is 26.25 molecules and in the second layer the number of neighbours is 56.55 molecules.

In the present calculations, similar radial distribution functions for the points SP₁, SP₄, SP₈ are compared with an identical location of the peak positions. The first maximum is R_max1 = 4.236 Å for the binary system, and R_max1 = 4.0.25 Å for the ternary system. The second maximum is R_max2 = 7.354 Å for the binary and R_max2 = 7.556 Å for the ternary systems. However, the first minimum is R_min1 = 5.456–5.542 Å for the three state points of binary system. The first minimum is R_min1 = 5.652–5.765 Å for the three state points of ternary system. Our simulation results of both the binary and ternary systems agree very well with experiment and theoretical works. All these results are regrouped in Table 3.

Table 3 Values of maximum and minimum for g(r) C⋯C.

State point	Rmax1 (Å)	Rmax2 (Å)	Rmin1 (Å)
SP1, SP4, SP8a	4.100	7.700	5.700−5.750
SP1, SP4, SP8b	4.050	7.750	5.750−5.850
SP1, SP4, SP8c	4.050	7.750	5.750−5.850
SP1, SP4 SP8d	1.050	1.750	–
SP1, SP4, SP8e MCRB	4.236	7.354	5.456−5.542
SP1, SP4, SP8e MCRT	4.025	7.556	5.652−5.765

MCRB: multi-component refrigerants for a binary mixture.

MCRT: multi-component refrigerants for a ternary mixture.

a Site–site interaction (Belkacem et al., 2005).

b Spherical approximation (Tchouar et al., 1998).

c Sesé work (Sesé, 1992, 1994).

d Experimental work (Sesé, 1992, 1993, 1994).

e Present work.

3.2

3.2 Transport properties

The properties of transport are studied by MD and determined by the formalism of Green–Kubo (Sesé, 1994). The last one describes the phenomenological coefficients mainly the coefficient of diffusion and the viscosity. The two parameters are calculated from temporal correlation functions.

The viscosity η is calculated from the integration of the stress autocorrelation function. For the precise values, the viscosity can be calculated by averaging these results from the stress tensor. We have calculated the viscosity as follow:

(6)

$$\eta =\frac{V}{10K_{B}T}{\int }_0^{\infty }<P(t)P(0)>\mathit{dt}$$ where V is the system volume, T the temperature, K_B the Boltzmann constant, and P the symmetric traceless pressure tensor (Dysthe et al., 1998, 1999; Cooper, 1965). The last factor allows obtaining the dynamic viscosity of environment.

The stress autocorrelation functions of binary and ternary mixtures are reported in Fig. 3. In the binary system and at the state point SP₁ (Fig. 3a), the curves provide a stable situation in the time because of the increased iteration number. The situation is completely reversed in the ternary (Fig. 3b). The stress autocorrelation function becomes negative for the molecular displacement around the YZ and ZZ stress tensor elements. We think that this situation presents an instable state at small period. We can confirm that the binary system gives good results with a high precision than the ternary. Therefore, the stress autocorrelation function sums all particles in MD box.

Close

Figure 3

Stress autocorrelation function at the state point SP₁. (a) Binary system. (b) Ternary system.

Export to PPT

In Fig. 4, we give the viscosity of CH₄–N₂ and CH₄–N₂CH₄–N₂–C₂H₆ systems at the state point SP₁. The viscosity of binary system is biggest than in the ternary system. The last one is more dense, and it contains numerous compounds ( ${\rho }_{{\text{CH}}_4-{\text{N}}_2}$  = 0.2109 g/cm³ and ${\rho }_{{\text{CH}}_4-{\text{N}}_2-{\text{C}}_2{\text{H}}_6}$ = 0.5439 g/cm³).

Close

Figure 4

Viscosity of CH₄–N₂ and CH₄–N₂–C₂H₆ in NVT system at the state point SP₁.

Export to PPT

The coefficient of distribution gives the information of diffusion of molecule for the liquid state. There values are positioned between 10⁻⁹ and 10⁻¹⁰ m²/s. At high temperatures, the molecules have an important kinetic energy which increases the distribution of the particles of methane in the mixture liquid.

The self-diffusion coefficient of liquid mixture D is obtained from the mean square displacement (MSD) at long times by the relation of Einstein (Allen, 1987) and the Green–Kubo method (Dysthe et al., 1998, 1999; Cooper, 1965)

(7)

$$D=\frac{1}{6N}\underset{t\to \infty }{\mathit{Lim}}\frac{d}{\mathit{dt}}\langle \sum _i^N[r_i(t)-r_i(0){\hspace{0.17em}]}^2\rangle$$ where the r_i(t) are the positions of particles at time t and N is the number of molecules in the system.

The mean square displacements of methane-nitrogen and methane-nitrogen-ethane mixtures at point SP₁ are reported in Fig. 5. The variation of the autocorrelation function is observed with linearly rising MSD until the equilibrium is reached. As time increases, an important gap appears between the two curves for the isothermal-canonical system (NVT).

Close

Figure 5

Mean square displacement of CH₄–N₂ and CH₄–N₂–C₂H₆ in MCR at the point SP₁.

Export to PPT

Through the binary CH₄–N₂CH₄–N₂ system, we observe a precise molecular displacement. We note also that this displacement is sluggish than one obtained with the ternary CH₄–N₂CH₄–N₂–C₂H₆ system. So, the simulation time is large and important for the calculations in CH₄N₂ system. In this case, the ternary model is very complex because it concerns of 2 sites for the nitrogen, 5 sites for the methane and 8 sites for the ethane. So, this system takes into account of all theses site interactions of molecule environments.

The formal Green–Kubo theory for D_ij in multi-component mixtures has been developed by Zhou and Miller (1996). We have used Eq. (8) in our calculations.

(8)

$$D_{\mathit{ij}}=\frac{V}{3{\mathit{NW}}_i{W}_j{m}_i{m}_j}{\int }_0^{\infty }<J_i(t)J_j(0)>\mathit{dt}$$ where V: system volume, N: total number of molecules, W_i: mass fraction, m_i: molecular mass, and J_i: mass flux of component i.

We have done the MD calculations for determining the evolution of self-diffusion coefficients of mixtures CH₄–N₂CH₄–N₂ and CH₄–N₂CH₄–N₂–C₂H₆ for the two NVT and NPT systems at the point SP₁ (sensitive state).

In the NPT ensemble, the values of diffusion coefficient are better than in NVT ensemble. In the first system, the volume fluctuates, the molecules diffuse easily with an important rate through the liquefied natural gas and the pressure is fixed along the simulation time.

Different results of nitrogen-n-pentane are established by two different models AUA (anisotropy united atom) and OPLS (optimised potentials for liquid simulation) (Rowlinson, 1969). There is a qualitative difference between the two models.

The OPLS model has a LJ spherical interaction centred on the carbon with a large diameter than the AUA model which has the site–site interaction displaced to the geometric centre of the valence electrons for the CH group (Erpenbeck, 1988).

For the mixture n-pentane–nitrogen, and under the following experimental conditions (T = 325.25 K, d = 0.711 g/cm³, P = 120.5 MPa), the values of coefficient of diffusion found by OPLS and AUA are respectively (3.44 ± 0.05)10⁻⁹m²/s and (3.87 ± 0.04)10⁻⁹m²/s, and those of viscosity are respectively (0.39 ± 0.02) MPa s and (0.31 ± 0.02) MPa s. We carried out our simulation at the same point states (temperature, pressure and density), for the two systems methane-nitrogen and methane-nitrogen-ethane in order to compare our result with the contribution quoted before.

Our calculations give us, the values of the coefficient of diffusion for the binary and ternary mixtures are respectively equal to (3.620 ± 0.030) 10⁻⁹m²/s and (3.705 ± 0.030) 10⁻⁹m²/s, and those of viscosity are correspondingly (0.257 ± 0.045) and (0.107 ± 0.033) MPa s. We note that ours results for the diffusion are intermediaries between the results of OPLS and AUA.

For the CH₄–N₂ mixture, the diffusion value decreases because it is more viscous than the CH₄–N₂–C₂H₆ mixture (Fig. 6). From these results, we confirmed that more the fluid is dense more it diffusion coefficient in the liquid will be strong. So, the diffusion of methane in the ternary mixture is much better than in the binary mixture. However, the situation for viscosity is reversed.

Close

Figure 6

Diffusion coefficient as a function of temperature for the methane–nitrogen, methane–nitrogen–ethane and n-pentane–nitrogen mixtures.

Export to PPT

4 Comparative study

In general, our calculated values are in perfect concordance with those obtained by Sesé (1992). For the points SP₄ and SP₈, our results have a big similarity with the path-integral method (PIMC) for the ternary system. Therefore in the binary system, the situation is comparable to the quantum approach of Wigner-Kirkwood (WK) (Sesé, 1991, 1992, 1993; Kincaid and Scheraga, 1982; Nicolas et al., 1979; Berne and Thirumalai, 1986; Singer and Smith, 1988; Powles and Abascal, 1983; Tchouar et al., 2003; Mesli and Mahboub, 2010; Catlow and Harker, 1975; Nagy et al., 1995; Maghri and Jalili, 2004; Evens and Morriss, 1984; Lennard-Jones, 1924; Mirsky et al., 1978; Buckingham, 1938; Evens and Hoover, 1986; Kihara, 1963; Shadman et al., 2009; Narten and Levy, 1972; Williams, 1967; Allen, 1987; Hansen and McDonald, 1986; Streett and Gubbins, 1977; Lowden and Chasndler, 1974; Sesé, 1994). However, we have examined the exception for the point SP₁ for the two systems. The point SP₁ is very sensitive state because it is near the triple point. Here, we can say that our molecule system is in a transition phase which corresponds to the passage from the gaseous to liquid state. Thus, this transition causes an unbalance into the system. Consequently, our SP₁ values differ from the others.

In addition, our thermodynamic properties especially reduced energy of configuration U∗ are in perfect agreement with the experimental data (Sesé, 1993, 1994; Singer and Smith, 1988). However, the reduced total energy E∗ agree very well with site–site interaction (Belkacem et al., 2005) and spherical approximation (Tchouar et al., 1998). The results of the dynamic properties: the coefficient of distribution and the viscosity are relative to the stream of particles. These dynamic properties are calculated from the functions of temporal correlation. The increase of the temperature induces a net increase of the total energy, the potential energy and the average enthalpy. Therefore, we notice that the values of the pressure are irregular and vary with the temperature and the density of the system. The fluctuations of these properties are weak if compared to classical model (Tchouar et al., 1998). (Table 5).

5 Conclusion

In this work, we have studied many properties for methane–nitrogen and methane–nitrogen–ethane mixtures by molecular dynamics. We have chosen the MD method to predict much better the LNG characteristics. The present calculations concern the site–site interaction study of 256 particles using the adapted Buckingham potential model. This adapted model allows determining the characteristics for each state point. We have applied this model to study the LNG properties for methane–nitrogen and methane–nitrogen–ethane mixtures by molecular dynamics.

We have calculated the thermodynamic, dynamic and structural properties for both the NVT and NPT ensembles of binary and ternary systems from the SP₁ to SP₉ points. Then, we have compared the results between binary and ternary systems. We have obtained a good prediction on transport properties for binary and ternary mixtures. The diffusion of methane in the ternary mixture is much better than in the binary mixture. However, the situation for viscosity is reversed.

From the calculated values of self-diffusion coefficient and viscosity, we confirmed that LNG mixture is in liquid state.

We conclude that our simulation model approach very well the experimental data. Our work presents the possibility to determine with a high precision the thermodynamic, dynamics, and structural properties of binary and ternary systems and to propose a good optimisation of the liquefaction process. We hope that this model could be an effectively starting material to study the properties of other complex systems in order to predict the transport phenomenon in the fluids.