Luchao Jin,Yongming He,Guobing Zhou,Qiuhao Chang,Liangliang Huang,Xingru Wu,*

1 Alchemy Sciences Inc.,6002 Rogerdale Rd,Ste 125,Houston,TX 77072,USA

2 College of Energy,Chengdu Science and Technology University,Chengdu 610059,China

3 University of Oklahoma,100 East Boyd st,Sarkeys Energy Center 1210,Norman,OK 73019,USA

Keywords:High-pressure high temperature Z-factor Molecular dynamics simulation Natural gas density Correlations

ABSTRACT This work applied molecular dynamics (MD) simulation to calculate densities of natural gas mixtures at extremely high pressure (＞138 MPa) and high temperature (＞200°C) conditions (xHPHT) to bridge the knowledge and technical gaps between experiments and classical theories.The experimental data are scarce at these conditions which are also out of assumptions for classical predictive correlations,such as the Dranchuk &Abou-Kassem (DAK) equation of state (EOS).Force fields of natural gas components were carefully chosen from literatures and the simulation results are validated with experimental data.The largest relative error is 2.67% for pure hydrocarbons,2.99% for C1/C3 mixture,7.85% for C1/C4 mixture,and 8.47% for pure H2 S.These satisfactory predictions demonstrate that the MD simulation approach is reliable to predict natural-and acid-gases thermodynamic properties.The validated model is further used to generate data for the study of the EOS with pressure up to 276 MPa and temperature up to 573 K.Our results also reveal that the Dranchuk&Abou-Kassem(DAK)EOS is capable of predicting natural gas compressibility to a satisfactory accuracy at xHPHT conditions,which extends the confidence range of the DAK EOS.

1.Introduction

The sharp increase in energy demand has pushed oil and gas acquisition to risky and extreme conditions.High pressure and high temperature (HPHT) environments are classified into a tier system with three main tiersas shown in Table 1.In this work,we honor this tier system and define extremely high pressure,high temperature (xHPHT) conditions as the temperature is more than 200°C and pressure is above 138 MPa.In some gas reservoirs in the United States,conditions of temperatures up to 260°C and pressures up to 241 MPa have been observed[1–3].Consequently,the knowledge of natural gas density under the xHPHT condition becomes critical in evaluating the risks of well drilling,the completion,and quantifying resources in place.

Measuring the gas densities has been an endeavor of many experimentalists.Over the last decades,thousands of density data of natural gases below 150°C and 70 MPa have been obtained [4–7].However,xHPHT data are scarce due to the great challenges of conducting laboratory experiments.Very often,the reported xHPHT measurement has poor accuracy in the low-density region[8].Recently,Liu,Wu[9]developed laboratory equipment with the capability of operating at around 280 MPa and 250°C(xHPHT)and measured large hydrocarbon molecules that do not commonly exist in natural gas reservoirs.In the literature,only a handful of pure components such as pentane and propane have been studied under xHPHT conditions [2,9,10].Other than hydrocarbons,density behaviors of impurities in natural gas including carbon dioxide and hydrogen sulfide under xHPHT conditions also should be studied as they are different from hydrocarbon in molecular polarity,weight,and association.There is a knowledge gap on the density of natural gas,due to limited understanding from those xHPHT conditions.

On the other hand,molecular dynamics (MD) simulation has been widely used to predict the physical properties of hydrocarbons under various conditions [15].If the force field is capable of describing the species and is validated by data in a range of temperatures and pressures,it can be extended to predict the gas physical properties in a higher temperature and pressure with high confidence,as long as there is no phase change and chemicalreactions caused by strong molecular interactions.Secondly,in the petroleum industry,corresponding states theory such as the Hall and Yarborough (HY) method [11]and the Dranchuk &Abou-Kassem (DAK) method [12]are easier to use than the MD simulation method.However,the DAK approach was developed by fitting the Standing and Katz [13]gas Z-factor chart to high accuracy in the regions 0.2 ≤Pr≤30;1.0 ≤Tr≤3.0 and Pr＜1.0;0.7 ＜Tr≤1.0[12].There is no study to date to demonstrate those EOS methods are reliable in xHPHT conditions.In this work,we carry out MD simulations to generate density points and supplement the experimental understandings from the xHPHT region.

Table 1 High pressure and high temperature tiers [3]

2.Simulation Methods and Details

2.1.Potential model

In this work,a series of MD simulations have been carried out to investigate the density of pure natural gases of n-alkanes as well as its mixtures with various compositions under xHPHT.Also,we have explored the effects of gas impurities,like dipolar hydrogen sulfide(H2S) and quadripolar carbon dioxide(CO2),on the density of n-alkanes studied here.For these compounds,a variety of force field parameters have been proposed in previous literature.Among them,the united-atom(UA)TraPPE model[14]was used for the nalkanes and the NERD[15]and EPM2[16]model was used for the impurities H2S and CO2,respectively.In the UA TraPPE model,each CH2or CH3group of n-alkanes was represented as an uncharged pseudo atom denoted as C.In the n-alkanes,the equilibrium C-C bond length and C-C-C angle were 1.54×10-10m and 114°,respectively.Concerning the impurity H2S,the NERD model developed by Nath was used and the equilibrium H-S bond length and H-S-H angle were 1.365 ? and 91.5°,respectively.For the n-alkanes and H2S mentioned above,all the bonds and angles were constraint to their equilibrium values by using the SHAKE algorithm [17].Besides,the EPM2 model was used for CO2and the equilibrium C-O length and the O-C-O angle were 1.15 ? and 180°,respectively.Here it is worth noticing that the SHAKE algorithm couldn’t be able to constrain an angle at 180°.As a result of this,a flexible EPM2 model with an additional bond stretching constant (kb=10,739 kJ·mol-1·?-2) and angle bending constant (kθ=1236 kJ·mol-1·rad-2) was applied from Nieto-Draghi,de Bruin [18].In all cases investigated here,the nonbonded interactions between species were described by the combination of electrostatic,and Lennard-Jones (L-J) potentials and the mixed L-J parameters were obtained from the Lorenz-Berthelot mixing rules,as shown in Eqs.(1)–(2).All L-J parameters and partial atomic charges used in this work were summarized and listed in Table 2.

The reason we chose this particular mixing ratio is based on the fact the most commercially viable natural gas reservoirs have relative low concentration of impurities including CO2,H2S,or N2.This mixing rule is simple in computation and has enough accuracy in characterizing the system of natural gas mixture.

2.2.MD simulation details

All simulations were performed using the Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS) software package[19].The isothermal-isobaric(NPT)ensemble is applied where the number of molecules(N),the pressure(P),and the temperature(T)was fixed.Both the temperature and pressure were maintained constant by using the Nose-Hoover method with coupling times of 0.1 and 1.0 ps,respectively.The periodic boundary conditions were applied in all three dimensions.Newton’s equation of motion was integrated by using the velocity-Verlet algorithm with a time step of 1.0 fs.A cutoff of 1.0 nm was applied for the nonbonded interactions and the electrostatic interactions were calculated by the particle-particle particle-mesh (PPPM) method (Hockney and Eastwood).For each MD simulation,a total of 10.0 ns was run,where the first 5.0 ns was for equilibrium,and the latter 5.0 ns was used for data analysis,in which the trajectory was updated every 100.0 fs.

3.Model and Force Field Verifications

The chosen force field in MD simulations is firstly validated using laboratory data from pure components and binary mixtures of hydrocarbon.The tuned MD models are then run to predict the properties under XHPHT conditions.

3.1.Pure components

During the past years,previous studies have used the TraPPE model to calculate the PVT properties of alkanes under low pressure and temperature,but only a few studies have paid attention to this model under the extreme conditions,which promotes us to do the verifications first.On the other hand,Setzmann and Wagner [8]proposed an EOS for methane properties covering temperatures up to 625 K and pressures up to 1000 MPa,which has been used as a reference EOS for many other studies [31–33].Consequently,in the HPHT region,we directly treat the density points from Setzmann and Wagner[8]as the testing data in the following description for comparison.For propane and n-pentane,we also selected density data experimentally measured by Miyamoto and Uematsu[10]and Gamwo,Burgess[2],respectively.The temperature and pressure conditions of these data points are from the HPHT region.To compare the simulated densities with the experimental ones,we define the relative error as follows:

where ρsimand ρexpare the simulated and the experimental densities,respectively.Here for propane and pentane,the simulation is carried out in the low-pressure region because of the limit of available experimental data.While for methane,the simulated pressure reaches up to 1000 MPa,as shown in Fig.1.We can observe from this figure that the simulated densities for three pure components with UA TraPPE model increase for the pressure and decrease with the increase of temperature,which is following the corresponding results from experiments or EOS.Meanwhile,it is clear from Fig.2 that all three pure components have a relative error of less than 3.0%,which is satisfactorily acceptable when considering the systematic deviation of the critical density.Especially in the case of CH4,the relative error is only about 1.0%when the pressure even reaches 1000 MPa.As a result of this,we believe that the benchmark of this part supports the choice of both models and force field parameters.

For CO2,Zhang and Duan (2005a) have compared their calculated results with the experimentally measured specific volumeof CO2for the pressure from 70 to 400 MPa and the temperature from 323.15 to 748.15 K.In their studies,they have compared several models,including MSM [20],EPM2 [16],TraPPE [21]and Errington [22],for CO2and finally they found that the Errington model shows the best accuracy.However,one of the drawbacks for Errington is that the exponent-6 potential framework hinders the extension to mixtures.Instead,the EPM2 model also could produce very satisfactory accuracy with a relative error of about 2.0%.Besides,Trinh,Vlugt [23]have verified the accuracy of flexible EPM2 model for CO2with the temperature ranging from 300 to 1000 K and pressure up to 200 MPa.Also the results show that the flexible EPM2 model gives similar accuracy at low pressures but significantly better accuracy at high pressure.In this case,we directly adopt the flexible EPM2 model for CO2to simulate the corresponding properties at xHPHT conditions.

Table 2 A list of force field parameters used in this work

Fig.1.Comparisons between the simulated and the experimental data of three pure hydrocarbons:CH4 ,CO2 and C5 H12 .

Fig.2.The relative errors of the simulated densities from the experimental data of the three studied pure hydrocarbons:CH4 ,C3 H8 ,and C5 H10 .

In the case of H2S,only a few data have been reported in previous studies.Sakoda and Uematsu[24]have developed an equation of state based on the Helmholtz free energy of the H2S fluid phase and the model has been used to predict the thermodynamic properties of H2S for temperature up to 760 K and pressure up to 170 MPa.Besides,this model has also been applied to evaluate the H2S properties in the National Institute of Standards and Technology(NIST)database.Based on the EOS results,we have studied two temperatures under four different pressures and the results are listed in Fig.3.It can be seen from this figure that the simulated results have the same variation trend compared with those of EOS,but the density values are slightly underestimated.The calculated relative error in Fig.4 shows the largest value of 8.0% under the condition of low temperature,but the value gradually becomes smaller when the pressures increase.This assures us that the model used for H2S in this study is reasonable.

3.2.Binary mixtures

In addition to the aforementioned benchmark calculations for the pure components,the MD simulations have also been employed to evaluate the accuracy of model and force field parameters describing the gas mixtures.Experimental densities of methane/propane (C1/C3) and methane/n-butane (C1/n-C4) mixtures have been measured by Lee [4]at temperatures from 310 to 510 K and pressures ranging from 0.1 to 69 MPa.Here it should be noted that the available experimental data for gas mixtures are all from low pressures.Thus,for the model validation,the simulations for gas mixtures of C1/C3 and C1/n-C4 with various molar fractions were carried out at T=510 K,which is the highest temperature of available data,and pressures only ranging from 31 to 69 MPa.

Fig.3.Comparisons between MD results and experimental data for pure H2 S densities.

Fig.4.Relative errors of the MD results and the experimental density of pure H2 S.

The relative error as a function of pressure for the simulated densities and the experimental data from Lee [4]is listed in Fig.5 which shows that all the systems have a higher relative error at low pressure(P ＜50 MPa)than that at high pressure.For example,the 0.6C1–0.4C3 mixtures have a relative error of about 3.0%at 32 MPa,but it is less than 1.0% at the pressure of 68 MPa.Overall the TraPPE model gives satisfactory predictions for C1/C3 mixtures with relative errors less than 3.0%,and the corresponding values for C1/n-C4 mixtures range from 2.0% to 8.0%.Although this deviation does not affect the choice of models and force fields in this study for xHPHT conditions,one has to pay closer attention if the low-pressure region is of interest.With the aforementioned benchmarks from Figs.2-5,we conclude that the choice of models and force fields,as shown in Table 2,provide satisfactory density predictions for both pure components and their mixtures.

4.Gas Density From MD Simulations

Based on the above validations and promising results,in this section,we have further extended the T–P conditions to the xHPHT region:up to 573.15 K(572°F)and 275.79 MPa(40,000 psi),which follows the classification of xHPHT condition in literature [3].The underlying philosophy of the prediction method is as long as there are no strong physical and chemical forces lead to phase change and chemical reactions,the density change is only a matter of intermolecular distances.In this study,we assume that there is no phase change or chemical reaction of natural gas components under the investigated temperatures and pressures.Therefore,the MD simulation method is reliable in this study to generate density values of natural gas in the xHPHT conditions.

Fig.5.Relative errors of the MD simulated densities against experimental data of the binary mixtures at 510 K [4].

Here the densities for pure components of C3 and C5,as well as the binary mixtures of C1/C3 and C1/n-C4,have been investigated and the results are shown in Figs.6-9,in which the verified results are also plotted in these figures for comparison.For the pure components of both C3 and n-C5,it can be seen from Figs.6 and 7 that the densities decrease with the increase of temperature under the same pressure,which should result from the increased intermolecular distance induced by the faster motion of the molecules under high temperature.On the other hand,the density values for each system with the same temperature increase concerning the enhancement of pressure due to the decrease of intermolecular distance under higher temperatures.Similarly,for C1/C3 and C1/n-C4 mixtures in Fig.8,the density curves display similar variation trends as we discussed for the pure components.Besides,we could also observe that the mixture density becomes larger when the molar fraction of C3 or C4 increases.This is because the C3 and C4 have larger molecular size and mass than that of C1,which leads to the mixtures with a higher fraction of C3 or C4 would have higher densities under the same volume condition.

Fig.6.Calculated densities of pure C3 from HPLT to xHPHT region(1 psi=6894.76 Pa).

Fig.7.Calculated densities of pure nC5 from HPLT to xHPHT region(1 psi=6894.76 Pa).

Fig.8.Calculated densities of hydrocarbon mixtures from LPLT to the xHPHT region (1 psi=6894.76 Pa).

Fig.9.Predicted densities of mixtures methane with impurities to xHPHT(1 psi=6894.76 Pa).

Apart from the alkane mixtures,we also have explored the effects of impurity on the density of pure alkanes.Generally,CO2and H2S are often treated as impurities of natural gases due to their relatively small quantities.The aim of simulating the nature gas mixtures containing these impurities is to evaluate their impact on EOS prediction.For each mixture,the mole fraction of C1 and impurity are fixed at 0.8 and 0.2,respectively,and the results are illustrated in Figs.9.It is obvious from this figure that the density of mixtures also increases with the pressure but decreases with the temperature as we discussed before.Besides,it could be noted that the density of 0.8 C1–0.2 CO2is higher than the corresponding of 0.8 C1–0.2 H2S.This is reasonable because the CO2has a larger molecular mass than that of H2S,which means the mixture with CO2would have a larger density than that with H2S when both these two impurities have the same mole fraction.

5.Gas Density From Equation of State Evaluation

In the petroleum industry,corresponding states theory has been widely used,for example,the Hall and Yarborough (HY) method[11]and the Dranchuk &Abou-Kassem (DAK) method [12].The DAK approach was developed by fitting the Standing and Katz[13]gas Z-factor chart in the regions 0.2 ≤Pr≤30;1.0 ≤Tr≤3.0 and Pr＜1.0;0.7 ＜Tr≤1.0 [12].However,the xHPHT conditions encountered in deep natural gas reservoirs are beyond the original training range.It is thus necessary to examine the reliability of those EOS methods before adopting them for xHPHT conditions.In this session,we compared our MD simulation results with the ones obtained from the EOS method.

A typical procedure of estimating Z-factor by the DAK EOS is(1)estimate pseudo critical pressure and temperature for given specific gas gravity using correlations.We have applied five methods as shown in Table 3;(2) remove the effects of impurities.The Wichert-Aziz correction method [25]was utilized in this work to adjust the pseudo critical pressure and temperature when the gas mixture contains H2S or CO2;(3)calculate the reduced properties by Eqs.(4) and (5) shown below;(4) apply the DAK equation(Eq.(7)) to iteratively solve Z-factor.

where Ppcand Tpcare pseudo critical pressure and pseudo critical temperature,respectively.

For comparison,the Z factor of 119 data points generated by our MD simulation was calculated as follows:

where Mgis the molecule weight of the gas mixture in g·mol-1;P is the pressure in MPa;ρ denotes the gas density in g·cm-3;R is the ideal gas constant,8.314 J·mol-1·K-1;T is the temperature in K.

The DAK equation is expressed as

Where the constants A1through A11are as follows:

To compare the MD simulation results with the estimations from the DAK method with impurity corrections,the absolute average error (AAE) is used:

Table 3 Evaluation results of DAK EOS at xHPHT conditions

in which the relative erroriis:

As shown in Table 3,the DAK EOS can provide a very good density estimation of natural gases studied.The AAE is less than 3%for various pseudo critical pressure and pseudo critical temperature methods.The largest deviation is produced from the Standing method [26],while the approach from Thomas,Hankinson [27]produces the best match with MD simulations.

The data points used in this comparison study are from Lee[4].Fig.10 shows the scatter plot off the 119 Z-factors determined by MD simulations and the DAK EOS.The data points fall on the diagonal line indicating the validity of the DAK EOS in predicting natural gas densities at xHPHT conditions up to 276 MPa and 573 K.Fig.11 plots the relative error calculated via Eq.(3),using MD results as the reference data to be compared with the DAK results.As shown in Fig.11,the relative error is larger in the lower pressure region but decreases with the increase of pressure.This is probably attributed to the large relative errors from MD simulations,as discussed in Section 3.2.Nevertheless,the DAK method can predict gas compressibility and density at xHPHT conditions to a satisfactory accuracy.With the support of MD simulations,in this work,we extend the DAK EOS to conditions where pressure is up to 276 MPa and temperature up to 573 K.

Fig.10.Comparisons of Z-factors from DAK and MD calculations.

Fig.11.Deviations provided by the DAK EOS from MD simulations.

6.Conclusions

In this paper,molecular dynamics simulations were conducted to predict natural-and acid-gases densities at pressure and temperature conditions,up to 276 MPa and 573 K,respectively.The TraPPE,NERD,and EPM2_Flex potential parameters were adopted to model the intermolecular interactions of n-alkanes,H2S,and CO2.The validity of DAK equation of state at the xHPHT conditions was also evaluated with the help of the MD simulation results.Conclusions that were drawn from this study also include:

1.The MD simulation model and force fields accurately reproduced available experimental data at low and intermediate HPHT conditions.The largest relative error is 2.67% for pure hydrocarbons,2.99% for C1/C3 mixture,7.85% for C1/C4 mixture,and 8.47% for pure H2S.These satisfactory predictions demonstrate that the MD simulation approach is reliable to predict natural-and acid-gases thermodynamic properties.The applicability of this approach is worth to be validated and extended to more complex gas mixtures when experimental data is available.

2.Overall 119 density points were generated by using MD simulation methods covering the xHPHT conditions with pressure ranging from 69 to 276 MPa and temperature ranging from 453 to 573 K.The studied systems include pure n-alkanes,C1/C3 and C1/C4 mixtures,and acid gases containing 80% C1 and 20% H2S (or 20% CO2).The developed database can be a basis for future studies of those gas mixture properties under xHPHT conditions.

3.The DAK EOS demonstrated a good match with the results from MD simulation even to the xHPHT region.The average absolute error (AAE) is 2.14% when using the Thomas,Hankinson and Phillips method to estimate pseudo critical temperature and pseudo critical pressure.With the results from this work,we extended the confidence range of the DAK EOS to pressure up to 276 MPa and temperature up to 573 K.

Acknowledgements

The authors are grateful to the partial financial support from Ballard Petroleum Holdings and Yangtze University,and the Schooner Supercomputing from the University of Oklahoma.L.Huang also wants to acknowledge the startup support from the University of Oklahoma.

Nomenclature

AAE average absolute error

Mgmolecule weight of natural gas,g·mol-1

kBBoltzmann constant,J·mol-1·K-1

P pressure,MPa

Ppcpseudo critical pressure,MPa

Pprreduced pseudo temperature,dimensionless

q point charge,e

R ideal gas constant,8.314 J·mol-1·K-1

T temperature,K

Tpcpseudo critical temperature,K

Tprreduced pseudo pressure,dimensionless

Z gas compressibility factor,dimensionless

∈ Lennard-Jones well depth,J·mol-1

ρ density,g·cm-3

σ Lennard-Jones molecular size,?