Molecular Dynamics approach for Fracture Simulation along a Weakly Bonded Interface

Opening mechanisms for fractures (joints) characterized by weak bonding (nano-micro scale cracks) in low-porosity rock play a key role in shale oil and shale gas development through staged hydraulic fracturing. We explore the cohesive interface strength of two weakly bonded slabs of Polymethylmethacrylate (PMMA) with nanoscale Molecular Dynamics (MD) methods, calibrated to experimental data and simulated by the Finite Element Method (FEM). The proper stress/strain state at the weakly bonded interface is determined because it is required for MD simulation. Then, we develop a representative PMMA structure as an input for LAMMPS software to simulate a tensile strength test. Results including per-atom stress values and pressures in different components are extracted, and we note that the fracture location and its behavior deduced from MD simulations follows experimental results. The total force for this simulation is close to 80 × 10 N, and we can use this parameter to reflect the fracture behavior of the weakly bonded PMMA slabs with an acceptable accuracy, considering the levels of uncertainty.

Fractured rocks, especially fractured carbonates and "shales", comprise hydrocarbon reservoirs of interest throughout the world. Enhanced Oil Recovery methods (EOR) or hydraulic fracturing procedures (e.g. for shale gas and shale oil stimulation) are needed to more economically access and release trapped hydrocarbons in most cases, leading to a need for better simulation of rock mass fabric effects (Hashemi et al. 2016). Most stimulation methods for tight sandstones and shales are based on induced hydraulic fracturing to form conductive pathways (Dmour & Shokir 2010) and to contact larger volumes of rock in the reservoir. Crack initiation occurs when the tensile stress at or close to a crack tip is greater than the "molecular cohesive strength of the material" (Sivakumar & Maji 2014). Static crack propagation in a reservoir domain relates to external stress/strain and internal pressurization conditions that can be explored to a degree in the laboratory by fracture mechanics methods and expressed quantitatively through analytical or numerical estimations of crack growth and specimen rupture (Veeranjaneyulu & Rao 2012).
Computational methods, including continuum and dis-continuum methods (Mitra 2006), now allow mathematical "experiments" of rock behavior under various states including local tension (Hart 2003). An ability to analyze fractured rock mass deformations after induced fracture is the major distinction between continuum and dis-continuum methods (Elmo et al. 2013). When there are few fractures, or when fracture opening is small, continuum approaches may suffice (Jing & Hudson 2002). Additionally, some experimental observations related to crack initiation and fracture propagation parameters at a microscopic scale cannot be explained by continuum methods (Huang et al. 2014). In such conditions, alternative methods are required.
Rock masses are characterized by flaws and fractures at all scalesfrom faults to grains to micro-cracks to crystal lattice defects. The main reason for fracture growth is the stress concentrations resulting from faults in grain boundaries and micro cracks (Mahbaz et al. 2013). The crack-formation and fracture behavior of polycrystalline materials are related to their atomic-bonds and grain boundary conditions (Abraham et al. 1997;Kim et al. 2001), which cannot be interpreted easily by continuum methods (Pei et al. 2013). However, useful atomic-scale knowledge about fracture propagation is obtained by discontinuum micro-methods (Rountree et al. 2002) such as Molecular Dynamics (MD) method, which gives insight into rock fracture behavior if implemented in appropriate states (Mosharafi et al. 2017).This well-known method determines the effects of external forces on particles with specific properties, based on calculating the reactions of every particle, such as those of the atoms in the crystalline structure of a material (Belak et al. 1995), or even the contact-law-based discrete element methods used in geomechanics.
In naturally fractured rock masses (Keshavarzi et al. 2014), joints form weak planar features in a strong rock matrix and hence are strongly favored to be the initial locations of fracture openings (Dehghan et al. 2015) and may also dominate the local tip pathway for induced fracture growth. The material model analogue to study in this MD investigation is therefore Mode I fracture of a weakly bonded planar interface in a strong matrix. Others have used transparent analogue specimens to emulate the fracture stimulation processes (Suarez-Rivera et al.  We apply the MD method to simulate an experimental Mode I fracture of a weak interface between two PMMA slabs. This interface was formed by controlled compressive stress-and temperature-induced partial bonding. Stable crack propagation is achieved through increasing the uniaxial compressive stress on the unconfined rectangular specimens, which had a central starter hole on the weak interface (Gomez et al. 2016) causing Mode I fracture propagation at the advancing tip. We simulate the experiment using continuum numerical models to give stress/strain conditions at the joint tip for implementation into a LAMMPS © MD model.

Simulation procedure and results
We , normal to the plane of the interface, leading to symmetric crack initiation in both directions. Considering the stress behavior from the FEM simulation, we focused on the SCZ (i.e., the bonding location of the two PMMA slabs (Fig. 3)), at molecular scale, and its behavior under an incrementally increasing tensile load.  To conduct the MD, a PMMA molecular structure was assembled through the following procedure:  One chain of PMMA ( Fig. 4.a) was built with 30 monomers, and the molecular structures optimized via a geometry optimization module with the Dreiding force field (Mayo et al. 1990).  The polymer structure was placed in an amorphous cell with dimensions 30×7×20.1 Å and a density of 1.18 g/cc, and again optimized with the Dreiding force field ( Fig.  4.b).
 To emulate the two lightly bonded PMMA slabs, two amorphous cells were positioned adjacent to each other, so the new cell thus created included both bonded PMMA slabs, and the geometrywas optimized once more.
The new cell was placed at a 422 K (149 ˚C) temperature and a pressure (stress) of 24 kPa applied to the surfaces normal to the joint trace with an NPT ensemble and Dreiding force field (Mayo et al. 1990) to create a partially annealed surface. However, the structure disintegrated because the impact of the temperature was greater than the restraining effect of the pressure. The annealing procedure was then changed to a NVT MD ensemble, the structure heated to 422 K (149 ˚C) (  The final PMMA structure was used by a LAMMPS © code to extend along the Z direction until separation. For the stretching process, a coarse-grained force field was used for proper interatomic interaction (Okada et al. 2000). An appropriate interatomic potential is essential for obtaining reasonable results in MD simulations, and arises from the bonded and non-bonded atomic interactions in the molecular structure of the system. The implemented interatomic potential can be expressed by the summation form shown in Eq. 1, which was used by Okada et al. (2000) in studying the polymer's molecular structure under tension.
Bond stretching and bending describe the interaction between two or three chemically bonded neighboring atoms, respectively (Kim & Lee 2002). Torsion angle (dihedral angle) expresses interactions among four connected atoms (Monasse & Boussinot 2014), and inversion (improper interaction) is applied to keep atoms in a plane (Baker 2017). The first four terms specify the potential among covalently bonded atoms, and the two final terms indicate the potential energy of non-bonded atoms.
The structure in Fig. 4.c was stretched in the Z direction, perpendicular to the bonded interface, at a velocity of 0.1 Å/femtosecond. The molecules had insufficient time to interact, and so were torn apart at the stretching jaw's location (Fig. 7.c & 7.d). The code was modified and used with an extensional velocity of 0.01 Å/femtosecond, with ambient temperature control described by the Nose-Hoover thermostat (Nose 1984;Hoover 1985). This thermostat worked well with the PMMA's molecular structure and the implemented atomic force field (Okada et al. 2000). At this velocity, all molecules had sufficient time to interact, so the structure successfully separated along the weakly bonded interface (Fig.  8). Upscaled stresses can be estimated through two approaches: (1) molecular formalism is based on the center of masses, velocities, positions, and integration of the total forces per molecule, and (2) atomistic calculation that is based on the masses, velocities, positions, and integration of the forces of every atom in the system (Cui et al. 1996). Per-atom stress cannot be directly interpreted from a physical perspective and the stresses should be considered only for a representative group of atoms, based on integration of forces acting on the atoms (Thompson et al. 2009). Alternatively, the atomistic stress tensor in MD can be calculated through the differentiation of strain energy of a specified volume of atoms by the strain tensors. Strain energy and strain tensors are evaluated through interatomic potentials and the distance between a pair of atoms, respectively (Saitoh et al. 2014).
Micromechanical behavior and interface tension are strongly influenced by local pressure (Heinz et al. 2005). The pressure of the system can be represented by a bulk pressure related to the stress values and considered as the average of the three normal principal stress components. In addition, the ideal gas law, using the density and temperature (Barisik & Beskok 2011), can be invoked to predict the pressure field. There are different methods to calculate the micromechanical local stress tensors, based on the same parameters listed above (Heinz et al. 2005).
The parameters computed within the LAMMPS © code at each time step during simulation of PMMA extensional straining include the per-atom stress in the stretch direction, total pressure, pressure component (stress) parallel to the extensional direction (Pzz), and one of the pressure's perpendicular components (Pyy) (Fig. 9). Valuable outcomes are obtained by considering the structure as a heterogeneous media in which stress can vary at different locations. Computing viral stresses on every atom is a common approach, which is conducted with regard to mass, velocity, characteristic volume, and force acting on individual atoms (Fenley et al. 2014). For calculating the per-atom stress, initially the stresses of all the atoms of the system were added together. The computed quantity was in units of pressure × volume (Li et al. 2015), and the volume over which the stress is averaged (simulation box volume) was considered to be the characteristic volume. Hence, the total stress was divided by system volume to obtain the stress per atom (Li et al. 2015) (Fig. 10).

Discussion
Stretching the structure affects the velocity, positions, and interactions among particles and thus affects the pressure (stress) values. Cartesian components of stress tensors are calculated based on the dimensions of the simulation cell. The components of these tensors can be different from the bulk pressure because of unbalanced stresses during the simulation (Miguel & Jackson 2006). The MD simulation shows total pressure and Pzz increasing slowly during the stretching process, while the Pyy was relatively constant and fluctuated only slightly around the zero line (Fig. 10). Fracture initiation of the PMMA slabs was observed at approximately 83 MPa at its maximum state (Gomez et al. 2016). The best MD outcome that could reflect this attribute was the total pressure, which was approximately 84 MPa at the opening time. Hence, the total pressure seems to reflect the real conditions in the experimental simulation more accurately than other parameters, and can be regarded as a proof for the validation of the simulation results. Validation of the MD simulation outcomes can be noted by considering some other issues such as implementing appropriate boundary conditions, ensembles, interatomic interaction functions, software for simulation, and relaxation time (Van Gunsteren & Mark 1998). For simulating the tensile test of PMMA as a homogenous polymer in a way that is analogous to a macroscopic experiment, the boundary conditions were considered periodic on every direction except the tensile direction (Makke 2011). Using these periodic boundary condition resulted in minimizing the system size effect (Shaw et al. 2016) and creating a larger system (Murray et al. 2010). The motion of the atoms in the loaded section was also integrated with an NVT ensemble, using the Nose-Hoover thermostat, which is the common approach in molecular simulations of uniaxial strength process (Zhao et al. 2018;Li et al. 2013;Makke 2011). Furthermore, the atomic force field was imported accurately from a valid research publication (Okada et al. 2000) into the LAMMPS software, which is a well-known and precise platform for conducting studies on the mechanical properties of polymers' molecular structures (Wu al. 2014;Yang et al. 2013). In addition, to ensure the appropriateness of the molecular structure, relaxation was continued for about 20 picoseconds after the balance time (Fig. 5 & Fig. 6). MD simulation results indicate an initial per-atom stress summation at zero strain, probably because of residual stresses from the annealing (bonding) process for the two PMMA slabs. This initial stress increases rapidly with strain until it reaches zero at a strain of about 0.02. Increasing the strain causes the total tension of all the atoms in Z direction (stretching direction) to increase, leading to instabilities, loss of the initial structure, and a fracture opening at the strain of about 0.81. As illustrated in Fig. 9, the MD simulation shows stable fracture propagation with axial stress (strain) increase, and of course, that is how the test procedure takes place. In addition, during the extensional process, a linear elastic behavior (with some negligible local fluctuations) is observed in the "peratom stress-strain" curve, and this result is logical considering the brittle nature of the bonds.
Atomic interactions such as those of Fig. 9 cannot be successfully resolved at high-speed extensional velocities. The lowest velocity that we considered was 0.01 Å/femtosecond, when all atoms had sufficient time to interact with each other. Reducing the velocity improved the simulation outcomes at some small computational expense. The real axial displacement velocity speed for the entire specimen test was 0.01 mm/min, equal to 1.66667×10 -12 Å/femtosecond. The extensional velocity is one of the most important parameters affecting the results, and it may even modify the final peak value for fracture opening, so different velocities should be explored to assess the effects of the stretching rate on per-atom stresses.
Two important issues related to this work should be noted. First, the molecular dynamics simulation was conducted considering periodic boundary conditions in the directions lateral to the stretching direction. Using these periodic boundary conditions allows covalent bonds to be formed between monomers inside the simulation cell and with their virtual images in other periodic boxes, so the polymer network continues to infinity (Yang et al. 2013). Although using periodic boundary conditions minimize the size effect (Shaw et al. 2016), changing the initial size of the molecular structure can still change its mechanical properties (Huang & Zhuo 2006). Determining the influence of changing the spatial scale of the initial molecular structure would be a useful focus for future research. Second, it would be helpful to develop such results into an upscaled continuum approach, such as a damage mechanics approach (Cheng & Dusseault 2004), in order to make problems with hundreds to thousands of intersecting natural interfaces computationally tractable (Yetisir et al. 2016). In other words, to be useful, results from the small scale must be homogenized in a suitable manner to give constitutive relationships that can be implemented at the large scale. For example, in the realm of rock mechanics, hydraulic fracturing and other processes at the reservoir scale are best analyzed with continuum models, in which the constitutive laws at the upper scale are derived from atomistic (discrete block) analyses at the lower scale, using the concepts of a representative volume.

Conclusion
We investigated the strength of a lightly bonded planar interface between two PMMA slabs under static load with the Molecular Dynamics method and compared outputs to experimental results. A Finite Element Method continuum model first provided stress concentrations and directions at the fracture initiation point as input for the MD simulation. The MD simulation was conducted to model the PMMA slabs via a coarse-grained force field. Different parameters were extracted from the simulation results, such as various components of pressure (stress), total pressure (stress), and per-atom forces of the structure. Our main conclusions are: 1. Extensional velocity (i.e. the loading rate) is one of the most significant parameters on fracture behavior emulation in MD. At high velocities, particles do not have enough time to interact with each other properly. 2. MD simulations indicated that the fracture initiated and propagated on the weakly cohesive interface of PMMA slabs, similar to the experimental results on annealed (lightly bonded) specimens. 3. Extracted total pressure from the MD simulations was 84 MPa at opening, reasonably close to experimental results. 4. The per-atom stress in PMMA kept increasing from the beginning of the extensional simulation until the time of fracture opening and even thereafter. This agrees with the stable fracture propagation experiments, where increases in axial compressive stress are needed to further extend the fracture in a stable configuration. 5. MD models can be used to simulate fracture opening and extract related parameters for brittle materials. 6. We conclude that extension of these results to brittle rock samples with polycrystalline structure and more heterogeneous properties is feasible to numerically study tensile crack initiation and propagation behavior, albeit with a substantial computational effort. 7. Other work we have done shows that this can lead to reasonable upscaled continuum mechanics simulations of a structure that at the small scale is comprised of individual interacting element. This also implies that MD simulations may help in making constitutive models at larger scales more rationally founded on the integrated behavior of the particles at the smaller scales.