- Email: [email protected]
C10 in bulk phase and nanoslits
International Journal of Heat and Mass Transfer 141 (2019) 643–650
Contents lists available at ScienceDirect
International Journal of Heat and Mass Transfer journal homepage: www.elsevier.com/locate/ijhmt
Molecular insight into the miscible mechanism of CO2/C10 in bulk phase and nanoslits Timing Fang a,b, Yingnan Zhang a, Jie Liu b, Bin Ding c, Youguo Yan a,⇑, Jun Zhang a,⇑ a
School of Materials Science and Engineering, China University of Petroleum, Qingdao, Shandong 266580, China Department of Chemical and Biomolecular Engineering, National University of Singapore, 117576 Singapore, Singapore c Research Institute of Petroleum Exploration & Development of PetroChina, Beijing 100083, China b
a r t i c l e
i n f o
Article history: Received 16 April 2019 Received in revised form 9 June 2019 Accepted 24 June 2019 Available online 9 July 2019 Keywords: Molecular dynamics simulation Carbon dioxide Miscibility Nanoslit
a b s t r a c t The determination of minimum miscibility pressure (MMP) of CO2/oil in tight oil reservoirs has aroused huge attention from researchers. In this study, we characterized the mass transfer process and determined the MMP of CO2/oil systems in nanoslits through molecular simulation. A bulk phase and five confined nanoslits with the pore size ranging from 2 to 15 nm were examined. It was found that the MMPs in nanoslits were significantly lower than those in bulk phase, thus revealing that rising the probability of molecular interaction is a vital factor affecting the miscible phase. Furthermore, due to the consideration of adsorption layer, the width of nanoslit with the lowest MMP was larger than the reported value, which will bring more uncertainty to the molecular migration. This simulation study highlighted the crucial role of the confinement effect in the miscible process of CO2 flooding, which might facilitate the development of low permeability reservoirs. Ó 2019 Elsevier Ltd. All rights reserved.
1. Introduction Both experimental and theoretical studies proved CO2 flooding an effective method for enhanced oil recovery (EOR), which is likely to enhance oil recovery by nearly 8–14% of the original oil in place (OOIP) through an efficient mass transfer [1–4]. For the tight and shale oil reservoirs with massive micropores, CO2 flooding plays a significant role in EOR because CO2 can easily transfer into the oil phase [5,6]. The minimum miscibility pressure (MMP) is defined as the lowest operating pressure, at which oil and CO2 are mutually miscible in the reservoir. Under a higher pressure over the MMP, oil recovery will be enhanced by CO2 injection. If the reservoir pressure is below the MMP, oil and CO2 will be incompletely miscible, thereby leading to a low oil recovery [7,8]. Accordingly, the accurate determination of the MMP in a CO2/oil system is critical for EOR. In the literature, several experimental techniques and theoretical models have been reported to determine the MMPs of crude oil/solvent systems. The slim-tube tests [9], the rising bubble apparatus [10] and the vanishing interfacial tension (VIT) technique [11] are considered the commonest experimental technique. Represented by the VIT, the method incorporates various research
⇑ Corresponding authors. E-mail address: [email protected] (J. Zhang). https://doi.org/10.1016/j.ijheatmasstransfer.2019.06.083 0017-9310/Ó 2019 Elsevier Ltd. All rights reserved.
methods (e.g. theoretical analysis [12] and molecular simulation [13]) to make the measurement object wider in range and more precise. In the VIT technique, the MMP is evaluated by linearly extrapolating the interfacial tension (IFT) at CO2/oil interface to zero, i.e., the miscibility between CO2 and oil is achieved at the zero-IFT condition [14]. Besides, by combining molecular dynamics (MD) simulations and the equation of state, Andrés et al. [13] and Wang et al. [15] estimated the MMPs and then explored the IFTs between injected gas and crude oil in bulk phase. The MMPs evaluated using the VIT technique were found close to those using slim-tube test. However, it has been estimated that nearly 40% of OOIP exists in unconventional oil reservoirs, where oil is adsorbed in nanosized pores [16]. Accordingly, it is necessary to investigate the MMPs of CO2/oil systems in nanopores. However, the exact MMPs should be experimentally measure in nanopores with sizes ranging over 2–100 nm [17]. Due to the confinement and surface interaction, the physical properties of fluids in nanopores significantly differ from those of conventional macroscopic systems [18,19]. These methods are primarily used to determine MMP in bulk rather than MMP in nanopores [20,21]. In recent years, the behavior of miscible in nanopores has gradually aroused huge attention. Considering the particularity of the phase change of the oil phase in confined spaces, instead of direct measurement, the researchers used theoretic modelling and deduction. The methods developed based on the physical and structural properties of oil and gas have been employed. For instance, the
644
T. Fang et al. / International Journal of Heat and Mass Transfer 141 (2019) 643–650
diminishing interface method (DIM) proposed by Zhang et al. [22] addressed the problem that the MMP of heterogeneous interface is difficult to measure accurately. The DIM determines the relationship between the variation of thickness and the pressure, based on the stable interface thickness of miscible phase, to measure the specific MMP value. Similar studies have been published on the proof and extension of the method [8,23]. Besides, from the millimeter to the nanometer scale, the researchers adopted the Peng-Robinson equation of state (PR-EOS) to measure MMPs in nanopores through the calculation of solubility parameters and capillary pressure. The critical properties dominating the phase behavior and miscibility of confined fluids were further revealed [24,25]. The above studies laid the solid foundation of theoretical models for further studies. To explore the microscopic mechanisms not readily available in experimental and theoretical methods, MD simulation was performed to study microstructural variations and dynamics features of oil molecules in nanopores [26,27]. To explore the effects of temperature and oil density on the diffusion coefficients of alkane in nanopores, Cao et al. utilized MD simulations [26]. Likewise, Wang et al. characterized the fluid flow of CO2/octane in nanopores affected by driving force and pore size [19]. Thu et al. examined the relation between interaction of CO2/oil mixtures and dynamics of oil molecules in silica pores. Their results demonstrated the effect of molecular interactions on the adsorption and migration behavior of alkanes, thus indicating the implication of microscopic analysis [28,29]. Note that molecular simulation is a more intuitive method to examine sophisticated interactions and unravel microscopic properties [30–32]. However, on molecular and atomic level, fewer studies have focused on the microscopic mechanism of miscible behavior in nanopores, which helps enhance utilization efficiency of CO2 in tight oil exploitation. Based on MD simulations, interfacial thickness method was employed in this study to characterize the miscibility of CO2/decane (C10) systems in nanoslits and then to evaluate the MMPs. To verify the methodology, CO2/C10 systems in bulk phase were also examined. In Section 2, the molecular models of CO2/C10 systems in bulk phase and nanoslit were described, followed by simulation methods. In Section 3, first, the conventional VIT method and the interfacial thickness method were used to estimate the MMPs of CO2/C10 systems in bulk phase. Next, the interfacial thickness method was extended to evaluate the MMPs in the nanoslits. Subsequently, the miscible dynamics in different nanoslits with widths of 2, 3, 5, 9 and 15 nm was examined, respectively. Lastly, the confinement effects on miscible behavior is discussed. In Section 4, the concluding remarks are summarized. 2. Methods 2.1. Models and simulations 2.1.1. In bulk phase All the systems were constructed Materials Studio (MS) [33]. Fig. 1a shows a typical CO2/C10 system in bulk phase. Initially,
Fig. 1. A CO2/C10 system in (a) bulk phase & (b) a nanoslit. CO2: red, C10: blue. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
the central region of the system contained solely C10 molecules, which were added into a simulation cell with random orientations; then, CO2 molecules were added on each side of C10. For each CO2/C10 system, the number of C10 was set to 300 and the number of CO2 molecules was calculated by nCO2 = xCO2nt, where xCO2 is the solubility of CO2 in C10 (in terms of mole fraction) and nt is the total number of molecules. To ensure sufficient CO2 molecules interacting with C10, the experimental xCO2 at 12.62 MPa was selected [34] and then 300 more CO2 molecules were added into the system. The thickness of oil film was set as 5 nm. 2.1.2. In nanoslits To examine the confinement effects, the CO2/C10 systems were placed into a nanoslit as shown in Fig. 1b. The nanoslit was represented by two hydroxylated silica surfaces, each was cleaved from a-quartz along the (1 0 0) plane with dimensions (xyz) of 39.98 nm 2.70 nm 1.507 nm. As the nanopore size in oil reservoirs normally ranges from 3 to 8 nm [35,36], most of the confined CO2/C10 systems were simulated in a nanoslit with width of 5-nm along the z-axis in this study. Nevertheless, other nanoslits (with width of 2, 3, 9 and 15 nm) were also considered to examine the effect of pore size. The initial C10 phase with dimensions of 8.50 nm 2.70 nm 5.00 nm consisted of 274 C10 molecules. 2.1.3. Determination of pressure in nanoslits To estimate the pressure inside the nanoslit, the system was divided into small bins each of 0.5 Å along the z-axis. The pressure within a bin was calculated by
Pbin ¼
Sxx þ Syy þ Szz 3V
ð1Þ
The symmetric stress tensor of atom i, Sab, was given by
2
Sab
3
Np P
ðr 1a F 1b þ r 2a F 2b Þ 7 6 mv a v b þ 7 6 n¼1 7 6 7 6 N b 7 6 1 P ðr1a F 1b þ r2a F 2b Þ 7 6þ2 7 6 n¼1 7 6 7 6 Na 7 6 1 P ðr1a F 1b þ r2a F 2b þ r3a F 3b Þ 7 6þ3 7 6 n¼1 7 ¼ 6 7 6 Nd 7 6 1 P ðr1a F 1b þ r2a F 2b þ r3a F 3b þ r4a F 4b Þ 7 6þ4 7 6 n¼1 7 6 7 6 Ni P 7 6 1 6þ4 ðr1a F 1b þ r2a F 2b þ r3a F 3b þ r4a F 4b Þ 7 7 6 n¼1 7 6 7 6 N p P 5 4 r ia F ib þKspaceðr ia ; F ib Þ þ 1 2
ð2Þ
n¼1
where a and b take values of x, y, z to produce six components of a symmetric tensor. In this equation, the 1st term on the right-hand side is the kinetic energy contribution from atom i; the 2nd term is a pairwise energy contribution, where n iterates on the Np neighbors of atom i, r1 and r2 are the positions of the two atoms, and F1 and F2 are the forces on the two atoms resulting from the pairwise interaction; the 3rd to 8th terms are contributions of the N b bonds, Na angle, Nd dihedral, Ni improper, long-range Coulombic interaction (Kspace term) and internal constraint force, respectively. This method has been widely employed to estimate the pressure in a confined system [37,38]. To demonstrate and compare with available simulation results in the literature, a system consisting of CH4 molecules in a graphitic nanoslit with width of 3.6 nm was constructed and simulated. Fig. S1 shows the density and pressure profiles of CH4 in the nanoslit along the z-axis (perpendicular to the pore surface). In the middle of the nanoslit, both density and pressure remain as a constant, representing ‘‘free gas” with negligible interaction with the pore
645
T. Fang et al. / International Journal of Heat and Mass Transfer 141 (2019) 643–650
surface. As listed in Table S1, the estimated density and pressure of free CH4 match very well with a reported simulation study [37]. For the current study, as shown in Fig. S2, a region of constant pressure is also observed for CO2 in a nanoslit. 2.1.4. Simulation details MD simulations were performed by large-scale atomic/molecular massively parallel simulator (LAMMPS) software ((Version 16 Mar 2018)) [39]. Based on practical oil reservoir conditions [40], the temperature was kept at 344.15 K and controlled by a NoséHoover thermostat [41]. The cutoff distance to evaluate the LJ interactions was 1.0 nm and the Ewald summation method was used for the electrostatic interactions. The periodic boundary conditions were applied in all the directions. In this work, we carried out all-atom MD simulation of each system, the simulation duration was 10 ns. The time step was 1 fs and every 1000 steps were saved for data analysis. After the simulation, all the snapshots were displayed by VMD software [42]. Bulk phase: the simulations of CO2/C10 systems in bulk phase were performed in a NPT ensemble (with constant number of molecules, pressure and temperature). Initially, the dimensions of each system were set to be 50 Å in the x and y directions. Nanoslits: the hydroxylated silica surfaces were described by the CLAYFF force field [43]. The simulations were conducted in a NVT ensemble (with constant number of molecules, volume and temperature). Identical to those in bulk phase. To obtain an initial equilibrated C10 phase, 5 ns simulation was first performed. Then, CO2 phase with different densities, were placed on each side of C10 phase for 10 ns simulation. The non-bonding interaction consists of the short-range van der Waals (vdW) and long-range electrostatic interaction. The vdW is represented by Lennard-Jones (LJ) potential and the electrostatic interaction is represented by Coulombic potential,
"
V non-bonded ¼ 4eij
rij
12
r ij
rij
rij
6 #
þ
qi qj 4pe0 r ij
ð3Þ
where eij and rij are the well depth and collision diameter, rij is the distance between atoms i and j, qi is the atomic charge of atom i, and e0 = 8.8542 1012 C2 N1 m2 is the permittivity of vacuum. For the bonding interaction, there are three components including bond stretching, angle bending and dihedral torsion. The specific parameters of force field are listed in Section S2 in the Supporting Information. 2.2. Interfacial tensions (IFTs) The interfacial tension was calculated by the formulation of the Gibbs interfacial tension in LAMMPS. For CO2/C10 systems in bulk phase, two interfaces are both perpendicular to the x-axis and parallel to the zy-plane, and the IFTs between CO2 and C10 phase were calculated by the expression of pressure tensor [44,45]:
1 c¼ 2
ZLx
Pyy þ Pzz 1 Lx Pxx ðPN ðxÞ PT ðxÞÞdx ¼ 2 2
ð4Þ
where PN and PT are the normal and tangential pressures, respectively; Paa (a = x, y, z) is the diagonal element of the pressure tensor; and Lx is the box length in x-axis. The pressure components were evaluated by the virial equation: N X i¼1
mi v i;a v i;b þ
N1 X N X i¼1
rij;a f ij;b
box. Finally, the average of IFTs during the last 2 ns was taken as the final value. 2.3. Interfacial thicknesses In both bulk phase and nanoslits, there is no distinct interface between CO2 and oil phase; moreover, the density of oil is not uniform along the x-axis. To estimate the CO2/oil interfacial thickness, a recently proposed characterization method for polymer/solvent systems [46] is adopted here. As illustrated in Fig. 2, based on the concept of Gibbs dividing surface (GDS), the density profile of C10 along the x-axis can be fitted by [46]:
1 2
qðxÞ ¼ q0 erf
x x1 pffiffiffi r 2
erf
x x2 pffiffiffi r 2
ð6Þ
where q0 is the bulk density away from the CO2/C10 interface, x1 and x2 represent the positions of Gibbs surfaces; r represents the standard deviation in length. d is the length (x2-x1) between two Gibbs surfaces. As the Gibbs surface is an idealized boundary separating two phases, d can be considered as the effective thickness of oil phase in the system. Similar methods were applied to determine the locations of fluid/solid interfaces [47–50]
3. Results and discussion 3.1. In bulk phase
0
Pab V ¼
Fig. 2. Density profile of C10 along the x-axis. The equilibrium snapshot for CO2/C10 in bulk phase is shown on the top.
ð5Þ
j>i
where vi, a the velocity of atom i in direction a; rij = ri rj; fij the force exerted on atom i by atom j; and V is the volume of simulation
3.1.1. MMP determination by VIT Fig. 3 shows the IFT between CO2 and oil versus pressure in bulk phase at 344.15 K. The predicted IFT in this study matches well with the experimental study by Shaver et al. [40]. With rising pressure, the IFT monotonically decreases and has a linear relationship with pressure. Through linear regression, we find IFT = 19.61– 1.62 P from our simulation and IFT = 19.14–1.63 P from the experiment [40]. The MMP can be evaluated using the VIT technique by extrapolating the IFT to zero [51]. From simulation, the MMP is evaluated to be 12.07 MPa, which is very close to experimentally measured 11.76 MPa [40]. This demonstrates the reliable force fields adopted in this study for CO2/C10 systems.
646
T. Fang et al. / International Journal of Heat and Mass Transfer 141 (2019) 643–650
Fig. 4b, the dd/dP exhibits three regions. With increasing P, it first increases slightly, then rapidly, and finally approaches a constant. The swelling and dissolution in bulk phase are not subject to space constraints and different from that in nanoslits to be analyzed later. After more than 13 MPa, the thickness of the oil phase cannot be accurately measured due to the complete dissolution of the oil molecules. Therefore, the MMP is determined to be 11.38 MPa from the intersection of the latter two regions. The value from the DIM is in close agreement with the results evaluated by the VIT method and experimental result [34]. Then, it is further extended to determine the MMPs in nanoslits. 3.2. In nanoslits
Fig. 3. IFT between CO2 and C10 versus pressure. h: in this study; : exp. by Shaver et al. [40].
3.1.2. MMP determination by DIM Experimentally, the MMP can be also evaluated using the DIM technique by extrapolating dd/dP to zero, where dd/dP is the derivative of interfacial thickness d with respect to pressure P. Physically, dd/dP = 0 implies that the d becomes constant and does not vary with P; thus, a stable d is achieved upon mutually miscible between CO2 and oil [22]. Fig. 4 plots the d for CO2/C10 system in bulk phase, where the d was determined by the GDS (see Fig. 2). Over the entire pressure range examined, the d increases monotonically and no constant d is observed. The phenomenon is in remarkable contrast to experimental observation, because the bulk phase in this study is at a microscopic scale; thus there is significant mass transfer during CO2 diffusion into oil phase and a noticeable amount of C10 dissolved in CO2 after reaching a miscible state [23]. As show in
Fig. 4. (a) Interfacial thickness d (b) dd/dP versus pressure in bulk phase.
3.2.1. MMP determination Fig. 5a shows the d versus P in 5-nm nanoslit. Due to the constraint of pore surface, the thickness cannot continuously increase, which was observed as in bulk phase (Fig. 4a). Thus, the dd/dP first increases with increasing P, and approaches zero at a high pressure after passing a maximal value. As shown in Fig. 5b, by extrapolating the dd/dP zero, the MMP is evaluated to be around 8.33 MPa. Apparently, the MMP in the nanoslit is lower than that (11– 12 MPa) predicted in bulk phase, indicating that oil is more likely to dissolve into CO2 in a confined space. To quantify the difference, the potentials of mean force (PMFs) were calculated to estimate the barriers for C10 to dissolve into CO2
PMFðxÞ ¼ kB Tln
qðxÞ q0
ð7Þ
where q(x) is the density of C10 along the x-axis, q0 is the density of C10 in CO2 phase, kB is the Boltzmann constant. Fig. 6 compares the PMFs in bulk phase and 5-nm nanoslit, both at 6.29 MPa. It should be noted that the pressure in the nanoslit is for the CO2 phase (see Fig. S2). With the CO2 phase as a reference, the PMF of C10 in the oil phase is negative (i.e. lower). Thus, there exists a free energy barrier
Fig. 5. (a) d; (b) dd/dP versus pressure in 5-nm nanoslit.
647
T. Fang et al. / International Journal of Heat and Mass Transfer 141 (2019) 643–650
PMF (kcal/mol)
4
4
2
2
6.29 MPa
0
0
-2
-4
6.29 MPa
3.58 kcal/mol 0
50
100
-2
(a)
150
200
x (Å)
250
300
350
-4
1.11 kcal/mol 0
50
100
150
(b) 200
x (Å)
250
300
350
Fig. 6. Potentials of mean force of C10 upon CO2 flooding at 6.29 MPa (a) in bulk phase and (b) in 5-nm nanoslit.
DF for C10 to transfer from the oil to CO2 phase. In bulk phase, the DF is about 3.58 kcal/mol and higher than 1.11 kcal/mol in the nanoslit. This reveals in the nanoslit, C10 can more readily cross the interface and dissolve into CO2 phase; thus the MMP in the nanoslit is lower to achieve mutually miscible. This difference could be further elucidated by the competitive adsorption and confined diffusion of CO2 and C10.
3.2.2. Competitive adsorption Fig. 7 gives the simulation snapshots of oil swelling by CO2 from 0.1 ns to 5 ns in nanoslit (Fig. 7a) and bulk phase (Fig. 7b). The pressure considered is 6.29 MPa. At 0.1 ns, the interfaces of CO2/C10 in both systems are clear. At 2.5 ns, a considerable amount of CO2 in nanoslit is adsorbed on the solid surface, causing the oil phase to be detached and surrounded by CO2 molecules. However, in the bulk phase, the two phases still dissolved in a single direction near the interface. At 5 ns, more C10 molecules are observed to dissolve in CO2 in nanoslit, whereas in the bulk phase C10 molecules still gathered together. The difference in the distribution of molecules quantificationally is shown in Fig. 7c. Unlike the uniform distribution in the bulk phase, the density distribution curve exhibits obvious peaks. From the peak position, CO2 is found to have detached and enclosed the oil phase from the surface, thus proving that the effect of CO2 on the oil distribution is obvious.
It was found that the distribution of molecules in nanoslits is sensitive to the change in pressure. The density distributions of oil and CO2 in nanoslits at different pressures are shown in Fig. 8. The density of CO2 at a higher pressure fluctuates significantly, especially the adsorbed CO2 near the pore surface. It is observed that the adsorption peak increases obviously, and the second adsorption layer gradually appears, as shown in Fig. 8c. In the oil phase, the oil density at the center of the pore varies slightly, whereas the density of the adsorbed oil decreases and the peak position shifts to the center of pores, thus leading to a significant desorption. Besides static structures, the dynamics of fluid molecules is also heavily affected by confinement. The mean-squared displacements (MSDs) of CO2 and C10 in nanoslits are shown in Fig. S3, revealing that the confined effect of nanoslits on molecular diffusion. Since such a large density change is difficult to observe in the bulk phase, the interface is in a uniform state for a long time without large fluctuations, which is one of the critical factors leading to a high energy barrier at the CO2/C10 interface.
3.2.3. Confinement effects on miscible process in nanoslits For the nanoslits with different widths (2, 3, 5, 9, 15 nm), Fig. S4 shows the d and dd/dP versus pressure in each nanoslit. Fig. 9a lists the results in different nanoslits and the MMP is plotted versus pore size. A low MMP value is considered to represent a better miscible state. The miscibility seems easier to achieve when the pore
Fig. 7. CO2/C10 system in (a) 5-nm nanoslit and (b) bulk phase at 6.29 MPa from 0.1 to 5 ns. CO2: red, C10: blue; (c) Densities of CO2 and C10 along z-axis at 5 ns.
648
T. Fang et al. / International Journal of Heat and Mass Transfer 141 (2019) 643–650
Fig. 8. Density evolution in 5-nm nanoslit at different pressure. (a) CO2 and (b) oil in nanoslit (c) CO2 and (d) oil near pore surface.
Fig. 9. (a) MMP versus pore size (b) Number density of C10 at 6.29 MPa (c) Free energy barrier DF of C10 in CO2 flooding versus pore size.
size is reduced because of stronger confinement effect. The trend of MMP value is in agreement with the published results by theoretical method [2,22,24]. In addition, it is noteworthy that in Fig. 9a, the MMP value starts to increase once the pore radius is smaller than 5-nm. It is because that the proportion of the adsorbed C10 increases with reducing pore size, while C10 molecules in the free state decrease. The size of C10 is schematically illustrated as an insert in Fig. 9b. The diameter of C10 is around 1.32 nm, so C10 could diffuse freely in a nanoslit with width larger than 5-nm. Obviously, a necessary condition for the fluid miscibility development is that both liquid and gas molecules can diffuse freely in the porous medium [2]. The above analysis proves that if the effect of the molecular adsorption layer is considered, the optimum pore width value will increase, and the stability of the adsorption layer will make the molecular migration more uncertain. In the latter study, the specific effects of the adsorption layer and different crude oil components on the miscible behavior need further analysis.
The free energy barrier DF for C10 in CO2 flooding is shown in Fig. 9c, presenting a similar pattern with the MMP of mixing. To be specific, the DF decreases with the reduction of pore size to 5nm, and it starts to increase afterwards. On the whole, a smaller DF leads to easier miscibility. Thus, the confinement effect facilitates the miscibility of CO2/oil. However, there is a bottom limit in the pore size reduction, which is a necessary condition for the occurrence of miscibility [2]. When the free-molecule area is larger than a single C10 molecule, the miscible state could be achieved. Otherwise, the miscible state would be hard to reach. 4. Conclusions In this study, molecular dynamic simulations were performed to study the miscible mechanism of CO2/oil phase in nanoslits. It has important economic benefits as the generation of all such properties using a rigorously experimental approach is formidable. This is especially true for CO2-EOR process involving the molecular
T. Fang et al. / International Journal of Heat and Mass Transfer 141 (2019) 643–650
microscopic interaction behavior at a nanoscale. In this study, first, the interface thickness of the oil phase by introducing the Gibbs dividing surface was determined. The method helps address the problem that the MMP value is hard to quantify through molecular simulation. In addition, the feasibility of DIM was also confirmed and extended to the calculation of MMP in nanoslits, thus helping to analyze the effect of molecular dynamics on MMP. Second, the results showed that the MMPs in nanoslits were significantly smaller than those in bulk phase, suggesting that oil and gas phases were more likely to attain a miscible state. To explain the phenomenon, the PMFs of C10 in CO2 flooding and competitive adsorption behavior of fluid molecules were analyzed, providing a microscopic explanation for the miscible mechanism in confined spaces. Rising the probability of interaction between C10 and CO2 is a vital factor affecting the miscible phase. Lastly, based on the effect of nanoslits on molecules, five confined nanoslits with the pore size ranging from 2 to 15 nm were examined. The simulation results revealed that the miscible state was easier to achieve with the reduction of pore size due to the confinement effect. Besides, it was found that the width of nanoslit with the lowest MMP was larger than the reported value due to the consideration of adsorption layer, thus making the molecular migration more uncertain. Our work provided the understanding of the miscible mechanism on a molecular level, and the results have significant promise for the oil exploitation and might show some promising applications in the future. Declaration of Competing Interest The authors declared that there is no conflict of interest. Acknowledgements This work is financially supported by the National Natural Science Foundation of China (51874332, U1663206, 51622405), Climb Taishan Scholar Program in Shandong Province (tspd20161004) and Fundamental Research Funds for the Central Universities (YCX2017067, 18CX06006A). Appendix A. Supplementary material Supplementary data to this article can be found online at https://doi.org/10.1016/j.ijheatmasstransfer.2019.06.083. References [1] S. Mondal, S. De, CO2 based power cycle with multi-stage compression and intercooling for low temperature waste heat recovery, Energy 90 (2015) 1132– 1143. [2] K. Zhang, N. Jia, S. Li, L. Liu, Thermodynamic phase behaviour and miscibility of confined fluids in nanopores, Chem. Eng. J. 351 (2018) 1115–1128. [3] W. Yu, H.R. Lashgari, K. Wu, CO2 injection for enhanced oil recovery in Bakken tight oil reservoirs, Fuel 159 (2015) 354–363. [4] M.G. Rezk, J. Foroozesh, Determination of mass transfer parameters and swelling factor of CO2-oil systems at high pressures, Int. J. Heat Mass Transf. 126 (2018) 380–390. [5] T. Lee, L. Bocquet, B. Coasne, Activated desorption at heterogeneous interfaces and long-time kinetics of hydrocarbon recovery from nanoporous media, Nat. Commun. 7 (2016). [6] B. Wei, L. Lu, W. Pu, R. Wu, X. Zhang, Y. Li, F. Jin, Production dynamics of CO2 cyclic injection and CO2 sequestration in tight porous media of Lucaogou formation in Jimsar sag, J. Petrol. Sci. Eng. 157 (2017) 1084–1094. [7] P. Gao, B.F. Towler, G. Pan, Strategies for evaluation of the CO2 miscible flooding process, Abu Dhabi International Petroleum Exhibition and Conference, Society of Petroleum Engineers, 2010. [8] K. Zhang, B. Kong, J. Zhan, R. He, T. Qin, K. Wu, S. Chen, J. Zhang, Effects of nanoscale pore confinement on CO2 immiscible and miscible processes, SPE Low Perm Symposium, Society of Petroleum Engineers, 2016. [9] P. Nguyen, D. Mohaddes, J. Riordon, H. Fadaei, P. Lele, D. Sinton, Fast fluorescence-based microfluidic method for measuring minimum miscibility pressure of CO2 in crude oils, Anal. Chem. 87 (2015) 3160–3164.
649
[10] K. Zhang, N. Jia, F. Zeng, Application of predicted bubble-rising velocities for estimating the minimum miscibility pressures of the light crude oil-CO2 systems with the rising bubble apparatus, Fuel 220 (2018) 412–419. [11] M. Escrochi, N. Mehranbod, S. Ayatollahi, The gas-oil interfacial behavior during gas injection into an asphaltenic oil reservoir, J. Chem. Eng. Data 58 (2013) 2513–2526. [12] O.N. Amézquita, S. Enders, P.T. Jaeger, R. Eggers, Interfacial properties of mixtures containing supercritical gases, J. Supercrit. Fluid 55 (2010) 724–734. [13] A. Mejía, M. Cartes, H. Segura, E.A. Müller, Use of equations of state and coarse grained simulations to complement experiments: describing the interfacial properties of carbon dioxide + decane and carbon dioxide + eicosane mixtures, J. Chem. Eng. Data 59 (2014) 2928–2941. [14] D.N. Rao, J.I. Lee, Application of the new vanishing interfacial tension technique to evaluate miscibility conditions for the Terra Nova Offshore Project, J. Petrol. Sci. Eng. 35 (2002) 247–262. [15] R. Wang, F. Peng, K. Song, G. Feng, Z. Guo, Molecular dynamics study of interfacial properties in CO2 enhanced oil recovery, Fluid Phase Equilibr. 467 (2018) 25–32. [16] T. Babadagli, Development of mature oil fields-a review, J. Petrol. Sci. Eng. 57 (2017) 221–246. [17] C.R. Clarkson, N. Solano, R.M. Bustin, A.M.M. Bustin, G.R.L. Chalmers, L. He, Y.B. Melnichenko, A.P. Radlin´ski, T.P. Blach, Pore structure characterization of North American shale gas reservoirs using USANS/SANS, gas adsorption, and mercury intrusion, Fuel 103 (2013) 606–616. [18] L. Bocquet, E. Charlaix, Nanofluidics, from bulk to interfaces, Chem. Soc. Rev. 39 (2010) 1073–1095. [19] S. Wang, Q. Feng, F. Javadpour, Q. Hu, K. Wu, Competitive adsorption of methane and ethane in montmorillonite nanopores of shale at supercritical conditions: a grand canonical Monte Carlo simulation study, Chem. Eng. J. 355 (2019) 76–90. [20] S.B. Hawthorne, D.J. Miller, L. Jin, C.D. Gorecki, Rapid and simple capillaryrise/vanishing interfacial tension method to determine crude oil minimum miscibility pressure: pure and mixed CO2, methane, and ethane, Energ. Fuel 30 (2016) 6365–6372. [21] A. Fazlali, M. Nikookar, A. Agha-Aminiha, A.H. Mohammadi, Prediction of minimum miscibility pressure in oil reservoirs using a modified SAFT equation of state, Fuel 108 (2013) 675–681. [22] K. Zhang, N. Jia, F. Zeng, P. Luo, A New diminishing interface method for determining the minimum miscibility pressures of light oil-CO2 systems in bulk phase and nanopores, Energ. Fuel 31 (2017) 12021–12034. [23] A. Sharbatian, A. Abedini, Z. Qi, D. Sinton, Full characterization of CO2-oil properties on-chip: solubility, diffusivity, extraction pressure, miscibility, and contact angle, Anal. Chem. 90 (2018) 2461–2467. [24] K. Zhang, S. Seetahal, D. Alexander, J. Lv, Y. Hu, X. Lu, D. Zhang, K. Wu, Z. Chen, Effect of confinement on gas and oil relative permeability during CO2 flooding in tight oil reservoirs, SPE Trinidad and Tobago Section Energy Resources Conference, Society of Petroleum Engineers, 2016. [25] Y. Zhang, H.R. Lashgari, Y. Di, K. Sepehrnoori, Capillary pressure effect on phase behavior of CO2/hydrocarbons in unconventional reservoirs, Fuel 197 (2017) 575–582. [26] H. Wang, X. Wang, X. Jin, D. Cao, Molecular dynamics simulation of diffusion of shale oils in montmorillonite, J. Phys. Chem. C 120 (2016) 8986–8991. [27] Y. Yan, Z. Dong, Y. Zhang, P. Wang, T. Fang, J. Zhang, CO2 activating hydrocarbon transport across nanopore throat: insights from molecular dynamics simulation, Phys. Chem. Chem. Phys. 19 (2017) 30439–30444. [28] T. Le, A. Striolo, D.R. Cole, Propane simulated in silica pores: Adsorption isotherms, molecular structure, and mobility, Chem. Eng. Sci. 121 (2015) 292– 299. [29] T. Le, A. Striolo, D.R. Cole, CO2-C4H10 mixtures simulated in silica slit pores: relation between structure and dynamics, J. Phys. Chem. C 119 (2015) 15274– 15284. [30] T. Fang, M. Wang, J. Li, B. Liu, Y. Shen, Y. Yan, J. Zhang, Study on the asphaltene precipitation in CO2 flooding: a perspective from molecular dynamics simulation, Ind. Eng. Chem. Res. 57 (2018) 1071–1077. [31] D. Hou, Z. Lu, X. Li, H. Ma, Z. Li, Reactive molecular dynamics and experimental study of graphene-cement composites: Structure, dynamics and reinforcement mechanisms, Carbon 115 (2017) 188–208. [32] D. Hou, T. Li, P. Wang, Molecular dynamics study on the structure and dynamics of NaCl solution transport in the nanometer channel of CASH gel, ACS Sustain. Chem. Eng. 6 (2018) 9498–9509. [33] Material Studio, Accelrys Software Inc.: San Diego, CA, 2010. [34] Z. Yang, M. Li, B. Peng, M. Lin, Z. Dong, Dispersion property of CO2 in oil. 1. Volume expansion of CO2 + alkane at near critical and supercritical condition of CO2, J. Chem. Eng. Data 57 (2012) 882–889. [35] Y. Yan, C. Li, Z. Dong, T. Fang, B. Sun, J. Zhang, Enhanced oil recovery mechanism of CO2 water-alternating-gas injection in silica nanochannel, Fuel 190 (2017) 253–259. [36] T. Dong, N.B. Harris, Pore size distribution and morphology in the horn river shale, middle and upper Devonian, northeastern British Columbia, Canada. AAPG Mem. 102 (2013) 67–79. [37] Z.Z. Li, T. Min, Q. Kang, Y.L. He, W.Q. Tao, Investigation of methane adsorption and its effect on gas transport in shale matrix through microscale and mesoscale simulations, Int. J. Heat. Mass. Transf. 98 (2016) 675–686. [38] S. Wang, Q. Feng, M. Zha, F. Javadpour, Q. Hu, Supercritical methane diffusion in shale nanopores: effects of pressure, mineral types, and moisture content, Energ. Fuel 32 (2017) 169–180.
650
T. Fang et al. / International Journal of Heat and Mass Transfer 141 (2019) 643–650
[39] S. Plimpton, Fast parallel algorithms for short-range molecular dynamics, J. Comput. Phys. 117 (1995) 1–19. [40] R.D. Shaver, R.L. Robinson Jr, K.A.M. Gasem, An automated apparatus for equilibrium phase compositions, densities, and interfacial tensions: data for carbon dioxide+ decane, Fluid Phase Equilibr. 179 (2001) 43–66. [41] S. Nosé, A molecular dynamics method for simulations in the canonical ensemble, Mol. Phys. 52 (1984) 255–268. [42] W. Humphrey, A. Dalke, K. Schulten, VMD: visual molecular dynamics, J. Mol. Graph. 14 (1996) 33–38. [43] R.T. Cygan, J.J. Liang, A.G. Kalinichev, Molecular models of hydroxide, oxyhydroxide, and clay phases and the development of a general force field, J. Phys. Chem. B 108 (2004) 1255–1266. [44] S. Iglauer, M.S. Mathew, F. Bresme, Molecular dynamics computations of brine-CO2 interfacial tensions and brine-CO2-quartz contact angles and their effects on structural and residual trapping mechanisms in carbon geosequestration, J. Colloid. Interf. Sci. 386 (2012) 405–414. [45] A.P. Thompson, S.J. Plimpton, W. Mattson, General formulation of pressure and stress tensor for arbitrary many-body interaction potentials under periodic boundary conditions, J. Chem. Phys. 131 (2009) 154107.
[46] Q. Xu, J. Jiang, Computational characterization of ultrathin polymer membranes in liquids, Macromolecules 51 (2018) 7169–7177. [47] L. Bocquet, J.L. Barrat, Hydrodynamic boundary conditions, correlation functions, and Kubo relations for confined fluids, Phys. Rev. E 49 (1994) 3079. [48] I.C. Bourg, C.I. Steefel, Molecular dynamics simulations of water structure and diffusion in silica nanopores, J. Phys. Chem. B 116 (2012) 11556– 11564. [49] T. Fang, M. Wang, Y. Gao, Y. Zhang, Y. Yan, J. Zhang, Enhanced oil recovery with CO2/N2 slug in low permeability reservoir: Molecular dynamics simulation, Chem. Eng. Sci. 197 (2019) 204–211. [50] S. Wang, F. Javadpour, Q. Feng, Molecular dynamics simulations of oil transport through inorganic nanopores in shale, Fuel 171 (2016) 74–86. [51] S. Li, K. Zhang, N. Jia, L. Liu, Evaluation of four CO2 injection schemes for unlocking oils from low-permeability formations under immiscible conditions, Fuel 234 (2018) 814–823.
Contents lists available at ScienceDirect
International Journal of Heat and Mass Transfer journal homepage: www.elsevier.com/locate/ijhmt
Molecular insight into the miscible mechanism of CO2/C10 in bulk phase and nanoslits Timing Fang a,b, Yingnan Zhang a, Jie Liu b, Bin Ding c, Youguo Yan a,⇑, Jun Zhang a,⇑ a
School of Materials Science and Engineering, China University of Petroleum, Qingdao, Shandong 266580, China Department of Chemical and Biomolecular Engineering, National University of Singapore, 117576 Singapore, Singapore c Research Institute of Petroleum Exploration & Development of PetroChina, Beijing 100083, China b
a r t i c l e
i n f o
Article history: Received 16 April 2019 Received in revised form 9 June 2019 Accepted 24 June 2019 Available online 9 July 2019 Keywords: Molecular dynamics simulation Carbon dioxide Miscibility Nanoslit
a b s t r a c t The determination of minimum miscibility pressure (MMP) of CO2/oil in tight oil reservoirs has aroused huge attention from researchers. In this study, we characterized the mass transfer process and determined the MMP of CO2/oil systems in nanoslits through molecular simulation. A bulk phase and five confined nanoslits with the pore size ranging from 2 to 15 nm were examined. It was found that the MMPs in nanoslits were significantly lower than those in bulk phase, thus revealing that rising the probability of molecular interaction is a vital factor affecting the miscible phase. Furthermore, due to the consideration of adsorption layer, the width of nanoslit with the lowest MMP was larger than the reported value, which will bring more uncertainty to the molecular migration. This simulation study highlighted the crucial role of the confinement effect in the miscible process of CO2 flooding, which might facilitate the development of low permeability reservoirs. Ó 2019 Elsevier Ltd. All rights reserved.
1. Introduction Both experimental and theoretical studies proved CO2 flooding an effective method for enhanced oil recovery (EOR), which is likely to enhance oil recovery by nearly 8–14% of the original oil in place (OOIP) through an efficient mass transfer [1–4]. For the tight and shale oil reservoirs with massive micropores, CO2 flooding plays a significant role in EOR because CO2 can easily transfer into the oil phase [5,6]. The minimum miscibility pressure (MMP) is defined as the lowest operating pressure, at which oil and CO2 are mutually miscible in the reservoir. Under a higher pressure over the MMP, oil recovery will be enhanced by CO2 injection. If the reservoir pressure is below the MMP, oil and CO2 will be incompletely miscible, thereby leading to a low oil recovery [7,8]. Accordingly, the accurate determination of the MMP in a CO2/oil system is critical for EOR. In the literature, several experimental techniques and theoretical models have been reported to determine the MMPs of crude oil/solvent systems. The slim-tube tests [9], the rising bubble apparatus [10] and the vanishing interfacial tension (VIT) technique [11] are considered the commonest experimental technique. Represented by the VIT, the method incorporates various research
⇑ Corresponding authors. E-mail address: [email protected] (J. Zhang). https://doi.org/10.1016/j.ijheatmasstransfer.2019.06.083 0017-9310/Ó 2019 Elsevier Ltd. All rights reserved.
methods (e.g. theoretical analysis [12] and molecular simulation [13]) to make the measurement object wider in range and more precise. In the VIT technique, the MMP is evaluated by linearly extrapolating the interfacial tension (IFT) at CO2/oil interface to zero, i.e., the miscibility between CO2 and oil is achieved at the zero-IFT condition [14]. Besides, by combining molecular dynamics (MD) simulations and the equation of state, Andrés et al. [13] and Wang et al. [15] estimated the MMPs and then explored the IFTs between injected gas and crude oil in bulk phase. The MMPs evaluated using the VIT technique were found close to those using slim-tube test. However, it has been estimated that nearly 40% of OOIP exists in unconventional oil reservoirs, where oil is adsorbed in nanosized pores [16]. Accordingly, it is necessary to investigate the MMPs of CO2/oil systems in nanopores. However, the exact MMPs should be experimentally measure in nanopores with sizes ranging over 2–100 nm [17]. Due to the confinement and surface interaction, the physical properties of fluids in nanopores significantly differ from those of conventional macroscopic systems [18,19]. These methods are primarily used to determine MMP in bulk rather than MMP in nanopores [20,21]. In recent years, the behavior of miscible in nanopores has gradually aroused huge attention. Considering the particularity of the phase change of the oil phase in confined spaces, instead of direct measurement, the researchers used theoretic modelling and deduction. The methods developed based on the physical and structural properties of oil and gas have been employed. For instance, the
644
T. Fang et al. / International Journal of Heat and Mass Transfer 141 (2019) 643–650
diminishing interface method (DIM) proposed by Zhang et al. [22] addressed the problem that the MMP of heterogeneous interface is difficult to measure accurately. The DIM determines the relationship between the variation of thickness and the pressure, based on the stable interface thickness of miscible phase, to measure the specific MMP value. Similar studies have been published on the proof and extension of the method [8,23]. Besides, from the millimeter to the nanometer scale, the researchers adopted the Peng-Robinson equation of state (PR-EOS) to measure MMPs in nanopores through the calculation of solubility parameters and capillary pressure. The critical properties dominating the phase behavior and miscibility of confined fluids were further revealed [24,25]. The above studies laid the solid foundation of theoretical models for further studies. To explore the microscopic mechanisms not readily available in experimental and theoretical methods, MD simulation was performed to study microstructural variations and dynamics features of oil molecules in nanopores [26,27]. To explore the effects of temperature and oil density on the diffusion coefficients of alkane in nanopores, Cao et al. utilized MD simulations [26]. Likewise, Wang et al. characterized the fluid flow of CO2/octane in nanopores affected by driving force and pore size [19]. Thu et al. examined the relation between interaction of CO2/oil mixtures and dynamics of oil molecules in silica pores. Their results demonstrated the effect of molecular interactions on the adsorption and migration behavior of alkanes, thus indicating the implication of microscopic analysis [28,29]. Note that molecular simulation is a more intuitive method to examine sophisticated interactions and unravel microscopic properties [30–32]. However, on molecular and atomic level, fewer studies have focused on the microscopic mechanism of miscible behavior in nanopores, which helps enhance utilization efficiency of CO2 in tight oil exploitation. Based on MD simulations, interfacial thickness method was employed in this study to characterize the miscibility of CO2/decane (C10) systems in nanoslits and then to evaluate the MMPs. To verify the methodology, CO2/C10 systems in bulk phase were also examined. In Section 2, the molecular models of CO2/C10 systems in bulk phase and nanoslit were described, followed by simulation methods. In Section 3, first, the conventional VIT method and the interfacial thickness method were used to estimate the MMPs of CO2/C10 systems in bulk phase. Next, the interfacial thickness method was extended to evaluate the MMPs in the nanoslits. Subsequently, the miscible dynamics in different nanoslits with widths of 2, 3, 5, 9 and 15 nm was examined, respectively. Lastly, the confinement effects on miscible behavior is discussed. In Section 4, the concluding remarks are summarized. 2. Methods 2.1. Models and simulations 2.1.1. In bulk phase All the systems were constructed Materials Studio (MS) [33]. Fig. 1a shows a typical CO2/C10 system in bulk phase. Initially,
Fig. 1. A CO2/C10 system in (a) bulk phase & (b) a nanoslit. CO2: red, C10: blue. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
the central region of the system contained solely C10 molecules, which were added into a simulation cell with random orientations; then, CO2 molecules were added on each side of C10. For each CO2/C10 system, the number of C10 was set to 300 and the number of CO2 molecules was calculated by nCO2 = xCO2nt, where xCO2 is the solubility of CO2 in C10 (in terms of mole fraction) and nt is the total number of molecules. To ensure sufficient CO2 molecules interacting with C10, the experimental xCO2 at 12.62 MPa was selected [34] and then 300 more CO2 molecules were added into the system. The thickness of oil film was set as 5 nm. 2.1.2. In nanoslits To examine the confinement effects, the CO2/C10 systems were placed into a nanoslit as shown in Fig. 1b. The nanoslit was represented by two hydroxylated silica surfaces, each was cleaved from a-quartz along the (1 0 0) plane with dimensions (xyz) of 39.98 nm 2.70 nm 1.507 nm. As the nanopore size in oil reservoirs normally ranges from 3 to 8 nm [35,36], most of the confined CO2/C10 systems were simulated in a nanoslit with width of 5-nm along the z-axis in this study. Nevertheless, other nanoslits (with width of 2, 3, 9 and 15 nm) were also considered to examine the effect of pore size. The initial C10 phase with dimensions of 8.50 nm 2.70 nm 5.00 nm consisted of 274 C10 molecules. 2.1.3. Determination of pressure in nanoslits To estimate the pressure inside the nanoslit, the system was divided into small bins each of 0.5 Å along the z-axis. The pressure within a bin was calculated by
Pbin ¼
Sxx þ Syy þ Szz 3V
ð1Þ
The symmetric stress tensor of atom i, Sab, was given by
2
Sab
3
Np P
ðr 1a F 1b þ r 2a F 2b Þ 7 6 mv a v b þ 7 6 n¼1 7 6 7 6 N b 7 6 1 P ðr1a F 1b þ r2a F 2b Þ 7 6þ2 7 6 n¼1 7 6 7 6 Na 7 6 1 P ðr1a F 1b þ r2a F 2b þ r3a F 3b Þ 7 6þ3 7 6 n¼1 7 ¼ 6 7 6 Nd 7 6 1 P ðr1a F 1b þ r2a F 2b þ r3a F 3b þ r4a F 4b Þ 7 6þ4 7 6 n¼1 7 6 7 6 Ni P 7 6 1 6þ4 ðr1a F 1b þ r2a F 2b þ r3a F 3b þ r4a F 4b Þ 7 7 6 n¼1 7 6 7 6 N p P 5 4 r ia F ib þKspaceðr ia ; F ib Þ þ 1 2
ð2Þ
n¼1
where a and b take values of x, y, z to produce six components of a symmetric tensor. In this equation, the 1st term on the right-hand side is the kinetic energy contribution from atom i; the 2nd term is a pairwise energy contribution, where n iterates on the Np neighbors of atom i, r1 and r2 are the positions of the two atoms, and F1 and F2 are the forces on the two atoms resulting from the pairwise interaction; the 3rd to 8th terms are contributions of the N b bonds, Na angle, Nd dihedral, Ni improper, long-range Coulombic interaction (Kspace term) and internal constraint force, respectively. This method has been widely employed to estimate the pressure in a confined system [37,38]. To demonstrate and compare with available simulation results in the literature, a system consisting of CH4 molecules in a graphitic nanoslit with width of 3.6 nm was constructed and simulated. Fig. S1 shows the density and pressure profiles of CH4 in the nanoslit along the z-axis (perpendicular to the pore surface). In the middle of the nanoslit, both density and pressure remain as a constant, representing ‘‘free gas” with negligible interaction with the pore
645
T. Fang et al. / International Journal of Heat and Mass Transfer 141 (2019) 643–650
surface. As listed in Table S1, the estimated density and pressure of free CH4 match very well with a reported simulation study [37]. For the current study, as shown in Fig. S2, a region of constant pressure is also observed for CO2 in a nanoslit. 2.1.4. Simulation details MD simulations were performed by large-scale atomic/molecular massively parallel simulator (LAMMPS) software ((Version 16 Mar 2018)) [39]. Based on practical oil reservoir conditions [40], the temperature was kept at 344.15 K and controlled by a NoséHoover thermostat [41]. The cutoff distance to evaluate the LJ interactions was 1.0 nm and the Ewald summation method was used for the electrostatic interactions. The periodic boundary conditions were applied in all the directions. In this work, we carried out all-atom MD simulation of each system, the simulation duration was 10 ns. The time step was 1 fs and every 1000 steps were saved for data analysis. After the simulation, all the snapshots were displayed by VMD software [42]. Bulk phase: the simulations of CO2/C10 systems in bulk phase were performed in a NPT ensemble (with constant number of molecules, pressure and temperature). Initially, the dimensions of each system were set to be 50 Å in the x and y directions. Nanoslits: the hydroxylated silica surfaces were described by the CLAYFF force field [43]. The simulations were conducted in a NVT ensemble (with constant number of molecules, volume and temperature). Identical to those in bulk phase. To obtain an initial equilibrated C10 phase, 5 ns simulation was first performed. Then, CO2 phase with different densities, were placed on each side of C10 phase for 10 ns simulation. The non-bonding interaction consists of the short-range van der Waals (vdW) and long-range electrostatic interaction. The vdW is represented by Lennard-Jones (LJ) potential and the electrostatic interaction is represented by Coulombic potential,
"
V non-bonded ¼ 4eij
rij
12
r ij
rij
rij
6 #
þ
qi qj 4pe0 r ij
ð3Þ
where eij and rij are the well depth and collision diameter, rij is the distance between atoms i and j, qi is the atomic charge of atom i, and e0 = 8.8542 1012 C2 N1 m2 is the permittivity of vacuum. For the bonding interaction, there are three components including bond stretching, angle bending and dihedral torsion. The specific parameters of force field are listed in Section S2 in the Supporting Information. 2.2. Interfacial tensions (IFTs) The interfacial tension was calculated by the formulation of the Gibbs interfacial tension in LAMMPS. For CO2/C10 systems in bulk phase, two interfaces are both perpendicular to the x-axis and parallel to the zy-plane, and the IFTs between CO2 and C10 phase were calculated by the expression of pressure tensor [44,45]:
1 c¼ 2
ZLx
Pyy þ Pzz 1 Lx Pxx ðPN ðxÞ PT ðxÞÞdx ¼ 2 2
ð4Þ
where PN and PT are the normal and tangential pressures, respectively; Paa (a = x, y, z) is the diagonal element of the pressure tensor; and Lx is the box length in x-axis. The pressure components were evaluated by the virial equation: N X i¼1
mi v i;a v i;b þ
N1 X N X i¼1
rij;a f ij;b
box. Finally, the average of IFTs during the last 2 ns was taken as the final value. 2.3. Interfacial thicknesses In both bulk phase and nanoslits, there is no distinct interface between CO2 and oil phase; moreover, the density of oil is not uniform along the x-axis. To estimate the CO2/oil interfacial thickness, a recently proposed characterization method for polymer/solvent systems [46] is adopted here. As illustrated in Fig. 2, based on the concept of Gibbs dividing surface (GDS), the density profile of C10 along the x-axis can be fitted by [46]:
1 2
qðxÞ ¼ q0 erf
x x1 pffiffiffi r 2
erf
x x2 pffiffiffi r 2
ð6Þ
where q0 is the bulk density away from the CO2/C10 interface, x1 and x2 represent the positions of Gibbs surfaces; r represents the standard deviation in length. d is the length (x2-x1) between two Gibbs surfaces. As the Gibbs surface is an idealized boundary separating two phases, d can be considered as the effective thickness of oil phase in the system. Similar methods were applied to determine the locations of fluid/solid interfaces [47–50]
3. Results and discussion 3.1. In bulk phase
0
Pab V ¼
Fig. 2. Density profile of C10 along the x-axis. The equilibrium snapshot for CO2/C10 in bulk phase is shown on the top.
ð5Þ
j>i
where vi, a the velocity of atom i in direction a; rij = ri rj; fij the force exerted on atom i by atom j; and V is the volume of simulation
3.1.1. MMP determination by VIT Fig. 3 shows the IFT between CO2 and oil versus pressure in bulk phase at 344.15 K. The predicted IFT in this study matches well with the experimental study by Shaver et al. [40]. With rising pressure, the IFT monotonically decreases and has a linear relationship with pressure. Through linear regression, we find IFT = 19.61– 1.62 P from our simulation and IFT = 19.14–1.63 P from the experiment [40]. The MMP can be evaluated using the VIT technique by extrapolating the IFT to zero [51]. From simulation, the MMP is evaluated to be 12.07 MPa, which is very close to experimentally measured 11.76 MPa [40]. This demonstrates the reliable force fields adopted in this study for CO2/C10 systems.
646
T. Fang et al. / International Journal of Heat and Mass Transfer 141 (2019) 643–650
Fig. 4b, the dd/dP exhibits three regions. With increasing P, it first increases slightly, then rapidly, and finally approaches a constant. The swelling and dissolution in bulk phase are not subject to space constraints and different from that in nanoslits to be analyzed later. After more than 13 MPa, the thickness of the oil phase cannot be accurately measured due to the complete dissolution of the oil molecules. Therefore, the MMP is determined to be 11.38 MPa from the intersection of the latter two regions. The value from the DIM is in close agreement with the results evaluated by the VIT method and experimental result [34]. Then, it is further extended to determine the MMPs in nanoslits. 3.2. In nanoslits
Fig. 3. IFT between CO2 and C10 versus pressure. h: in this study; : exp. by Shaver et al. [40].
3.1.2. MMP determination by DIM Experimentally, the MMP can be also evaluated using the DIM technique by extrapolating dd/dP to zero, where dd/dP is the derivative of interfacial thickness d with respect to pressure P. Physically, dd/dP = 0 implies that the d becomes constant and does not vary with P; thus, a stable d is achieved upon mutually miscible between CO2 and oil [22]. Fig. 4 plots the d for CO2/C10 system in bulk phase, where the d was determined by the GDS (see Fig. 2). Over the entire pressure range examined, the d increases monotonically and no constant d is observed. The phenomenon is in remarkable contrast to experimental observation, because the bulk phase in this study is at a microscopic scale; thus there is significant mass transfer during CO2 diffusion into oil phase and a noticeable amount of C10 dissolved in CO2 after reaching a miscible state [23]. As show in
Fig. 4. (a) Interfacial thickness d (b) dd/dP versus pressure in bulk phase.
3.2.1. MMP determination Fig. 5a shows the d versus P in 5-nm nanoslit. Due to the constraint of pore surface, the thickness cannot continuously increase, which was observed as in bulk phase (Fig. 4a). Thus, the dd/dP first increases with increasing P, and approaches zero at a high pressure after passing a maximal value. As shown in Fig. 5b, by extrapolating the dd/dP zero, the MMP is evaluated to be around 8.33 MPa. Apparently, the MMP in the nanoslit is lower than that (11– 12 MPa) predicted in bulk phase, indicating that oil is more likely to dissolve into CO2 in a confined space. To quantify the difference, the potentials of mean force (PMFs) were calculated to estimate the barriers for C10 to dissolve into CO2
PMFðxÞ ¼ kB Tln
qðxÞ q0
ð7Þ
where q(x) is the density of C10 along the x-axis, q0 is the density of C10 in CO2 phase, kB is the Boltzmann constant. Fig. 6 compares the PMFs in bulk phase and 5-nm nanoslit, both at 6.29 MPa. It should be noted that the pressure in the nanoslit is for the CO2 phase (see Fig. S2). With the CO2 phase as a reference, the PMF of C10 in the oil phase is negative (i.e. lower). Thus, there exists a free energy barrier
Fig. 5. (a) d; (b) dd/dP versus pressure in 5-nm nanoslit.
647
T. Fang et al. / International Journal of Heat and Mass Transfer 141 (2019) 643–650
PMF (kcal/mol)
4
4
2
2
6.29 MPa
0
0
-2
-4
6.29 MPa
3.58 kcal/mol 0
50
100
-2
(a)
150
200
x (Å)
250
300
350
-4
1.11 kcal/mol 0
50
100
150
(b) 200
x (Å)
250
300
350
Fig. 6. Potentials of mean force of C10 upon CO2 flooding at 6.29 MPa (a) in bulk phase and (b) in 5-nm nanoslit.
DF for C10 to transfer from the oil to CO2 phase. In bulk phase, the DF is about 3.58 kcal/mol and higher than 1.11 kcal/mol in the nanoslit. This reveals in the nanoslit, C10 can more readily cross the interface and dissolve into CO2 phase; thus the MMP in the nanoslit is lower to achieve mutually miscible. This difference could be further elucidated by the competitive adsorption and confined diffusion of CO2 and C10.
3.2.2. Competitive adsorption Fig. 7 gives the simulation snapshots of oil swelling by CO2 from 0.1 ns to 5 ns in nanoslit (Fig. 7a) and bulk phase (Fig. 7b). The pressure considered is 6.29 MPa. At 0.1 ns, the interfaces of CO2/C10 in both systems are clear. At 2.5 ns, a considerable amount of CO2 in nanoslit is adsorbed on the solid surface, causing the oil phase to be detached and surrounded by CO2 molecules. However, in the bulk phase, the two phases still dissolved in a single direction near the interface. At 5 ns, more C10 molecules are observed to dissolve in CO2 in nanoslit, whereas in the bulk phase C10 molecules still gathered together. The difference in the distribution of molecules quantificationally is shown in Fig. 7c. Unlike the uniform distribution in the bulk phase, the density distribution curve exhibits obvious peaks. From the peak position, CO2 is found to have detached and enclosed the oil phase from the surface, thus proving that the effect of CO2 on the oil distribution is obvious.
It was found that the distribution of molecules in nanoslits is sensitive to the change in pressure. The density distributions of oil and CO2 in nanoslits at different pressures are shown in Fig. 8. The density of CO2 at a higher pressure fluctuates significantly, especially the adsorbed CO2 near the pore surface. It is observed that the adsorption peak increases obviously, and the second adsorption layer gradually appears, as shown in Fig. 8c. In the oil phase, the oil density at the center of the pore varies slightly, whereas the density of the adsorbed oil decreases and the peak position shifts to the center of pores, thus leading to a significant desorption. Besides static structures, the dynamics of fluid molecules is also heavily affected by confinement. The mean-squared displacements (MSDs) of CO2 and C10 in nanoslits are shown in Fig. S3, revealing that the confined effect of nanoslits on molecular diffusion. Since such a large density change is difficult to observe in the bulk phase, the interface is in a uniform state for a long time without large fluctuations, which is one of the critical factors leading to a high energy barrier at the CO2/C10 interface.
3.2.3. Confinement effects on miscible process in nanoslits For the nanoslits with different widths (2, 3, 5, 9, 15 nm), Fig. S4 shows the d and dd/dP versus pressure in each nanoslit. Fig. 9a lists the results in different nanoslits and the MMP is plotted versus pore size. A low MMP value is considered to represent a better miscible state. The miscibility seems easier to achieve when the pore
Fig. 7. CO2/C10 system in (a) 5-nm nanoslit and (b) bulk phase at 6.29 MPa from 0.1 to 5 ns. CO2: red, C10: blue; (c) Densities of CO2 and C10 along z-axis at 5 ns.
648
T. Fang et al. / International Journal of Heat and Mass Transfer 141 (2019) 643–650
Fig. 8. Density evolution in 5-nm nanoslit at different pressure. (a) CO2 and (b) oil in nanoslit (c) CO2 and (d) oil near pore surface.
Fig. 9. (a) MMP versus pore size (b) Number density of C10 at 6.29 MPa (c) Free energy barrier DF of C10 in CO2 flooding versus pore size.
size is reduced because of stronger confinement effect. The trend of MMP value is in agreement with the published results by theoretical method [2,22,24]. In addition, it is noteworthy that in Fig. 9a, the MMP value starts to increase once the pore radius is smaller than 5-nm. It is because that the proportion of the adsorbed C10 increases with reducing pore size, while C10 molecules in the free state decrease. The size of C10 is schematically illustrated as an insert in Fig. 9b. The diameter of C10 is around 1.32 nm, so C10 could diffuse freely in a nanoslit with width larger than 5-nm. Obviously, a necessary condition for the fluid miscibility development is that both liquid and gas molecules can diffuse freely in the porous medium [2]. The above analysis proves that if the effect of the molecular adsorption layer is considered, the optimum pore width value will increase, and the stability of the adsorption layer will make the molecular migration more uncertain. In the latter study, the specific effects of the adsorption layer and different crude oil components on the miscible behavior need further analysis.
The free energy barrier DF for C10 in CO2 flooding is shown in Fig. 9c, presenting a similar pattern with the MMP of mixing. To be specific, the DF decreases with the reduction of pore size to 5nm, and it starts to increase afterwards. On the whole, a smaller DF leads to easier miscibility. Thus, the confinement effect facilitates the miscibility of CO2/oil. However, there is a bottom limit in the pore size reduction, which is a necessary condition for the occurrence of miscibility [2]. When the free-molecule area is larger than a single C10 molecule, the miscible state could be achieved. Otherwise, the miscible state would be hard to reach. 4. Conclusions In this study, molecular dynamic simulations were performed to study the miscible mechanism of CO2/oil phase in nanoslits. It has important economic benefits as the generation of all such properties using a rigorously experimental approach is formidable. This is especially true for CO2-EOR process involving the molecular
T. Fang et al. / International Journal of Heat and Mass Transfer 141 (2019) 643–650
microscopic interaction behavior at a nanoscale. In this study, first, the interface thickness of the oil phase by introducing the Gibbs dividing surface was determined. The method helps address the problem that the MMP value is hard to quantify through molecular simulation. In addition, the feasibility of DIM was also confirmed and extended to the calculation of MMP in nanoslits, thus helping to analyze the effect of molecular dynamics on MMP. Second, the results showed that the MMPs in nanoslits were significantly smaller than those in bulk phase, suggesting that oil and gas phases were more likely to attain a miscible state. To explain the phenomenon, the PMFs of C10 in CO2 flooding and competitive adsorption behavior of fluid molecules were analyzed, providing a microscopic explanation for the miscible mechanism in confined spaces. Rising the probability of interaction between C10 and CO2 is a vital factor affecting the miscible phase. Lastly, based on the effect of nanoslits on molecules, five confined nanoslits with the pore size ranging from 2 to 15 nm were examined. The simulation results revealed that the miscible state was easier to achieve with the reduction of pore size due to the confinement effect. Besides, it was found that the width of nanoslit with the lowest MMP was larger than the reported value due to the consideration of adsorption layer, thus making the molecular migration more uncertain. Our work provided the understanding of the miscible mechanism on a molecular level, and the results have significant promise for the oil exploitation and might show some promising applications in the future. Declaration of Competing Interest The authors declared that there is no conflict of interest. Acknowledgements This work is financially supported by the National Natural Science Foundation of China (51874332, U1663206, 51622405), Climb Taishan Scholar Program in Shandong Province (tspd20161004) and Fundamental Research Funds for the Central Universities (YCX2017067, 18CX06006A). Appendix A. Supplementary material Supplementary data to this article can be found online at https://doi.org/10.1016/j.ijheatmasstransfer.2019.06.083. References [1] S. Mondal, S. De, CO2 based power cycle with multi-stage compression and intercooling for low temperature waste heat recovery, Energy 90 (2015) 1132– 1143. [2] K. Zhang, N. Jia, S. Li, L. Liu, Thermodynamic phase behaviour and miscibility of confined fluids in nanopores, Chem. Eng. J. 351 (2018) 1115–1128. [3] W. Yu, H.R. Lashgari, K. Wu, CO2 injection for enhanced oil recovery in Bakken tight oil reservoirs, Fuel 159 (2015) 354–363. [4] M.G. Rezk, J. Foroozesh, Determination of mass transfer parameters and swelling factor of CO2-oil systems at high pressures, Int. J. Heat Mass Transf. 126 (2018) 380–390. [5] T. Lee, L. Bocquet, B. Coasne, Activated desorption at heterogeneous interfaces and long-time kinetics of hydrocarbon recovery from nanoporous media, Nat. Commun. 7 (2016). [6] B. Wei, L. Lu, W. Pu, R. Wu, X. Zhang, Y. Li, F. Jin, Production dynamics of CO2 cyclic injection and CO2 sequestration in tight porous media of Lucaogou formation in Jimsar sag, J. Petrol. Sci. Eng. 157 (2017) 1084–1094. [7] P. Gao, B.F. Towler, G. Pan, Strategies for evaluation of the CO2 miscible flooding process, Abu Dhabi International Petroleum Exhibition and Conference, Society of Petroleum Engineers, 2010. [8] K. Zhang, B. Kong, J. Zhan, R. He, T. Qin, K. Wu, S. Chen, J. Zhang, Effects of nanoscale pore confinement on CO2 immiscible and miscible processes, SPE Low Perm Symposium, Society of Petroleum Engineers, 2016. [9] P. Nguyen, D. Mohaddes, J. Riordon, H. Fadaei, P. Lele, D. Sinton, Fast fluorescence-based microfluidic method for measuring minimum miscibility pressure of CO2 in crude oils, Anal. Chem. 87 (2015) 3160–3164.
649
[10] K. Zhang, N. Jia, F. Zeng, Application of predicted bubble-rising velocities for estimating the minimum miscibility pressures of the light crude oil-CO2 systems with the rising bubble apparatus, Fuel 220 (2018) 412–419. [11] M. Escrochi, N. Mehranbod, S. Ayatollahi, The gas-oil interfacial behavior during gas injection into an asphaltenic oil reservoir, J. Chem. Eng. Data 58 (2013) 2513–2526. [12] O.N. Amézquita, S. Enders, P.T. Jaeger, R. Eggers, Interfacial properties of mixtures containing supercritical gases, J. Supercrit. Fluid 55 (2010) 724–734. [13] A. Mejía, M. Cartes, H. Segura, E.A. Müller, Use of equations of state and coarse grained simulations to complement experiments: describing the interfacial properties of carbon dioxide + decane and carbon dioxide + eicosane mixtures, J. Chem. Eng. Data 59 (2014) 2928–2941. [14] D.N. Rao, J.I. Lee, Application of the new vanishing interfacial tension technique to evaluate miscibility conditions for the Terra Nova Offshore Project, J. Petrol. Sci. Eng. 35 (2002) 247–262. [15] R. Wang, F. Peng, K. Song, G. Feng, Z. Guo, Molecular dynamics study of interfacial properties in CO2 enhanced oil recovery, Fluid Phase Equilibr. 467 (2018) 25–32. [16] T. Babadagli, Development of mature oil fields-a review, J. Petrol. Sci. Eng. 57 (2017) 221–246. [17] C.R. Clarkson, N. Solano, R.M. Bustin, A.M.M. Bustin, G.R.L. Chalmers, L. He, Y.B. Melnichenko, A.P. Radlin´ski, T.P. Blach, Pore structure characterization of North American shale gas reservoirs using USANS/SANS, gas adsorption, and mercury intrusion, Fuel 103 (2013) 606–616. [18] L. Bocquet, E. Charlaix, Nanofluidics, from bulk to interfaces, Chem. Soc. Rev. 39 (2010) 1073–1095. [19] S. Wang, Q. Feng, F. Javadpour, Q. Hu, K. Wu, Competitive adsorption of methane and ethane in montmorillonite nanopores of shale at supercritical conditions: a grand canonical Monte Carlo simulation study, Chem. Eng. J. 355 (2019) 76–90. [20] S.B. Hawthorne, D.J. Miller, L. Jin, C.D. Gorecki, Rapid and simple capillaryrise/vanishing interfacial tension method to determine crude oil minimum miscibility pressure: pure and mixed CO2, methane, and ethane, Energ. Fuel 30 (2016) 6365–6372. [21] A. Fazlali, M. Nikookar, A. Agha-Aminiha, A.H. Mohammadi, Prediction of minimum miscibility pressure in oil reservoirs using a modified SAFT equation of state, Fuel 108 (2013) 675–681. [22] K. Zhang, N. Jia, F. Zeng, P. Luo, A New diminishing interface method for determining the minimum miscibility pressures of light oil-CO2 systems in bulk phase and nanopores, Energ. Fuel 31 (2017) 12021–12034. [23] A. Sharbatian, A. Abedini, Z. Qi, D. Sinton, Full characterization of CO2-oil properties on-chip: solubility, diffusivity, extraction pressure, miscibility, and contact angle, Anal. Chem. 90 (2018) 2461–2467. [24] K. Zhang, S. Seetahal, D. Alexander, J. Lv, Y. Hu, X. Lu, D. Zhang, K. Wu, Z. Chen, Effect of confinement on gas and oil relative permeability during CO2 flooding in tight oil reservoirs, SPE Trinidad and Tobago Section Energy Resources Conference, Society of Petroleum Engineers, 2016. [25] Y. Zhang, H.R. Lashgari, Y. Di, K. Sepehrnoori, Capillary pressure effect on phase behavior of CO2/hydrocarbons in unconventional reservoirs, Fuel 197 (2017) 575–582. [26] H. Wang, X. Wang, X. Jin, D. Cao, Molecular dynamics simulation of diffusion of shale oils in montmorillonite, J. Phys. Chem. C 120 (2016) 8986–8991. [27] Y. Yan, Z. Dong, Y. Zhang, P. Wang, T. Fang, J. Zhang, CO2 activating hydrocarbon transport across nanopore throat: insights from molecular dynamics simulation, Phys. Chem. Chem. Phys. 19 (2017) 30439–30444. [28] T. Le, A. Striolo, D.R. Cole, Propane simulated in silica pores: Adsorption isotherms, molecular structure, and mobility, Chem. Eng. Sci. 121 (2015) 292– 299. [29] T. Le, A. Striolo, D.R. Cole, CO2-C4H10 mixtures simulated in silica slit pores: relation between structure and dynamics, J. Phys. Chem. C 119 (2015) 15274– 15284. [30] T. Fang, M. Wang, J. Li, B. Liu, Y. Shen, Y. Yan, J. Zhang, Study on the asphaltene precipitation in CO2 flooding: a perspective from molecular dynamics simulation, Ind. Eng. Chem. Res. 57 (2018) 1071–1077. [31] D. Hou, Z. Lu, X. Li, H. Ma, Z. Li, Reactive molecular dynamics and experimental study of graphene-cement composites: Structure, dynamics and reinforcement mechanisms, Carbon 115 (2017) 188–208. [32] D. Hou, T. Li, P. Wang, Molecular dynamics study on the structure and dynamics of NaCl solution transport in the nanometer channel of CASH gel, ACS Sustain. Chem. Eng. 6 (2018) 9498–9509. [33] Material Studio, Accelrys Software Inc.: San Diego, CA, 2010. [34] Z. Yang, M. Li, B. Peng, M. Lin, Z. Dong, Dispersion property of CO2 in oil. 1. Volume expansion of CO2 + alkane at near critical and supercritical condition of CO2, J. Chem. Eng. Data 57 (2012) 882–889. [35] Y. Yan, C. Li, Z. Dong, T. Fang, B. Sun, J. Zhang, Enhanced oil recovery mechanism of CO2 water-alternating-gas injection in silica nanochannel, Fuel 190 (2017) 253–259. [36] T. Dong, N.B. Harris, Pore size distribution and morphology in the horn river shale, middle and upper Devonian, northeastern British Columbia, Canada. AAPG Mem. 102 (2013) 67–79. [37] Z.Z. Li, T. Min, Q. Kang, Y.L. He, W.Q. Tao, Investigation of methane adsorption and its effect on gas transport in shale matrix through microscale and mesoscale simulations, Int. J. Heat. Mass. Transf. 98 (2016) 675–686. [38] S. Wang, Q. Feng, M. Zha, F. Javadpour, Q. Hu, Supercritical methane diffusion in shale nanopores: effects of pressure, mineral types, and moisture content, Energ. Fuel 32 (2017) 169–180.
650
T. Fang et al. / International Journal of Heat and Mass Transfer 141 (2019) 643–650
[39] S. Plimpton, Fast parallel algorithms for short-range molecular dynamics, J. Comput. Phys. 117 (1995) 1–19. [40] R.D. Shaver, R.L. Robinson Jr, K.A.M. Gasem, An automated apparatus for equilibrium phase compositions, densities, and interfacial tensions: data for carbon dioxide+ decane, Fluid Phase Equilibr. 179 (2001) 43–66. [41] S. Nosé, A molecular dynamics method for simulations in the canonical ensemble, Mol. Phys. 52 (1984) 255–268. [42] W. Humphrey, A. Dalke, K. Schulten, VMD: visual molecular dynamics, J. Mol. Graph. 14 (1996) 33–38. [43] R.T. Cygan, J.J. Liang, A.G. Kalinichev, Molecular models of hydroxide, oxyhydroxide, and clay phases and the development of a general force field, J. Phys. Chem. B 108 (2004) 1255–1266. [44] S. Iglauer, M.S. Mathew, F. Bresme, Molecular dynamics computations of brine-CO2 interfacial tensions and brine-CO2-quartz contact angles and their effects on structural and residual trapping mechanisms in carbon geosequestration, J. Colloid. Interf. Sci. 386 (2012) 405–414. [45] A.P. Thompson, S.J. Plimpton, W. Mattson, General formulation of pressure and stress tensor for arbitrary many-body interaction potentials under periodic boundary conditions, J. Chem. Phys. 131 (2009) 154107.
[46] Q. Xu, J. Jiang, Computational characterization of ultrathin polymer membranes in liquids, Macromolecules 51 (2018) 7169–7177. [47] L. Bocquet, J.L. Barrat, Hydrodynamic boundary conditions, correlation functions, and Kubo relations for confined fluids, Phys. Rev. E 49 (1994) 3079. [48] I.C. Bourg, C.I. Steefel, Molecular dynamics simulations of water structure and diffusion in silica nanopores, J. Phys. Chem. B 116 (2012) 11556– 11564. [49] T. Fang, M. Wang, Y. Gao, Y. Zhang, Y. Yan, J. Zhang, Enhanced oil recovery with CO2/N2 slug in low permeability reservoir: Molecular dynamics simulation, Chem. Eng. Sci. 197 (2019) 204–211. [50] S. Wang, F. Javadpour, Q. Feng, Molecular dynamics simulations of oil transport through inorganic nanopores in shale, Fuel 171 (2016) 74–86. [51] S. Li, K. Zhang, N. Jia, L. Liu, Evaluation of four CO2 injection schemes for unlocking oils from low-permeability formations under immiscible conditions, Fuel 234 (2018) 814–823.