- Email: [email protected]
Molecular dynamics study on oil migration inside silica nanopore
Chemical Physics Letters 678 (2017) 186–191
Contents lists available at ScienceDirect
Chemical Physics Letters journal homepage: www.elsevier.com/locate/cplett
Research paper
Molecular dynamics study on oil migration inside silica nanopore Jichao Sun a, Heng Zhang a, Mei Hu b, Xiangjia Meng c, Shiling Yuan a,⇑ a
Key Laboratory of Colloid and Interface Chemistry, Shandong University, Jinan 250100, China Shandong Institute for Food and Drug Control, Jinan 250100, China c School of Information Engineering, Shandong Youth University of Political Science, Jinan 250103, China b
a r t i c l e
i n f o
Article history: Received 19 March 2017 In final form 17 April 2017 Available online 18 April 2017
a b s t r a c t Two suggested systems of oil in the porosity of the reservoir rock after water flooding were built with an oil cylinder inside a silica nanopore. A series of MD simulations were performed to obtain the aggregation structure of oil phase. The results revealed that heavy oil components showed different distribution position inside silica nanopore in the two oil systems. Heavy oil molecules in heptane can precipitate and adsorbed on silica surface. Steered MD simulation was used to study the oil displacement. By analyzing solvent accessible surface area (SASA), we demonstrated the migration of heavy oil components at molecular level. Ó 2017 Elsevier B.V. All rights reserved.
1. Introduction Recently, oil still plays an important role in the development of society. Many researchers pay attention to how to improve oil recovery. However, with water flooding, reservoir environment gets more and more complex. Injection of excess water leads to the increase of costs in oil/water separation [1,2]. Moreover, residual oil, which is trapped into the pores of the reservoir rock under capillary interaction, is difficult to exploit. To extract residual oil out of the pores, knowledge of the aggregation properties of oil molecules and displacement of oil droplets in the pores is urgently required. A lot of work has been reported in exploring the interfacial properties and interaction between oil/substrate [3,4], oil/oil displacement agent [5], or oil displacement/substrate [6,7]. It is expected to useful in selection of efficient oil displacement agent in enhanced oil recovery. Although much work has been done using experimental and theoretical methods in enhanced oil recovery, the significant problem of how oil molecules transport in the pore of reservoir rock remains unanswered. Petroleum is mixed with organic and inorganic materials. Heavy components, generally contain asphaltenes [8] and resins [9], is readily form aggregate structures in bulk crude oil phase due to the stacking of polyaromatic rings [10]. The aggregation leads to many serious problems during oil recovery and transport; for instance, high viscosity [11], deposition [12], reservoir plugging and production pipelines fouled [13]. Therefore, it is significant to understand the displacement of oil molecules clearly. Mozaffari et al. reported that lower value of volume fraction of bitumen and ⇑ Corresponding author. E-mail address: [email protected] (S. Yuan). http://dx.doi.org/10.1016/j.cplett.2017.04.057 0009-2614/Ó 2017 Elsevier B.V. All rights reserved.
asphaltenes was observed in heptane [14]. It correlates well with the insolubility of asphaltenes and higher precipitation rate in heptane. However, asphaltenes and resins can highly solubilize in toluene [15] due to their polyaromatic structures. Although the behavior of asphaltenes and resins aggregated at oil/water interface are well investigated [16–18], little attention has been given to exploring the displacement of heavy oil molecules inside a nanopore [19], which can be suggested as a model for oil recovery. It is well-known that, because of polyaromatic associations in crude oil, the condensed polyaromatic fraction of crude oil is easily precipitated and deposited during production and transportation, which leads to blocking of reservoir rocks and transport pipes [12,20,21]. Moreover, the polyaromatic association phenomena in oil phase strongly effect solubility, viscosity, density, and other physical properties of crude oil [20,21]. As a result, that will provide problems in oil production and displacement. In this work, we performed molecular dynamics (MD) simulation and steered molecular dynamics (SDM) [22] simulation to study oil aggregation and displacement inside a cylindrical silica nanopore using the GROMACS package [23]. Silica was selected as capillary pore because of it is a main component of rock minerals in many geological environments. And, it is used as micro model to investigate residual oil detachment in many experimental studies [5,19,24]. Two simulated systems were built with different main components. Toluene was selected as main component for system A. Heptane was selected for system B. Asphaltenes and resins were added as heavy oil components both in the two systems. We focused on the adsorption of different oil molecules on the pore surface with equilibrium MD simulation. Steered MD simulation was employed to study the displacement process of
J. Sun et al. / Chemical Physics Letters 678 (2017) 186–191
an oil drop which was applied external force. The observation was analyzed in detail and may be useful in understanding oil detachment mechanism for enhanced oil recovery.
2. Simulation details All the MD simulations were performed with GROMACS 4.5.4 software package [23]. The united-atom optimized performance for the liquid systems (GROMOS 53a6) force field was adopted for all of the potential function terms to calculate the interatomic interactions. The simple point charge (SPC) model was used to describe water molecules [25]. The force field parameters of silica nanopore were referred to the previous publications, which had obtained good results in the investigation of dynamics properties of amorphous silica [26,27]. The size and force field parameters of silica nanopore was derived from previous work [24]. Firstly, the energies of each of the initial configurations were minimized by steepest descent method for each of the systems. Then, a 10 ns MD simulation under canonical ensemble (NVT) was carried out for every system with periodic boundary conditions applied in all directions. The position of silica block was restrained without surface atoms during the minimization and equilibration MD simulation. In all simulations, the V-rescale thermostat algorithm [28] was used to keep the temperature constant at 298 K. LINCS algorithm [29] was adopted to constrain bond lengths. The Lennard-Jones interactions were cut off at 1.2 nm for the nonbonded potential. Electrostatics interactions calculated using the particle mesh Ewald (PME) method [30].
3. Results and discussion 3.1. MD simulation The cylindrical oil phase was put into silica pore with water filled up. After 10 ns equilibrium MD simulation, the oil phases of two systems showed very different configuration. Fig. 1 shows the two-dimensional distribution maps of heavy oil components which were calculated on x-y plane (normal to the direction of silica pore). It can be seen that there is no obvious difference compared to initial configuration after MD simulation in system A. Interestingly, a complete different two-dimensional distribution map was obtained in system B with 10 ns MD simulation. The circular distribution map means that heavy oil molecules moved toward silica surface and adsorbed on it under strong nonbonded interactions. Another reason is the insolubility and higher precipitation rate of heavy components in heptane. To have a better understanding of changes of oil phase, radial relative density profile distribution (Fig. 2) was calculated every 0.1 nm. The area of central circle, which radius is 0.1 nm, was set as 1. On this basis, the area of other annuluses was computed respectively. Then, the density profile was obtained with the achievement of atom number of heavy oil. Due to the small area, error of the density profile within 0.3 nm is too large to be considered. We can see that the peaks of radial relative density profile distribution curve become higher with its number decreased after MD simulation in system A (Fig. 2a). However, the main distribution region was not changed. That indicates heavy oil molecules aggregated but not adsorbed on silica surface. In contrast, the heavy oil distribution region of system B obviously moved away from silica nanopore center (Fig. 2b). That indicates heavy oil molecules can supplant heptane molecules and adsorb on silica surface. The reduction of main distribution region means heavy oil molecules also aggregated in system B.
187
3.2. Steered MD simulation To simulate the driving force from the displacing fluid during the oil recovery process, an external force was applied on the center of mass (COM) of oil phase along the silica nanopore after 10 ns equilibrium MD simulation. As a result, Oil molecules moved gradually along the silica pore under applied force. After 3 ns SMD simulation, the oil phase was completely separated into two sections, residual oil and separate oil. Fig. 3 shows the snapshots of simulation systems at the moment of oil phase almost separated into two sections. With the separation of oil phase, separate oil located at the pore center and can move along the silica nanopore under applied force. Because of large space between separate oil and silica surface, separate oil moved readily under weak adsorption interaction of silica surface. Residual oil still adsorbed on silica surface and kept initial position due to strong interaction between them. Interestingly, residual oil showed different conformation in the pore (Fig. 3). In system A, residual oil molecules almost had not moved under external force and the impact of water and separate oil molecules. They remained initial position and adsorbed on silica surface after 3 ns SMD simulation. However, in spite of adsorption interaction caused by silica nanopore, residual oil moved along the pore in system B and showed a large distribution area in the pore. This indicates that residual oil can be displaced gradually with the impact of external force, water and separate molecules in system B. To illustrate the movement of oil molecules inside the pore, the distance between the COM of oil phase and the nonmoving reference with time evolution (Fig. 4) was calculated for the two oil sections. We can see from Fig. 4a that the slops of the curves increased with time evolution before 1.7 ns. It should be noticed that the increase of the slop became faster after 1 ns. Combining Figs. 3 and 4, the pulling process can be separated into 3 stages. The first stage (0–1.0 ns) named detachment process. Because oil phase was unseparated, strong interactions between silica surface and oil molecules lead to difficult displacement of oil phase. Compared to oil molecules near silica surface, those oil molecules which located at pore center moved slowly under nonbonded interactions of other adsorbed components. The second stage (1.0–1.7 ns) was called quick displacement process. Oil phase was separated and separate oil broke away from initial position gradually. The impact of residual oil reduced and eliminated gradually with time evolution. Thus, separate oil was displaced quickly under external force and water scouring. The third stage (1.7–3.0 ns) named adsorption blockage process. Due to the attractive interaction of silica nanopore, separate oil got close to silica surface while moving along the pore. Finally, oil molecules displaced water layer and adsorbed on silica surface. As a result, the movement of separate oil became slow down. From Fig. 4a, we also found that the distance-time curves of the two systems showed similar variation before 1.7 ns. However, during the third stage of system B, an inflection point and a peak appeared at 2.0 ns and 2.75 ns respectively. That illustrates the adsorption of oil molecules in system B was not very strong, but with obvious adsorption-desorption process. As described above, we believe that there are differences between the two simulation systems, especially for residual oil. But, it is not obviously pronounced in Fig. 4a because of the large spacing and weak interactions among separate oil and silica surface. Herein, the distancetime curves of residual oil are depicted in Fig. 4b. Different results under same external force reveal different displacement ability. Compared to system B, a relatively smooth curve and very short distance of system A means it is difficult to displace residual oil from silica surface with external force. In system B, the distance increased faster and larger during the whole simulation process indicates that residual oil is more easily pulled away from
188
J. Sun et al. / Chemical Physics Letters 678 (2017) 186–191
72
72
(a)
28.00 24.00 20.00 16.00 12.00 8.000 4.000 0
y
48 36 24
60
(b)
28.00 24.00 20.00 16.00 12.00 8.000 4.000 0
48
y
60
36 24
12
12 24
36
48
60
72
84
24
36
48
x
72
84
x
72 60
60
72
(c)
60
24.00
(d)
24.00 20.00
20.00
48
48
16.00
16.00 12.00
y
36
8.000
y
12.00
36
8.000 4.000
4.000
24
24
0
12
0
12 24
36
48
60
72
84
24
36
48
x
60
72
84
x
Fig. 1. Two-dimensional distributions of heavy oil molecules on x-y plane. (a) System A at 0 ns; (b) system A at 10 ns; (c) system B at 0 ns; (d) system B at 10 ns.
5
(a)
(b)
0ns 10ns
Relative density profile
Relative density profile
5 4 3 2 1
0ns 10ns
4 3 2 1 0
0 0.0
0.5
1.0
1.5
2.0
2.5
3.0
r (nm)
0.0
0.5
1.0
1.5
2.0
2.5
3.0
r (nm)
Fig. 2. Radial relative density profile distribution of heavy oil molecules. (a) System A; (b) system B.
silica surface. In addition, residual oil migrated with adsorptiondesorption process which can be evident by the two peaks appeared in the third stage. According to the discussion above, we think that oil in system A migrated difficultly inside the silica nanopore and readily resulted in pore blockage. In MD simulation, different distribution of heavy oil inside the silica nanopore resulted in changes of light components distribution. Actually, the adsorption and migration of heavy oil will impact on light components, especially during SMD simulation. The two systems have same heavy oil molecules. Thus, to investigate the displacement of oil phase in silica nanopore, it is important to focus on light oil components. Fig. 5 shows the RDFs of light oil molecules near the silica nanopore surface. The RDFs indicate how the light oil molecules associated with the pore surface. We can see that there is distinctively different light oil molecules
distribution between the two systems. Weak interaction between oil molecules and silica nanopore surface leads to an increase in distance among the two components, which is caused by the competitive adsorption of heavy oil constituents. The toluene distribution near silica nanopore surface in Fig. 5 shows a much higher g(r) and its prominent peak appears closer to silica surface as compared to heptane. And, the distribution curve of toluene appears earlier than heptane. To have a better understanding of the transportation process, the nonbonded interaction energy between silica surface and light oil molecules with time evolution was calculated and shown in Fig. S1. In MD simulation, light oil components in system A showed stronger nonbonded interactions with silica surface. Moreover, most of light oil molecules in system B were separated from silica surface by heavy oil components which were adsorbed on silica
189
J. Sun et al. / Chemical Physics Letters 678 (2017) 186–191
Fig. 3. The snapshots of oil distribution at different simulation time. (a) System A; (b) system B. Asphaltenes, red; resins, green; light components, brown. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
16
Distance (nm)
4
(a) system A system B
(b) system A system B
3
Distance (nm)
20
12 8 4
2
1
0
0 0.0
0.5
1.0
1.5
2.0
2.5
3.0
t (ns)
0.0
0.5
1.0
1.5
2.0
2.5
3.0
t (ns)
Fig. 4. The distance between the COMs of oil phase and the nonmoving reference with time evolution. (a) Separate oil; (b) residual oil.
more notable (Fig. S3). Because heavy oil molecules have larger surface area than light oil molecules, the effects of water molecules was more pronounced. Herein, residual oil in system B can move gradually under external force and water scouring.
3.0 system A system B
2.4
3.3. Solvent accessible surface area
g (r)
1.8 1.2 0.6 0.0 0.0
0.4
0.8
1.2
1.6
2.0
r (nm) Fig. 5. RDFs of light oil molecules near silica nanopore surface averaged over the last 1 ns.
surface. They can move firstly under applied force. As a result, residual oil which adsorbed on silica surface was exposed to bulk water. In spite of much stronger nonbonded interaction between silica surface and heavy oil molecules in system B (Fig. S2), with time evolution during SMD simulation, effects of water molecules on heavy oil components in residual oil increased and became
Fig. 6a shows change of solvent accessible surface area (SASA), which is only effected by aggregation of heavy oil molecules (without consideration of influence of other components) while calculation, with time evolution. At the beginning of MD simulation in system A, heavy oil molecules aggregated together resulted in quick reduction of SASA. With most of the molecules had been assembled, the amplitude of SASA changes declined after 2.5 ns. Without the obvious adsorption of heavy oil on silica nanopore wall, it is negligible that the changes of SASA, which is caused by the influence of silica surface on the conformation of aggregates. In system B, same SASA-time curve can be found before 6 ns in MD simulation. With time evolution, heavy oil molecules move closer to silica surface and adsorbed on it under nonbonded interactions. As a result, part of the aggregates was destroyed. SASA increased in the last 4 ns because of the dispersion of heavy oil molecules. However, the heavy components were not dispersed completely. SASA at 10 ns is less than its initial value. In SMD simulation, the aggregates dispersed as individual molecules quickly under external force. Oil molecules distributed larger range along silica pore (Fig. S4). That led to rapid increase
2
2
Solvent accessible surface area (nm )
J. Sun et al. / Chemical Physics Letters 678 (2017) 186–191
Solvent accessible surface area (nm )
190
MD
(a)
190
SMD
system A system B
180 170 160 150 140 130
0
2
4
6
8
0
2
(b)
35
MD
SMD
system A system B
30 25 20 15 10 0
2
4
t (ns)
6
8
0
2
t (ns)
Fig. 6. Solvent accessible surface area of heavy oil molecules over the simulation time of 10 ns in MD simulation and 3 ns in SMD simulation. (a) Without affection of silica pore surface; (b) Reduction of solvent accessible surface area caused by adsorbing to silica pore surface.
of SASA in both of the simulation systems. The notable difference is that an extremal point appeared at 1.4 ns in system B. This is because of relatively strong interaction among heavy oil molecules resulted in reaggregation occurred in separate oil phase. Therefore, although separate oil had been pulled away from residual oil, SASA decreased after 1.4 ns. Fig. 6b reveals reduction of SASA caused by adsorption of heavy oil molecules onto silica pore surface. In system A, we found that SASA remained stable after slow increase in 2.5 ns MD simulation. That means only small part of heavy oil molecules, which located close to silica surface in initial configuration, was adsorbed on pore wall. Quite differently, SASA increased within initial 7.5 ns in system B. Then, it followed a slight decrease in last 2.5 ns. That denotes most of heavy oil molecules were adsorbed on silica surface. Finally, with the saturated adsorption, conformation change or desorption of individual heavy oil molecule led to the decrease of SASA. During the SMD simulation process, we can see that SASA increased invariably in system A. Separate oil was pulled away from initial position under external force. Its heavy oil components were adsorbed on silica surface gradually because of losing protection from residual oil. As a result, SASA that occupied by silica pore wall increased. In system B, a sudden decrease of SASA appeared at 1 ns and maintained about 0.4 ns, which indicates the desorptionadsorption process of heavy oil molecules in residual oil. Residual oil molecules desorbed from silica surface under external force and water scouring. With short distance displaced, they were adsorbed on silica surface again. Although heavy oil components were adsorbed on silica surface readily in system B, residual oil can be easier detached by external force compared to system A. Therefore, to enhance oil-displacement efficiency, it is significant to reasonably exploit the transportation of oil molecules inside silica nanopore, which will help to design efficient oil displacement agent. 4. Conclusion In this work, a series of MD simulations were performed to research the mechanism of oil migration inside a silica nanopore. With equilibrium MD simulations, the distributions of heavy oil components of crude oil inside silica nanopore were revealed. The results showed that heavy oil molecules dispersed well in toluene. However, due to insolubility in heptane, heavy oil components deposited and adsorbed on silica surface. The steered MD simulations were conducted to investigate oil transport process inside silica pore. Oil molecules in the central part of oil phase migrated firstly under external force. With pulling away from initial position and exposure to silica surface, they were adsorbed on silica pore wall due to nonbonded interactions. Residual oil in sys-
tem B is more easily pulled away from silica surface under external force than system A. The results provide insight vision of oil migration and were expected to promote the enhanced oil recovery. Acknowledgements We gratefully appreciate the financial support from the National Natural Science Foundation of China, China (21573130 and 11647013) and the Shandong Provincial Natural Science Foundation, China (ZR2016BP08).
Appendix A. Supplementary material Supplementary data associated with this article can be found, in the online version, at http://dx.doi.org/10.1016/j.cplett.2017.04. 057. References [1] M.O. Elsharafi, B. Bai, Effect of weak preformed particle gel on unswept oil zones/areas during conformance control treatments, Ind. Eng. Chem. Res. 51 (2012) 11547–11554. [2] R. Seright, Washout of Cr (III)-Acetate-HPAM Gels from Fractures, Presented at the SPE international symposium on oilfield chemistry, Houston, Texas, 5–7 February 2003; Paper 80200. [3] Q. Xue, Y. Tao, Z. Liu, S. Lu, X. Li, T. Wu, Y. Jin, X. Liu, Mechanism of oil molecules transportation in nano-sized shale channel: MD simulation, RSC Adv. 5 (2015) 25684–25692. [4] J. Zhong, X. Wang, J. Du, L. Wang, Y. Yan, J. Zhang, Combined molecular dynamics and quantum mechanics study of oil droplet adsorption on different self-assembly monolayers in aqueous solution, J. Phys. Chem. C 117 (2013) 12510–12519. [5] Q. Liu, S. Yuan, H. Yan, X. Zhao, Mechanism of oil detachment from a silica surface in aqueous surfactant solutions: molecular dynamics simulations, J. Phys. Chem. B 116 (2012) 2867–2875. [6] X. Hu, Y. Li, H. Sun, X. Song, Q. Li, X. Cao, Z. Li, Effect of divalent cationic ions on the adsorption behavior of zwitterionic surfactant at silica/solution interface, J. Phys. Chem. B 114 (2010) 8910–8916. [7] N.R. Tummala, L. Shi, A. Striolo, Molecular dynamics simulations of surfactants at the silica-water interface: anionic vs nonionic headgroups, J. Colloid Interface Sci. 362 (2011) 135–143. [8] T. Takanohashi, S. Sato, R. Tanaka, Structural relaxation behaviors of three different asphaltenes using MD calculations, Pet. Sci. Technol. 22 (2004) 901– 914. [9] O. Castellano, R. Gimon, C. Canelon, Y. Aray, H. Soscun, Molecular interactions between Orinoco belt resins, Energy Fuels 26 (2012) 2711–2720. [10] O. Mullins, S. Betancourt, M. Cribbs, F. Dubost, J. Creek, B. Andrews, L. Venkataramanan, The colloidal structure of crude oil and the structure of oil reservoirs, Energy Fuels 21 (2007) 2785–2794. [11] C. Pierre, L. Barré, A. Pina, M. Moan, Composition and heavy oil rheology, Oil Gas Sci. Technol. 59 (2006) 489–501. [12] P.M. Spiecker, K.L. Gawrys, C.B. Trail, P.K. Kilpatrick, Effects of petroleum resins on asphaltene aggregation and water-in-oil emulsion formation, Colloids and Surf. A 220 (2003) 9–27.
J. Sun et al. / Chemical Physics Letters 678 (2017) 186–191 [13] V.A.M. Branco, G.A. Mansoori, L.C.D.A. Xavier, S.J. Park, H. Manafi, Asphaltene flocculation and collapse from petroleum fluids, J. Pet. Sci. Eng. 32 (2001) 217– 230. [14] S. Mozaffari, P. Tchoukov, J. Atias, J. Czarnecki, N. Nazemifard, Effect of asphaltene aggregation on rheological properties of diluted athabasca bitumen, Energy Fuels 29 (2015) 5595–5599. [15] L. Goual, A. Firoozabadi, Measuring asphaltenes and resins, and dipole moment in petroleum fluids, AIChE J. 48 (2002) 2646–2663. [16] F. Gao, Z. Xu, G. Liu, S. Yuan, Molecular dynamics simulation: the behavior of asphaltene in crude oil and at the oil/water interface, Energy Fuels 28 (2014) 7368–7376. [17] G. Lv, F. Gao, G. Liu, S. Yuan, The properties of asphaltene at the oil-water interface: a molecular dynamics simulation, Colloids and Surf. A 515 (2017) 34–40. [18] J. Zhang, D. Tian, M. Lin, Z. Yang, Z. Dong, Effect of resins, waxes and asphaltenes on water-oil interfacial properties and emulsion stability, Colloids and Surf. A 507 (2016) 1–6. [19] H. Yan, S. Yuan, Molecular dynamics simulation of the oil detachment process within silica nanopores, J. Phys. Chem. C 120 (2016) 2667–2674. [20] J.G. Speight, Petroleum asphaltenes – Part 1 – Asphaltenes, resins and the structure of petroleum, Oil Gas Sci. Technol. 59 (2004) 467–477. [21] C. Pierre, L. Barre, A. Pina, M. Moan, Composition and heavy oil rheology, Oil Gas Sci. Technol. 59 (2004) 489–501.
191
[22] S. Izrailec, S. Stepaniants, B. Isralewitz, D. Kosztin, H. Lu, F. Molnar, W. Wriggers, K. Schulten, Steered Molecular Dynamics, Springer-Verlag, Berlin, 1998. [23] D. van der Spoel, E. Lindahl, B. Hess, The GROMCAS development team, GROMACS User Manual version 4.6.3, Royal Institute of Technology and Uppsala University, Stockholm and Uppsala, 2013. [24] H. Zhang, Y. Ma, Q. Hao, H. Wang, G. Liu, S. Yuan, Molecular dynamics study on mechanism of preformed particle gel transporting through nanopores: surface hydration, RSC Adv. 6 (2016) 7172–7180. [25] H.J.C. Berendsen, J.P.M. Postma, W.F. Van Gunsteren, J. Hermans. In: B. Pullman (Ed.), Intermolecular Forces, D. Reidel Publishing Company, Dordrecht, 1981. [26] M.D. Elola, J. Rodriguez, D. Laria, Liquid methanol confined within functionalized silica nanopores. 2. Solvation dynamics of coumarin 153, J. Phys. Chem. B 115 (2011) 12859–12867. [27] M.D. Elola, J. Rodriguez, D. Laria, Structure and dynamics of liquid methanol confined within functionalized silica nanopores, J. Chem. Phys. 133 (2010) 154707. [28] G. Bussi, D. Donadio, M. Parrinello, Canonical sampling through velocity rescaling, J. Chem. Phys. 126 (2007) 014101. [29] B. Hess, H. Bekker, H.J.C. Berendsen, J. Fraaije, Lincs: a linear constraint solver for molecular simulations, J. Comput. Chem. 18 (1997) 1463–1472. [30] U. Essmann, L. Perera, M.L. Berkowitz, T. Darden, H. Lee, L. Pedersen, A smooth particle mesh ewald method, J. Chem. Phys. 103 (1995) 8577–8593.
Contents lists available at ScienceDirect
Chemical Physics Letters journal homepage: www.elsevier.com/locate/cplett
Research paper
Molecular dynamics study on oil migration inside silica nanopore Jichao Sun a, Heng Zhang a, Mei Hu b, Xiangjia Meng c, Shiling Yuan a,⇑ a
Key Laboratory of Colloid and Interface Chemistry, Shandong University, Jinan 250100, China Shandong Institute for Food and Drug Control, Jinan 250100, China c School of Information Engineering, Shandong Youth University of Political Science, Jinan 250103, China b
a r t i c l e
i n f o
Article history: Received 19 March 2017 In final form 17 April 2017 Available online 18 April 2017
a b s t r a c t Two suggested systems of oil in the porosity of the reservoir rock after water flooding were built with an oil cylinder inside a silica nanopore. A series of MD simulations were performed to obtain the aggregation structure of oil phase. The results revealed that heavy oil components showed different distribution position inside silica nanopore in the two oil systems. Heavy oil molecules in heptane can precipitate and adsorbed on silica surface. Steered MD simulation was used to study the oil displacement. By analyzing solvent accessible surface area (SASA), we demonstrated the migration of heavy oil components at molecular level. Ó 2017 Elsevier B.V. All rights reserved.
1. Introduction Recently, oil still plays an important role in the development of society. Many researchers pay attention to how to improve oil recovery. However, with water flooding, reservoir environment gets more and more complex. Injection of excess water leads to the increase of costs in oil/water separation [1,2]. Moreover, residual oil, which is trapped into the pores of the reservoir rock under capillary interaction, is difficult to exploit. To extract residual oil out of the pores, knowledge of the aggregation properties of oil molecules and displacement of oil droplets in the pores is urgently required. A lot of work has been reported in exploring the interfacial properties and interaction between oil/substrate [3,4], oil/oil displacement agent [5], or oil displacement/substrate [6,7]. It is expected to useful in selection of efficient oil displacement agent in enhanced oil recovery. Although much work has been done using experimental and theoretical methods in enhanced oil recovery, the significant problem of how oil molecules transport in the pore of reservoir rock remains unanswered. Petroleum is mixed with organic and inorganic materials. Heavy components, generally contain asphaltenes [8] and resins [9], is readily form aggregate structures in bulk crude oil phase due to the stacking of polyaromatic rings [10]. The aggregation leads to many serious problems during oil recovery and transport; for instance, high viscosity [11], deposition [12], reservoir plugging and production pipelines fouled [13]. Therefore, it is significant to understand the displacement of oil molecules clearly. Mozaffari et al. reported that lower value of volume fraction of bitumen and ⇑ Corresponding author. E-mail address: [email protected] (S. Yuan). http://dx.doi.org/10.1016/j.cplett.2017.04.057 0009-2614/Ó 2017 Elsevier B.V. All rights reserved.
asphaltenes was observed in heptane [14]. It correlates well with the insolubility of asphaltenes and higher precipitation rate in heptane. However, asphaltenes and resins can highly solubilize in toluene [15] due to their polyaromatic structures. Although the behavior of asphaltenes and resins aggregated at oil/water interface are well investigated [16–18], little attention has been given to exploring the displacement of heavy oil molecules inside a nanopore [19], which can be suggested as a model for oil recovery. It is well-known that, because of polyaromatic associations in crude oil, the condensed polyaromatic fraction of crude oil is easily precipitated and deposited during production and transportation, which leads to blocking of reservoir rocks and transport pipes [12,20,21]. Moreover, the polyaromatic association phenomena in oil phase strongly effect solubility, viscosity, density, and other physical properties of crude oil [20,21]. As a result, that will provide problems in oil production and displacement. In this work, we performed molecular dynamics (MD) simulation and steered molecular dynamics (SDM) [22] simulation to study oil aggregation and displacement inside a cylindrical silica nanopore using the GROMACS package [23]. Silica was selected as capillary pore because of it is a main component of rock minerals in many geological environments. And, it is used as micro model to investigate residual oil detachment in many experimental studies [5,19,24]. Two simulated systems were built with different main components. Toluene was selected as main component for system A. Heptane was selected for system B. Asphaltenes and resins were added as heavy oil components both in the two systems. We focused on the adsorption of different oil molecules on the pore surface with equilibrium MD simulation. Steered MD simulation was employed to study the displacement process of
J. Sun et al. / Chemical Physics Letters 678 (2017) 186–191
an oil drop which was applied external force. The observation was analyzed in detail and may be useful in understanding oil detachment mechanism for enhanced oil recovery.
2. Simulation details All the MD simulations were performed with GROMACS 4.5.4 software package [23]. The united-atom optimized performance for the liquid systems (GROMOS 53a6) force field was adopted for all of the potential function terms to calculate the interatomic interactions. The simple point charge (SPC) model was used to describe water molecules [25]. The force field parameters of silica nanopore were referred to the previous publications, which had obtained good results in the investigation of dynamics properties of amorphous silica [26,27]. The size and force field parameters of silica nanopore was derived from previous work [24]. Firstly, the energies of each of the initial configurations were minimized by steepest descent method for each of the systems. Then, a 10 ns MD simulation under canonical ensemble (NVT) was carried out for every system with periodic boundary conditions applied in all directions. The position of silica block was restrained without surface atoms during the minimization and equilibration MD simulation. In all simulations, the V-rescale thermostat algorithm [28] was used to keep the temperature constant at 298 K. LINCS algorithm [29] was adopted to constrain bond lengths. The Lennard-Jones interactions were cut off at 1.2 nm for the nonbonded potential. Electrostatics interactions calculated using the particle mesh Ewald (PME) method [30].
3. Results and discussion 3.1. MD simulation The cylindrical oil phase was put into silica pore with water filled up. After 10 ns equilibrium MD simulation, the oil phases of two systems showed very different configuration. Fig. 1 shows the two-dimensional distribution maps of heavy oil components which were calculated on x-y plane (normal to the direction of silica pore). It can be seen that there is no obvious difference compared to initial configuration after MD simulation in system A. Interestingly, a complete different two-dimensional distribution map was obtained in system B with 10 ns MD simulation. The circular distribution map means that heavy oil molecules moved toward silica surface and adsorbed on it under strong nonbonded interactions. Another reason is the insolubility and higher precipitation rate of heavy components in heptane. To have a better understanding of changes of oil phase, radial relative density profile distribution (Fig. 2) was calculated every 0.1 nm. The area of central circle, which radius is 0.1 nm, was set as 1. On this basis, the area of other annuluses was computed respectively. Then, the density profile was obtained with the achievement of atom number of heavy oil. Due to the small area, error of the density profile within 0.3 nm is too large to be considered. We can see that the peaks of radial relative density profile distribution curve become higher with its number decreased after MD simulation in system A (Fig. 2a). However, the main distribution region was not changed. That indicates heavy oil molecules aggregated but not adsorbed on silica surface. In contrast, the heavy oil distribution region of system B obviously moved away from silica nanopore center (Fig. 2b). That indicates heavy oil molecules can supplant heptane molecules and adsorb on silica surface. The reduction of main distribution region means heavy oil molecules also aggregated in system B.
187
3.2. Steered MD simulation To simulate the driving force from the displacing fluid during the oil recovery process, an external force was applied on the center of mass (COM) of oil phase along the silica nanopore after 10 ns equilibrium MD simulation. As a result, Oil molecules moved gradually along the silica pore under applied force. After 3 ns SMD simulation, the oil phase was completely separated into two sections, residual oil and separate oil. Fig. 3 shows the snapshots of simulation systems at the moment of oil phase almost separated into two sections. With the separation of oil phase, separate oil located at the pore center and can move along the silica nanopore under applied force. Because of large space between separate oil and silica surface, separate oil moved readily under weak adsorption interaction of silica surface. Residual oil still adsorbed on silica surface and kept initial position due to strong interaction between them. Interestingly, residual oil showed different conformation in the pore (Fig. 3). In system A, residual oil molecules almost had not moved under external force and the impact of water and separate oil molecules. They remained initial position and adsorbed on silica surface after 3 ns SMD simulation. However, in spite of adsorption interaction caused by silica nanopore, residual oil moved along the pore in system B and showed a large distribution area in the pore. This indicates that residual oil can be displaced gradually with the impact of external force, water and separate molecules in system B. To illustrate the movement of oil molecules inside the pore, the distance between the COM of oil phase and the nonmoving reference with time evolution (Fig. 4) was calculated for the two oil sections. We can see from Fig. 4a that the slops of the curves increased with time evolution before 1.7 ns. It should be noticed that the increase of the slop became faster after 1 ns. Combining Figs. 3 and 4, the pulling process can be separated into 3 stages. The first stage (0–1.0 ns) named detachment process. Because oil phase was unseparated, strong interactions between silica surface and oil molecules lead to difficult displacement of oil phase. Compared to oil molecules near silica surface, those oil molecules which located at pore center moved slowly under nonbonded interactions of other adsorbed components. The second stage (1.0–1.7 ns) was called quick displacement process. Oil phase was separated and separate oil broke away from initial position gradually. The impact of residual oil reduced and eliminated gradually with time evolution. Thus, separate oil was displaced quickly under external force and water scouring. The third stage (1.7–3.0 ns) named adsorption blockage process. Due to the attractive interaction of silica nanopore, separate oil got close to silica surface while moving along the pore. Finally, oil molecules displaced water layer and adsorbed on silica surface. As a result, the movement of separate oil became slow down. From Fig. 4a, we also found that the distance-time curves of the two systems showed similar variation before 1.7 ns. However, during the third stage of system B, an inflection point and a peak appeared at 2.0 ns and 2.75 ns respectively. That illustrates the adsorption of oil molecules in system B was not very strong, but with obvious adsorption-desorption process. As described above, we believe that there are differences between the two simulation systems, especially for residual oil. But, it is not obviously pronounced in Fig. 4a because of the large spacing and weak interactions among separate oil and silica surface. Herein, the distancetime curves of residual oil are depicted in Fig. 4b. Different results under same external force reveal different displacement ability. Compared to system B, a relatively smooth curve and very short distance of system A means it is difficult to displace residual oil from silica surface with external force. In system B, the distance increased faster and larger during the whole simulation process indicates that residual oil is more easily pulled away from
188
J. Sun et al. / Chemical Physics Letters 678 (2017) 186–191
72
72
(a)
28.00 24.00 20.00 16.00 12.00 8.000 4.000 0
y
48 36 24
60
(b)
28.00 24.00 20.00 16.00 12.00 8.000 4.000 0
48
y
60
36 24
12
12 24
36
48
60
72
84
24
36
48
x
72
84
x
72 60
60
72
(c)
60
24.00
(d)
24.00 20.00
20.00
48
48
16.00
16.00 12.00
y
36
8.000
y
12.00
36
8.000 4.000
4.000
24
24
0
12
0
12 24
36
48
60
72
84
24
36
48
x
60
72
84
x
Fig. 1. Two-dimensional distributions of heavy oil molecules on x-y plane. (a) System A at 0 ns; (b) system A at 10 ns; (c) system B at 0 ns; (d) system B at 10 ns.
5
(a)
(b)
0ns 10ns
Relative density profile
Relative density profile
5 4 3 2 1
0ns 10ns
4 3 2 1 0
0 0.0
0.5
1.0
1.5
2.0
2.5
3.0
r (nm)
0.0
0.5
1.0
1.5
2.0
2.5
3.0
r (nm)
Fig. 2. Radial relative density profile distribution of heavy oil molecules. (a) System A; (b) system B.
silica surface. In addition, residual oil migrated with adsorptiondesorption process which can be evident by the two peaks appeared in the third stage. According to the discussion above, we think that oil in system A migrated difficultly inside the silica nanopore and readily resulted in pore blockage. In MD simulation, different distribution of heavy oil inside the silica nanopore resulted in changes of light components distribution. Actually, the adsorption and migration of heavy oil will impact on light components, especially during SMD simulation. The two systems have same heavy oil molecules. Thus, to investigate the displacement of oil phase in silica nanopore, it is important to focus on light oil components. Fig. 5 shows the RDFs of light oil molecules near the silica nanopore surface. The RDFs indicate how the light oil molecules associated with the pore surface. We can see that there is distinctively different light oil molecules
distribution between the two systems. Weak interaction between oil molecules and silica nanopore surface leads to an increase in distance among the two components, which is caused by the competitive adsorption of heavy oil constituents. The toluene distribution near silica nanopore surface in Fig. 5 shows a much higher g(r) and its prominent peak appears closer to silica surface as compared to heptane. And, the distribution curve of toluene appears earlier than heptane. To have a better understanding of the transportation process, the nonbonded interaction energy between silica surface and light oil molecules with time evolution was calculated and shown in Fig. S1. In MD simulation, light oil components in system A showed stronger nonbonded interactions with silica surface. Moreover, most of light oil molecules in system B were separated from silica surface by heavy oil components which were adsorbed on silica
189
J. Sun et al. / Chemical Physics Letters 678 (2017) 186–191
Fig. 3. The snapshots of oil distribution at different simulation time. (a) System A; (b) system B. Asphaltenes, red; resins, green; light components, brown. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
16
Distance (nm)
4
(a) system A system B
(b) system A system B
3
Distance (nm)
20
12 8 4
2
1
0
0 0.0
0.5
1.0
1.5
2.0
2.5
3.0
t (ns)
0.0
0.5
1.0
1.5
2.0
2.5
3.0
t (ns)
Fig. 4. The distance between the COMs of oil phase and the nonmoving reference with time evolution. (a) Separate oil; (b) residual oil.
more notable (Fig. S3). Because heavy oil molecules have larger surface area than light oil molecules, the effects of water molecules was more pronounced. Herein, residual oil in system B can move gradually under external force and water scouring.
3.0 system A system B
2.4
3.3. Solvent accessible surface area
g (r)
1.8 1.2 0.6 0.0 0.0
0.4
0.8
1.2
1.6
2.0
r (nm) Fig. 5. RDFs of light oil molecules near silica nanopore surface averaged over the last 1 ns.
surface. They can move firstly under applied force. As a result, residual oil which adsorbed on silica surface was exposed to bulk water. In spite of much stronger nonbonded interaction between silica surface and heavy oil molecules in system B (Fig. S2), with time evolution during SMD simulation, effects of water molecules on heavy oil components in residual oil increased and became
Fig. 6a shows change of solvent accessible surface area (SASA), which is only effected by aggregation of heavy oil molecules (without consideration of influence of other components) while calculation, with time evolution. At the beginning of MD simulation in system A, heavy oil molecules aggregated together resulted in quick reduction of SASA. With most of the molecules had been assembled, the amplitude of SASA changes declined after 2.5 ns. Without the obvious adsorption of heavy oil on silica nanopore wall, it is negligible that the changes of SASA, which is caused by the influence of silica surface on the conformation of aggregates. In system B, same SASA-time curve can be found before 6 ns in MD simulation. With time evolution, heavy oil molecules move closer to silica surface and adsorbed on it under nonbonded interactions. As a result, part of the aggregates was destroyed. SASA increased in the last 4 ns because of the dispersion of heavy oil molecules. However, the heavy components were not dispersed completely. SASA at 10 ns is less than its initial value. In SMD simulation, the aggregates dispersed as individual molecules quickly under external force. Oil molecules distributed larger range along silica pore (Fig. S4). That led to rapid increase
2
2
Solvent accessible surface area (nm )
J. Sun et al. / Chemical Physics Letters 678 (2017) 186–191
Solvent accessible surface area (nm )
190
MD
(a)
190
SMD
system A system B
180 170 160 150 140 130
0
2
4
6
8
0
2
(b)
35
MD
SMD
system A system B
30 25 20 15 10 0
2
4
t (ns)
6
8
0
2
t (ns)
Fig. 6. Solvent accessible surface area of heavy oil molecules over the simulation time of 10 ns in MD simulation and 3 ns in SMD simulation. (a) Without affection of silica pore surface; (b) Reduction of solvent accessible surface area caused by adsorbing to silica pore surface.
of SASA in both of the simulation systems. The notable difference is that an extremal point appeared at 1.4 ns in system B. This is because of relatively strong interaction among heavy oil molecules resulted in reaggregation occurred in separate oil phase. Therefore, although separate oil had been pulled away from residual oil, SASA decreased after 1.4 ns. Fig. 6b reveals reduction of SASA caused by adsorption of heavy oil molecules onto silica pore surface. In system A, we found that SASA remained stable after slow increase in 2.5 ns MD simulation. That means only small part of heavy oil molecules, which located close to silica surface in initial configuration, was adsorbed on pore wall. Quite differently, SASA increased within initial 7.5 ns in system B. Then, it followed a slight decrease in last 2.5 ns. That denotes most of heavy oil molecules were adsorbed on silica surface. Finally, with the saturated adsorption, conformation change or desorption of individual heavy oil molecule led to the decrease of SASA. During the SMD simulation process, we can see that SASA increased invariably in system A. Separate oil was pulled away from initial position under external force. Its heavy oil components were adsorbed on silica surface gradually because of losing protection from residual oil. As a result, SASA that occupied by silica pore wall increased. In system B, a sudden decrease of SASA appeared at 1 ns and maintained about 0.4 ns, which indicates the desorptionadsorption process of heavy oil molecules in residual oil. Residual oil molecules desorbed from silica surface under external force and water scouring. With short distance displaced, they were adsorbed on silica surface again. Although heavy oil components were adsorbed on silica surface readily in system B, residual oil can be easier detached by external force compared to system A. Therefore, to enhance oil-displacement efficiency, it is significant to reasonably exploit the transportation of oil molecules inside silica nanopore, which will help to design efficient oil displacement agent. 4. Conclusion In this work, a series of MD simulations were performed to research the mechanism of oil migration inside a silica nanopore. With equilibrium MD simulations, the distributions of heavy oil components of crude oil inside silica nanopore were revealed. The results showed that heavy oil molecules dispersed well in toluene. However, due to insolubility in heptane, heavy oil components deposited and adsorbed on silica surface. The steered MD simulations were conducted to investigate oil transport process inside silica pore. Oil molecules in the central part of oil phase migrated firstly under external force. With pulling away from initial position and exposure to silica surface, they were adsorbed on silica pore wall due to nonbonded interactions. Residual oil in sys-
tem B is more easily pulled away from silica surface under external force than system A. The results provide insight vision of oil migration and were expected to promote the enhanced oil recovery. Acknowledgements We gratefully appreciate the financial support from the National Natural Science Foundation of China, China (21573130 and 11647013) and the Shandong Provincial Natural Science Foundation, China (ZR2016BP08).
Appendix A. Supplementary material Supplementary data associated with this article can be found, in the online version, at http://dx.doi.org/10.1016/j.cplett.2017.04. 057. References [1] M.O. Elsharafi, B. Bai, Effect of weak preformed particle gel on unswept oil zones/areas during conformance control treatments, Ind. Eng. Chem. Res. 51 (2012) 11547–11554. [2] R. Seright, Washout of Cr (III)-Acetate-HPAM Gels from Fractures, Presented at the SPE international symposium on oilfield chemistry, Houston, Texas, 5–7 February 2003; Paper 80200. [3] Q. Xue, Y. Tao, Z. Liu, S. Lu, X. Li, T. Wu, Y. Jin, X. Liu, Mechanism of oil molecules transportation in nano-sized shale channel: MD simulation, RSC Adv. 5 (2015) 25684–25692. [4] J. Zhong, X. Wang, J. Du, L. Wang, Y. Yan, J. Zhang, Combined molecular dynamics and quantum mechanics study of oil droplet adsorption on different self-assembly monolayers in aqueous solution, J. Phys. Chem. C 117 (2013) 12510–12519. [5] Q. Liu, S. Yuan, H. Yan, X. Zhao, Mechanism of oil detachment from a silica surface in aqueous surfactant solutions: molecular dynamics simulations, J. Phys. Chem. B 116 (2012) 2867–2875. [6] X. Hu, Y. Li, H. Sun, X. Song, Q. Li, X. Cao, Z. Li, Effect of divalent cationic ions on the adsorption behavior of zwitterionic surfactant at silica/solution interface, J. Phys. Chem. B 114 (2010) 8910–8916. [7] N.R. Tummala, L. Shi, A. Striolo, Molecular dynamics simulations of surfactants at the silica-water interface: anionic vs nonionic headgroups, J. Colloid Interface Sci. 362 (2011) 135–143. [8] T. Takanohashi, S. Sato, R. Tanaka, Structural relaxation behaviors of three different asphaltenes using MD calculations, Pet. Sci. Technol. 22 (2004) 901– 914. [9] O. Castellano, R. Gimon, C. Canelon, Y. Aray, H. Soscun, Molecular interactions between Orinoco belt resins, Energy Fuels 26 (2012) 2711–2720. [10] O. Mullins, S. Betancourt, M. Cribbs, F. Dubost, J. Creek, B. Andrews, L. Venkataramanan, The colloidal structure of crude oil and the structure of oil reservoirs, Energy Fuels 21 (2007) 2785–2794. [11] C. Pierre, L. Barré, A. Pina, M. Moan, Composition and heavy oil rheology, Oil Gas Sci. Technol. 59 (2006) 489–501. [12] P.M. Spiecker, K.L. Gawrys, C.B. Trail, P.K. Kilpatrick, Effects of petroleum resins on asphaltene aggregation and water-in-oil emulsion formation, Colloids and Surf. A 220 (2003) 9–27.
J. Sun et al. / Chemical Physics Letters 678 (2017) 186–191 [13] V.A.M. Branco, G.A. Mansoori, L.C.D.A. Xavier, S.J. Park, H. Manafi, Asphaltene flocculation and collapse from petroleum fluids, J. Pet. Sci. Eng. 32 (2001) 217– 230. [14] S. Mozaffari, P. Tchoukov, J. Atias, J. Czarnecki, N. Nazemifard, Effect of asphaltene aggregation on rheological properties of diluted athabasca bitumen, Energy Fuels 29 (2015) 5595–5599. [15] L. Goual, A. Firoozabadi, Measuring asphaltenes and resins, and dipole moment in petroleum fluids, AIChE J. 48 (2002) 2646–2663. [16] F. Gao, Z. Xu, G. Liu, S. Yuan, Molecular dynamics simulation: the behavior of asphaltene in crude oil and at the oil/water interface, Energy Fuels 28 (2014) 7368–7376. [17] G. Lv, F. Gao, G. Liu, S. Yuan, The properties of asphaltene at the oil-water interface: a molecular dynamics simulation, Colloids and Surf. A 515 (2017) 34–40. [18] J. Zhang, D. Tian, M. Lin, Z. Yang, Z. Dong, Effect of resins, waxes and asphaltenes on water-oil interfacial properties and emulsion stability, Colloids and Surf. A 507 (2016) 1–6. [19] H. Yan, S. Yuan, Molecular dynamics simulation of the oil detachment process within silica nanopores, J. Phys. Chem. C 120 (2016) 2667–2674. [20] J.G. Speight, Petroleum asphaltenes – Part 1 – Asphaltenes, resins and the structure of petroleum, Oil Gas Sci. Technol. 59 (2004) 467–477. [21] C. Pierre, L. Barre, A. Pina, M. Moan, Composition and heavy oil rheology, Oil Gas Sci. Technol. 59 (2004) 489–501.
191
[22] S. Izrailec, S. Stepaniants, B. Isralewitz, D. Kosztin, H. Lu, F. Molnar, W. Wriggers, K. Schulten, Steered Molecular Dynamics, Springer-Verlag, Berlin, 1998. [23] D. van der Spoel, E. Lindahl, B. Hess, The GROMCAS development team, GROMACS User Manual version 4.6.3, Royal Institute of Technology and Uppsala University, Stockholm and Uppsala, 2013. [24] H. Zhang, Y. Ma, Q. Hao, H. Wang, G. Liu, S. Yuan, Molecular dynamics study on mechanism of preformed particle gel transporting through nanopores: surface hydration, RSC Adv. 6 (2016) 7172–7180. [25] H.J.C. Berendsen, J.P.M. Postma, W.F. Van Gunsteren, J. Hermans. In: B. Pullman (Ed.), Intermolecular Forces, D. Reidel Publishing Company, Dordrecht, 1981. [26] M.D. Elola, J. Rodriguez, D. Laria, Liquid methanol confined within functionalized silica nanopores. 2. Solvation dynamics of coumarin 153, J. Phys. Chem. B 115 (2011) 12859–12867. [27] M.D. Elola, J. Rodriguez, D. Laria, Structure and dynamics of liquid methanol confined within functionalized silica nanopores, J. Chem. Phys. 133 (2010) 154707. [28] G. Bussi, D. Donadio, M. Parrinello, Canonical sampling through velocity rescaling, J. Chem. Phys. 126 (2007) 014101. [29] B. Hess, H. Bekker, H.J.C. Berendsen, J. Fraaije, Lincs: a linear constraint solver for molecular simulations, J. Comput. Chem. 18 (1997) 1463–1472. [30] U. Essmann, L. Perera, M.L. Berkowitz, T. Darden, H. Lee, L. Pedersen, A smooth particle mesh ewald method, J. Chem. Phys. 103 (1995) 8577–8593.