A critical study of the parameters governing molecular dynamics simulations of nanostructured materials

Computational Materials Science 153 (2018) 183–199

Contents lists available at ScienceDirect

Computational Materials Science journal homepage: www.elsevier.com/locate/commatsci

A critical study of the parameters governing molecular dynamics simulations of nanostructured materials A.R. Aliana,1, S.A. Meguida, a

T

⁎

Mechanics and Aerospace Design Laboratory, Department of Mechanical and Industrial Engineering, University of Toronto, Toronto, Ontario M5S 3G8, Canada

A R T I C LE I N FO

A B S T R A C T

Keywords: Interatomic potential function Time steps Cut-off function Strain rate Edge effect Carbon nanotube Graphene Boron nitride sheet

Molecular dynamics (MD) simulations have been used extensively over the past two decades to determine the mechanical and physical properties of nanomaterials. However, the discrepancy between the reported results from these atomistic studies shadows the reliability of this computationally efficient technique. This inconsistency is attributed to the misuse and incorrect application of MD as evidenced by the arbitrary use of interatomic potentials, cut-off function parameters, strain rate, time increment, and domain size in the conducted simulations. In this paper, we highlight erroneous simulations by investigating the influence of these parameters on the elastic and fracture properties of nanostructured materials; including carbon nanotubes, graphene, and boron nitride (BN) sheets subject to direct and contact loads. The effect of interatomic potential type was investigated by comparing the predicted properties from AIREBO, Tersofff, CVFF, and ReaxFF potentials with those obtained with experimental and DFT techniques. The cut-off function parameters were also investigated to determine the optimum inner and outer cut-off radii selected to capture the actual physical behavior and avoid the reported strain hardening phenomena. Furthermore, MD simulations with strain rates spanning several orders of magnitudes and time increments ranging from 0.1 to 20 fs were performed to define the maximum allowable parameters for each material and loading scheme. Additionally, graphene and BN sheets with side length up to 500 Å were modeled to determine the size and edge effects on the mechanical properties. Finally, a set of parameters is recommended in each investigation to help guiding future atomistic studies obtaining reliable results using the available computational resources.

1. Introduction Molecular dynamics (MD) simulations have been used extensively over the last two decades for modeling nanomaterials [1]. However, several parameters influence the obtained properties from MD simulations including interatomic potential type, cut-off function radii, time step, strain rate, and domain size [2,3]. Determining quantitatively the effect of each parameter on material properties is very critical in achieving reliable results. Different interatomic potentials and force fields have been used over the past two decades to study CNTs, graphene, and their composites [4–7]. As a result, the reinforcement effect of CNTs and graphene in their nanocomposite might be underestimated or overestimated depending on the interatomic potential used to model the system. Furthermore, all force fields that use harmonic function to model atomic bond such as CVFF force field do not permit bond formation and breakage, which limits the simulation to only small deformation [8–10]. On the other hand, other reactive interatomic potentials such as AIREBO and ReaxFF ⁎

1

potentials can be used to model bond breakage and hence allowing studying large deformation and fracture mechanisms [11–14]. Another source of the contradictions in the reported mechanical properties the literature is the selection of the cut-off function parameters in Tersoff and AIREBO potentials [15,16]. For example, the original minimum and maximum radii in the cut-off function for the covalent bonds in AIREBO potential are 1.7 Å and 2.0 Å, respectively. However, these values cause unphysical strain hardening in the structure at high strains. To address this problem, researchers varied the cut-off function parameters to achieve properties that are in good agreement with those obtained with experimental and DFT techniques [17,18]. Selecting the appropriate time step in MD simulations is very critical for obtaining reliable results with a reasonable computational cost, achieving energy conservation, and preventing the MD system from exploding [19]. Different time steps ranging from 0.1 to 15 fs were used in previous MD simulation of CNTs, graphene, and BN structures [20–22]. Budarapu et al. [23] have studied the effect of time step on crack initiation and growth mechanisms in Graphene using MD

Corresponding author. E-mail address: [email protected] (S.A. Meguid). On Leave from Mechanical Design and Production Dept., Faculty of Engineering, Cairo University, Giza 12613, Egypt.

https://doi.org/10.1016/j.commatsci.2018.06.028 Received 6 March 2018; Received in revised form 14 June 2018; Accepted 18 June 2018 0927-0256/ © 2018 Elsevier B.V. All rights reserved.

Computational Materials Science 153 (2018) 183–199

A.R. Alian, S.A. Meguid

MD simulations were compared with recently obtained numerical and experimental findings [15,39,41–51]. Decidedly, the findings of this work should help in understanding the source of anomalies and inconsistency among the reported data in the literature and also in finding the appropriate parameters for each specific loading and material type in future studies. This paper is organized as follows. Following this introduction, the methodology and the objectives of this study are outlined in Section 2. The results of the conducted MD simulations for the considered materials are presented and discussed in Section 3. Finally, the main contributions and conclusions are summarized in Section 4.

simulations. Their results showed a significant dependence of post yielding behavior of Graphene on time step and a time step of 0.1 fs was required to accurately monitor crack initiation and propagation in the fracturing sheet. However, such a small time step will require an enormous computational cost to model large systems for longer periods. Conducting MD simulations at very high strain rates means that the molecular structure does not have enough time to respond (equilibrate) to the applied load [24]. The effect of strain rate in MD simulation on the predicted mechanical properties were investigated in the literature to determine a suitable range of loading rates that gives reliable results. However, the reported results of these investigations were contradicting, with some papers indicating a significant effect of the strain rate on all mechanical properties, while others showing slight to no effect [25]. For example, Liu et al. [26] found that the tensile strength of CNTs is strongly affected by the applied strain rate [27]. On the other hand, Wei et al. [28] found that the tensile yield strain of a SWCNT was slightly influenced by the strain rate. Zhang et al. [29] reported that the strain rate used in MD simulations has a negligible impact on the predicted Young’s modulus of graphene. Zhao and Aluru [30] showed that the influence of strain rate on the fracture stress of graphene was only significant at elevated temperatures. Han et al. [31] studied the effect of strain rate on the mechanical properties of BN sheets via MD simulations. Their results showed that Young’s modulus decreases while fracture stress and strain increase with increasing the strain rate. A contrast relationship between mechanical properties of BN sheets and strain rate was reported by Mortazavi et al. [32]. Their results showed no dependence between Young’s modules and the strain rate and a very slight dependence between tensile strength and the loading rate. The effect of graphene and BN sheet size on the mechanical properties were also investigated by several research groups [33–37]. However, the reported results of these atomistic studies are contradicting, with several studies showing a growing weakening effect with increasing the size [37,38], while other concluding the opposite [39]. The influence of the sheet size originates from the presence of free edges and the associated compressive forces due to dangling bonds in edge atoms [40]. For example, Zhao et al. [39] studied the size effect of the elastic properties of graphene using MD simulations. Their results showed that Young's modulus of graphene increases with the increase of the sheet size. Another MD study concluded that the elastic properties of graphene increase with the sheet size until reaching a stable point (bulk level) at a width of 4 nm [22]. On the other hand, an MD study by Qiang et al. [38] indicated that the elastic and fracture properties of graphene decrease with the increase in the sheet size and the properties reach their bulk level at a sheet width of 6 nm. This discrepancy in the reported results were attributed to the difference in the interatomic potential and temperature in each study [22]. In this paper, a comprehensive MD investigation is carried out to examine the effects of interatomic potential type, cut-off function parameters, time step, strain rate, and domain size on the elastic and fracture properties of CNTs, graphene, and BN sheets. Atomistically, the impact of these parameters on the properties depends on the testing technique adopted. This is due to the distinct behavior and fracture mechanism associated with each loading type. To address this issue, numerical compression and tension tests were performed on SWCNT, while numerical tensile and nanoindentation tests were carried out on graphene and BN sheets. A thorough investigation of the effect of the type of the interatomic potential on the mechanical properties of graphene and CNTs is presented by considering the following interatomic potentials: CVFF, ReaxFF, Tersoff, and AIREBO. A wide range of strain rates (ε ̇ = 2 × 10−4 − 16 ns−1) and time steps (t = 0.1–20 fs) have also been investigated to determine the critical values for each material and testing technique; ensuring that we capture the material behavior without altering the physical nature of the process. Graphene and BN sheets with dimensions up to 500 Å × 500 Å have been modeled to study the effect of domain size (i.e. edge effect) on the material strength and fracture mechanism. The predicted properties from the conducted

2. Methodology Molecular dynamics is a well-established method typically used to compute the time evolution of a system of interacting atoms by calculating the interaction force between them and numerically integrating their kinetics using Newton’s law of motion [52]. The interaction forces between the atoms are calculated from the gradient of the interatomic potential. Several parameters influence the ability of MD simulations to accurately predict the material behavior, such as interatomic potential type, cut-off function radii, strain rate, time step, and domain size (edge effect). We conducted comprehensive MD simulations to study the effect of the aforementioned parameters on the elastic and fracture properties of CNTs, graphene, and BN sheets. Fig. 1(a)–(f) show the MD models and the applied boundary conditions used in our simulations. Both uniaxial tensile and indentation simulations were performed for graphene (Fig. 1(a) and (c)) and BN sheet (Fig. 1(b) and (d)), while uniaxial tension and compression simulations were conducted for SWCNTs (Fig. 1(f)). A pyramid-shaped diamond tip was used for the nanoindentation simulations, as shown in Fig. 1(e). The height and base dimensions of the indenter were selected to be 34 Å and 24 Å × 24 Å, respectively. The indenter tip contained 9 carbon atoms (3 atoms × 3 atoms). All MD simulations were performed with large-scale atomic/molecular massively parallel simulator LAMMPS [53] using one of the following interatomic potentials: the adaptive intermolecular reactive bond order (AIREBO) potential [34], Tersoff-type empirical interatomic potential [54], the consistent valence force field (CVFF) [55], and the reactive force field (ReaxFF) [13]. A detailed description of each potential function is presented in the Supplementary Material. The cut-off function parameters in AIREBO and Tersoff were varied to investigate their effects on the mechanical properties. The equations of motion was integrated using the velocity–Verlet algorithm [56] with time steps ranging from 0.1 fs to 20 fs. The energy of the initial structure was minimized using the conjugate gradient algorithm to obtain CNTs, graphene, and BN sheets with optimized configurations. The minimized structures were considered to be optimized once the change in the total potential energy (PE) of the system between subsequent steps is less than 1.0 × 10−10 kcal/mol [57,58]. Subsequently, MD simulations were performed in the NVT ensemble for 105 time-steps to equilibrate the minimized structures. The tensile simulations were then performed with strain rates ranging from ε ̇ ranging from 2 × 10−4to ∼ 16 ns−1, while the nanoindentation simulations were performed at different loading rates (indenter speed ranging from 0.001 to 10 Å/ps). The procedures of the uniaxial tensile and nanoindentation simulations are presented in detailed in our previous studies [18,59,60]. Table 1 summarizes the investigated parameters and the conducted MD simulations in the current study. A detailed discussion of each investigation is presented in the following sections. 3. Results and discussion In this Section, we present and discuss the results of the conducted MD simulations of the effect of the aforementioned parameters on the 184

Computational Materials Science 153 (2018) 183–199

A.R. Alian, S.A. Meguid

Fig. 1. MD models and applied boundary conditions used in the current study: (a) Graphene sheet of 100 Å × 100 Å under uniaxial tensile, (b) BN sheet of 100 Å × 100 Å under uniaxial tensile, (c) Graphene sheet 100 Å ×100 Å under nanoindentation, (d) BN sheet 100 Å × 100 Å under nanoindentation, (e) Pyramidshaped diamond indenter, and (f) Armchair (6,6) CNT of length of 65 Å under tensile loading.

Table 1 List of the considered parameters and the corresponding MD simulations conducted in the current study. Study type

Investigated parameter

Testing method

Cut-off distance

AIREBO potential

Tensile test of CNT Tensile test of graphene Indentation of graphene Tensile test of BN sheet Indentation of BN sheet

Tersoff potential Strain rate

Quasi static vs dynamic analysis (ε ̇ = 2 × 10−4−16ns−1)

Tensile test of CNT Tensile test of graphene Indentation of graphene Tensile test of BN sheet Indentation of BN sheet

Domain size

Lgraphene_side = 50–500 Å LBN_side = 50–500 Å

Tensile test of graphene Tensile test of BN sheet

Time step

τ = 0.1–20 fs

Tensile test of a CNT Tensile test of graphene Tensile test of BN sheet

Interatomic potential type

CVFF, AIREBO, Tersoff, ReaxFF

Tensile test of a CNT Tensile test of graphene

185

Computational Materials Science 153 (2018) 183–199

A.R. Alian, S.A. Meguid

Fig. 2. Stress-strain curves obtained from the conducted MD simulation of a uniaxial tensile test of an (5,5) armchair SWCNT of length 65 Å using different interatomic potential functions.

Fig. 3. Stress-strain predictions obtained from the conducted MD simulation of a uniaxial tensile test in the armchair direction of a 100 Å × 100 Å graphene sheet using different potential functions.

elastic and fracture properties of CNT, graphene, and BN sheet. The obtained properties are compared with those obtained in the literature using both experimental and numerical techniques. Additionally, a recommended set of parameters is provided in each investigation to help guiding future atomistic studies obtaining reliable results using the available computational resources.

were defined by the large sudden drop in the stress. The results indicate that Young’s modulus of the CNT moderately depends on the potential type. On the other hand, the fracture stress and strain are significantly impacted by that selection. For example, Young’s modulus ranged from ∼864 GPa for AIREBO potential to ∼1111 GPa for Tersoff potential. The predicted Young’s modulus for ReaxFF and CVFF force fields are 857 and 922 GPa, respectively. It is very clear from the corresponding stress-strain curves of Fig. 2 that CVFF force field fails to simulate the fracture stage of CNTs due to the harmonic assumption of covalent bonds. However, this force field is very efficient in modeling large polymeric systems and can be used to model the elastic behavior of CNT-reinforced composites [9,57]. The stress-strain curve obtained with ReaxFF force field is almost linear until fracture, which predicts a higher fracture stress (∼157.6 GPa) than those obtained experimentally [41,42] and numerically [43,44]. We recommend using ReaxFF force field if the interest is the prediction of the elastic behavior of the nanostructured material. Table 2 summarizes all results obtained in the course of this investigation and previous numerical and experimental investigations [41–44].

3.1. Effect of potential function type Interatomic potential is a mathematical formulation for calculating the potential energy of interacting atoms at given positions in space [17,61]. Potential functions are either empirical in nature with parameters obtained experimentally using X-ray and electron diffraction or numerical using ab initio, or semi-empirical, using quantum mechanics calculations [17,62]. Different potential functions and force fields were developed over the past decades to simulate different materials [63]. We studied the effect of potential function type on Young’s modulus, fracture stress, and strain of CNT and graphene by comparing the output results of MD simulations based on AIREBO potential, Tersoff potential, ReaxFF force field, and CVFF force field.

3.1.2. Effect of potential function type on mechanical properties of graphene The effect of potential function type on the predicted mechanical properties of graphene was also investigated. Fig. 3 shows the stressstrain curves obtained from uniaxial tensile MD simulations of a 100 Å × 100 Å graphene using different interatomic potentials. Table 3 summarizes the obtained elastic and fracture properties for all cases considered. The results indicate that the elastic modulus moderately

3.1.1. Effect of potential function type on mechanical properties of a SWCNT The stress-strain curves for an armchair SWCNT of (6,6) chirality and length of 65 Å are presented in Fig. 2. We determined the elastic modulus from each simulation from the slope of the initial part of the stress-strain curve (ε = 2%), while the fracture strain and fracture stress

Table 2 Mechanical properties of a SWCNT obtained from the literature and the present work for different interatomic potentials. The elastic and fracture properties of the (5,5) armchair nanotube of a length of 65 Å were obtained in this work using uniaxial tensile MD simulation. Interatomic potential type

Young’s modulus (GPa)

Fracture stress (GPa)

Fracture strain (%)

CVFF ReaxFF AIREBO (R1 = R2 = 2.0 Å) Tersoff (R = 2.0 Å, D = 0.0 Å) Ab initio [43] Finite element method [44] Experimental [41] Experimental [42]

921.9 856.7 863.8 1101.4 750–1180 939.1 320–1470 270–950

– 157.6 99.8 109.3 – 94.0 10–52 11–63

– 16.46 18.55 17.67 – 16.4 5.3 12.0

186

Computational Materials Science 153 (2018) 183–199

A.R. Alian, S.A. Meguid

Table 3 Mechanical properties of graphene obtained from the literature and the present work for different interatomic potentials. The elastic and fracture of a 100 Å × 100 Å graphene sheet were obtained in this work using a uniaxial tensile MD simulation. Interatomic potential type

Young’s modulus (GPa)

Fracture stress (GPa)

Fracture strain (%)

CVFF ReaxFF AIREBO (R1 = R2 = 2.0 Å) Tersoff (R = 2.0 Å, D = 0.0 Å) MD simulations [45] MD simulations [39] Ab initio [46] Ab initio [47] Experimental [48]

968.9 948.8 847.1 969.7 890.0 980–1040 1024.0 1050.0 900.0–1100.0

– 81.8 86.3 87.0 105.0 96.0 103.0 107.0–121.0 120.0–140.0

– 8.9 11.95 11.90 17.0 13.0 18.5 19.4 25.0

cut-off function parameters (R1 = R2 = 2.0 Å) predict very comparable results (Fracture stress = ∼87 GPa and fracture strain = ∼11.9%) and are in a good agreement with the reported values in the literature (see Table 3). The predicted fracture strain and fracture stress with ReaxFF force field are 25% and 5% smaller than those obtained with AIREBO and Tersoff potentials.

depend on the potential function with AIREBO potential predicting the lowest elastic modulus of 847 GPa and Tersoff and CVFF predicting the highest elastic modulus of ∼970 GPa. On the other hand, the fracture properties are significantly influenced by the potential type. Furthermore, all potentials succeeded in capturing the fracture behavior except CVFF force field. Both AIREBO and Tersoff potentials with modified

Fig. 4. Effect of cut-off radii upon behavior of a graphene sheet: (a) a reduced structure of graphene with inner and outer cut-off radii, (b) variation of the bond strength with the radial distance, and (c) strain hardening of a graphene sheet under uniaxial tensile simulation.

Fig. 5. Effect of cut-off radii in BN sheet: (a) A reduced structure of BN sheet with inner and outer cut-off radii. (b) Variation of the bond strength with the radial distance. (c) Strain hardening of BN sheet under uniaxial tensile simulation. 187

Computational Materials Science 153 (2018) 183–199

A.R. Alian, S.A. Meguid

250

3.2. Effect of cut-off radii in AIREBO and Tersoff potentials The purpose of the cut-off function in AIREBO and Tersoff potentials is to limit the bonded interactions to the nearest neighbouring atom (first-neighbour shell) [64,65], as shown in Figs. 4 and 5. Selecting the appropriate cut-off function parameters is very critical in obtaining reliable results and avoiding unrealistic hardening behavior observed at higher strains near the fracture point [47,48,66] (see Fig. 4(c) and (c)). We carried out comprehensive MD simulations to determine the effect of cut-off function parameters in AIREBO and Tersoff potentials on the mechanical properties of CNT, graphene, and BN sheet. The obtained properties are compared with the experimental data in the literature to determine the optimum cut-off function parameters for each material and loading condition. The original cut-off function f (rij ) of REBO potential component in the AIREBO potential considers a gradual reduction in the interaction between atoms as follows:

⎧ ⎪

rij < R(1)

1 (1)

Uniaxial stress (GPa)

200

100

50

0 0

(1)

0.5

0.6

230

(2)

where R = R(1) = R(2) and its value ranges from 1.7 to 2.0 to determine the optimum cut-off function parameters for graphene and CNTs. The cut-off function fc in Tersoff potential also has a significant effect on the predicted fracture properties of BN sheets. Fig. 4 shows the stress-strain curve of BN sheet under uniaxial strain and it is very clear that the predicted material behavior using the original cut-off function leads to unphysical strain hardening at higher strains [67]. The original cut-off function in Tersoff potential is given by

1 rij < (R−D) ⎧ ⎫ ⎪1 1 ⎪ π (r − R) fc (rij ) = 2 − 2 sin ⎡ 2D ⎤ (R−D) < rij < (R + D) ⎣ ⎦ ⎨ ⎬ ⎪ ⎪ 0 rij > (R + D) ⎩ ⎭

0.2 0.3 0.4 Tensile strain

250

Fracture stress (GPa)

1 rij < R ⎫ f (rij ) = ⎧ ⎨ ⎩ 0 rij > R ⎬ ⎭

0.1

(a)

where R(1) is the inner cut-off radius and is equal to 1.7 Å and R(2) is the outer cut-off radius and is equal to 2 Å, as shown in Fig. 4. In order to eliminate the strain hardening near fracture, we considered R(1) to be equal to R(2) which yields a modified cut-off function in the form of a step function:

210 190 170 150

130 110 90

1.7

1.75 1.8 1.85 1.9 1.95 Cut-off minimum radius R1 (Å)

2

(b) (3)

Fig. 6. Tensile MD simulations of (19,19) armchair CNT of length 77 Å using different R1 values in the cut-off function in AIREBO potential: (a) Stress-strain curves and (b) Variation of the fracture stress with R1.

The parameters R and D are carefully chosen to only include interactions with first-neighbour and their values range from 1.9 to 2.3 Å, and from 0.05 to 0.1 Å, respectively [68]. In order to eliminate the strain hardening near fracture, we considered D = 0 which yields a modified cut-off function in the form of a step function:

1 rij < R ⎫ fc (rij ) = ⎧ ⎨ ⎩ 0 rij > R ⎬ ⎭

150

⎫ ⎪

π (rij − R ) f (rij ) = 1 + cos ⎡ (2) (1) ⎤ 2 R(1) < rij < R(2) ⎬ ⎨ ⎣ R −R ⎦ ⎪ ⎪ (2) < r 0 R ij ⎭ ⎩

R1=1.7 Å R1=1.80 Å R1=1.9 Å R1=1.92 Å R1=1.95 Å R1=1.97 Å R1=2.0 Å

CNT experiences a strain hardening behavior for R1 values below 1.97 Å. The fracture stress and strain decrease with increasing the cutoff inner radius, as shown in Fig. 6(b). Fig. 7(a)–(c) show snapshots of the deformed CNT before and after fracture. It is very clear that the maximum elongation and the fracture mechanism is significantly different in each case. We recommend conducting MD tensile simulations of CNTs and CNT reinforced composites with R1 equal to 2.0 Å. We also found out that the cut-off function parameters do not have any impact on the CNT behavior during compression test.

(4)

We varied R in our analysis from 1.9 to 2.3 Å to determine the optimum cut-off function parameters for BN sheets. 3.2.1. Effect of cut-off function on mechanical properties of CNT Uniaxial tension and compression simulations of an armchair SWCNT of (19,19) chirality and length of 77 Å were conducted for inner radius R1 ranging from 1.7 Å to 2.0 Å in the cut-off function in AIREBO potential. Fig. 6(a) shows the stress-strain curves of the CNT under tension for different values of R1. The elastic portion of the curves up to 20% strain is identical for all R1 values. However, at higher strains, the

3.2.2. Effect of cut-off function on mechanical properties of graphene A series of MD simulations were carried out to determine the effect of cut-off function parameters in AIREBO potential on the elastic and fracture properties of a 100 Å × 100 Å graphene sheet subjected to

188

Computational Materials Science 153 (2018) 183–199

A.R. Alian, S.A. Meguid

Fig. 7. Snapshots of an (19,19) armchair CNT of length 77 Å under tension before and after fracture for different inner cut-off radii: (a) R1 = 1.7 Å, (b) R1 = 1.92 Å, and (c) R1 = 2.0 Å.

inner cut-off radius. Fig. 11(a)–(c) show that the cut-off distance influences the fracture behavior of graphene.

uniaxial tensile and nanoindentation tests. Fig. 8(a) shows the indentation force vs the indentation depth for cut-off minimum radii R1 ranging from 1.7 Å to 2.0 Å. The initial segments of the force-indentation depth curves are identical for all cut-off distances. However, the maximum indentation force at the onset of fracture is significantly influenced by R1, as shown in Fig. 8(b). The maximum indentation force increases from 300 nN to 645 nN when R1 is increased from 1.7 Å to 1.96 Å before dropping to 137 nN at R1 = 1.97 Å and dropping again to 38 nN when R1 reached 2.0 Å. Fig. 9(a) and (b) show snapshots of the MD simulations post the complete fracture of the graphene sheet for R1 equal to 1.92 Å and 2.0 Å, respectively. The graphene sheet exhibits significant elongation prior to fracture in case of R1 = 1.92 Å, while showing a brittle fracture behavior for R1 = 2.0 Å. It is very clear that using R1 = 1.92 Å which was recommended in several studies [69,70] leads to 20 times higher residence force when compared to R1 = 2.0 Å which was also recommended in several studies in the literature [71–73]. The effect of cut-off function in AIREBO potential on the material behavior of graphene under uniaxial tensile simulation is also investigated. Fig. 10(a) shows the stress-strain curves obtained from tensile simulations for different values of R1 ranging from 1.7 to 2.0 Å. The graphene sheet experiences unrealistic strain hardening and very high fracture stress (140 GPa) and strain (40%) for R1 values lower than 1.97 Å (see Fig. 10(b)). The fracture stress and strain fluctuate for R1 values ranging from 1.97 to 2.0 Å. For example, the fracture strain reaches 18.4% for R1 = 1.97 Å prior to dropping to 10.4% for R1 = 1.99 Å, and finally increases to 12.5% for R1 = 2.0%. It is very clear that the fracture behavior of graphene is very sensitive to the

3.2.3. Effect of cut-off function on mechanical properties of BN sheet A series of MD simulations were carried out to determine the effect of cut-off function parameters in Tersoff potential on the elastic and fracture properties of a 100 Å × 100 Å BN sheet under nanoindentation and uniaxial tensile tests. Fig. 12(a) shows the variation of the indentation forces with the indentation depth for cut-off radii R ranging from 1.9 Å to 2.3 Å. The initial segments of the force-indentation depth curves are almost identical for all cut-off distances. However, the maximum indentation is significantly influenced by the R value force at the onset of fracture, as shown in Fig. 12(b). The maximum indentation force fluctuates around 45 nN for R ranging from 1.9 Å to 2.2 Å before dropping to 31 nN at D = 2.25 Å. We also investigated the effect of the cut-off function on BN sheet subjected to tensile loading. Fig. 13(a) shows the stress-strain curves of the BN sheet for different cut-off radii ranging from 1.9 Å to 2.3 Å. Our results reveal that the BN sheet elastic modulus is independent of the cut-off distance and is ∼610 GPa. However, the fracture strain and the fracture stress are significantly influenced by R value. Our research also shows that the respective fracture strain and strength fluctuate around 11.5% and 62 GPa. The strain hardening due to the cut-off function parameters can alter the BN sheet behavior and leads to significantly different fracture properties. The large variations of the fracture properties (35% in the fracture strain and 28% in the fracture stress) with the cut-off radius in Tersoff potential clearly demonstrate the critical role that this parameter plays in achieving reliable results.

189

Computational Materials Science 153 (2018) 183–199

A.R. Alian, S.A. Meguid

3.3.1. Effect of graphene sheet size on the mechanical properties The potential energy of the graphene sheet increases with increasing the sheet size due to the decreasing percentage of atoms located at the sheet edges. This increase in potential energy indicates the structural stability of the graphene sheet. Fig. 14(a) shows that the potential energy/atom rapidly increases with increasing the sheet size from 25 Å to 100 Å and then stabilizes to reach its bulk level value 7.37 eV/atom. Fig. 14(b) shows the variation of Young’s modulus and the fracture stress with the sheet size obtained from uniaxial loading along armchair direction. Unlike previous studies [22,38] that investigated the size effect by studying relatively small graphene sheets with side length up to 50 Å, we continued our investigation to include larger sheets. We found out that the elastic and fracture properties of graphene start to become independent of the sheet size when the side length reached 200 Å. For example, the elastic modulus increased from 730 GPa to 845 GPa, when the sheet side length increased from 25 Å to 300 Å. On the contrary, the fracture stress dropped from 93 GPa to 76 GPa, when the sheet side length increased from 25 Å to 300 Å.

Dcc=1.70

600 Indentation Force (nN)

Dcc=1.80 Dcc=1.90

500

Dcc=1.97 Dcc=2.00

400 300 200 100 0

0

10

20 30 40 50 Indentation time (ps)

60

70

(a)

3.3.2. Effect of BN sheet size on the mechanical properties The potential energy/atom and Young’s modulus of BN sheet have similar dependency on the sheet size, as shown in Fig. 15(a). They increase sharply with the initial increase in the BN sheet length from 25 Å to 100 Å. The potential energy is increased by 3.6% from 7.13 to 7.39 eV/atom, while Young’s modulus increased by 48% from 396 GPa to 586 GPa. This large increase in Young’s modulus is caused by the corresponding reduction in the percentage of edge atoms to the total number of atoms in the sheet from 15.8% to 4.2%. The potential energy and Young’s modules continued increasing but with a lower rate and finally reached the bulk level: 7.46 eV/atom and 660 GPa, respectively, for a sheet length of 200 Å. On the contrary, the fracture stress and the fracture strain decrease with increasing sheet size (see Fig. 15(b)). For example, the fracture stress decreased from 82 GPa to 51.5 GPa when the sheet length increased from 25 Å to 500 Å, while the fracture strain decreased from19.5% to 8.6%. These results indicate that the BN sheet behavior and fracture mechanism are significantly influenced by the sheet size due to the edge effect. Table 4 summarizes all results obtained in the course of this investigation as well as previous numerical and experimental studies of the mechanical properties of BN sheets [41–44]. The large discrepancy between the predicted fracture properties in our study and those predicted by MD simulations in Ref. [15] is caused by the difference in the selected cut-off parameters in the Tersoff potential and the sheet size. These results clearly indicate that the sheet size plays a significant role in obtaining reliable results. We recommend modeling BN sheets with width over 200 Å to eliminate the edge effect and enable the accurate prediction of the bulk properties.

Maximum Indentation Force (nN)

700

600 500

400 300

200 100

0 1.7

1.75

1.8 1.85 1.9 Cut-off distance (Å)

1.95

2

(b) Fig. 8. MD indentation simulations of a 100 Å × 100 Å graphene sheet with pyramid-shaped diamond tip using different R1 values in the cut-off function in AIREBO potential: (a) Stress-strain curves and (b) Variation of the maximum indentation force with R1.

3.4. Effect of strain rate/indentation speed on mechanical properties The strain rate at which the MD simulation is conducted significantly affects the output results [19]. Due to the computational limitations, MD simulations are always conducted at several orders of magnitude higher strain rates than those used in actual experimental tests [24]. Both quasi static and dynamics loading techniques were used in previous MD simulations [59,60]. In the quasi static simulations, an incremental strain (deformation) is applied to the system and the whole structure is then left to relax to reach equilibrium condition and allow the deformation to uniformly distributed though the MD system [75]. This loading step is then repeated until reaching the targeted strain and the strain rate is calculated as the change of strain over the MD

3.3. Effect of domain size on mechanical properties Finally, we investigated the effect of the sheet size and edge effect in MD simulations on the deformation and fracture mechanisms of graphene and BN during uniaxial tensile tests. The percentage of dangling bonds in the unpaired atoms located at the edges of these structures will impact the material strength its structural stability. A square graphene and BN monolayers with diagonal lengths up to 700 Å were considered in our study of the sheet size effect on the elastic and fracture properties. We adjusted the loading rate in each case in order to conduct the simulation with the same strain rate.

190

Computational Materials Science 153 (2018) 183–199

A.R. Alian, S.A. Meguid

(a) R1= 1.92 Å

(b) R1= 1.92 Å Fig. 9. Snapshots of the graphene sheet after fracture: (a) R1 = 1.92 Å and (b) R1 = 2.0 Å.

equilibrated during the loading process. Fig. 16(a) shows the stressstrain curves under uniaxial loading for different strain rates. The elastic modulus was found to decrease with increasing the strain rate. The fracture stress and strain are significantly influenced by the strain rate of the test, as shown in Fig. 17(a) and (b). The fracture stress and the fracture strain decreased from 101.7 GPa and 19.4% to 95.8 GPa and 17.3%, respectively, when the strain rate decreased from 16 ns−1 to 0.016 ns−1. The relationship between the fracture stress and the strain rate obtained with the dynamics scheme can be described using a logarithmic function (Fig. 16(b)). Higher strain rates lead to higher fracture stress and strain because the system atoms at lower loading rates have enough time to equilibrate and hence have a higher probability to exceed the bond breaking barrier [29]. Furthermore, at very strain rates higher than 10 ns−1, the CNT atoms cannot properly equilibrate and the potential energy has large fluctuation amplitudes. Table 5 summarizes the obtained fracture properties from all MD simulations for all considered strain rates.

relaxation/equilibration time [28]. In the dynamic analysis, a continuous load is applied on the molecular structure by deforming the structure with a constant strain rate. The potential energy, stresses, and forces are calculated throughout the simulation and used to determine the material properties. The key factor in obtaining reliable results with the dynamics/continuous deformation is loading the system with a strain rate low enough to allow the system to respond to the applied load (i.e. conduct the heat that is formed due to the extremely high loading rate [20]). In this investigation, deformation simulations were performed with strain rates spanning three orders of magnitude to determine the maximum allowable strain rate for CNT, graphene, and BN sheet. 3.4.1. Effect of strain rate on mechanical properties of CNT Uniaxial tension simulations of an armchair SWCNT of (6,6) chirality and length of 66 Å were conducted with an inner cut-off radius 2.0 Å in the AIREBO potential using 0.5 fs time step. Both dynamic and quasi static loading conditions were investigated to determine the appropriate loading method in MD simulations. The quasi static loading scheme was applied by displacing one end of the CNT by 0.125 Å while fixing the other end, followed by 10 ps (20 × 103 steps) of equilibration process in the NVT ensemble to allow the nanotube atoms between the two ends to deform and reach equilibrium positions. This loading was repeated until the CNT was completely fractured. In the dynamics loading scheme, one end of the CNT is displaced at a constant speed (strain rates ranging from 0.016 to 16 ns−1) while the system atoms are

3.4.2. Effect of indentation speed on resistance of graphene We investigated the effect of the deformation rate on the graphene resistance to indentation by conducting nanoindentation simulations of 100 Å × 100 Å graphene sheet with indentation speeds ranging from 0.25 Å/ps to 10 Å/ps using a pyramid shaped diamond indenter. All MD simulations in this investigation were conducted with an inner cut-off radius 2.0 Å in the AIREBO potential using 0.5 fs time step. Fig. 17(a) shows the resistance force vs the indentation depth for all cases

191

Computational Materials Science 153 (2018) 183–199

A.R. Alian, S.A. Meguid

150 R1=1.70 Å R1=1.80 Å R1=1.90 Å R1=1.96 Å R1=1.97 Å R1=2.00 Å

Tensile stress (GPa)

125

100 75 50

(a) R1= 1.7 Å

25 0

0

0.2

0.4 Tensile strain

0.6

0.8

145

Fracture stress (GPa)

135 125 115

(b) R1= 1.92 Å

105 95 85 75 1.7

1.75 1.8 1.85 1.9 1.95 Cut-off minimum radius R1 (Å)

2

Fig. 10. Tensile MD simulation of a 100 Å × 100 Å graphene sheet using different R1 values in the cut-off function in AIREBO potential: (a) Stress-strain curves and (b) Variation of the fracture stress with R1.

(c) R1= 2.0 Å considered. It is very clear that the sheet resistance varies significantly with the speed/strain rate. Fig. 17(b) shows the variation of the indentation force measured at a 30 Å depth with the indentation speed. The indentation force almost doubled (from 152.2 to 248.2 nN) when the speed increased from 0.25 to 10 A/ps. The relationship between the indentation force and indenter speed can be described with a linear function (Fig. 17(b)).

Fig. 11. A 100 Å × 100 Å Graphene sheet at the onset of fracture due to uniaxial tensile loading along the armchair direction at different values of R1 in the cutoff function: (a) R1 = 1.7 Å, (b) R1 = 1.92 Å, and (c) R1 = 2.0 Å.

simulations in this investigation were conducted with cut-off function parameters R = 2.0 Å and D = 0 Å in the Tersoff potential using 0.5 fs time step. Fig. 18(a) shows the force-depth curves for all speed considered. It is very clear that the BN sheet response is totally altered by the change in indentation speed. The fracture force and fracture depth increased from 12.11 nN and 13.4% to 23.8 nN and 21.8%, respectively, when the indenter speed increased from 0.001 Å/ps to 1 Å/ps.

3.4.3. Effect of indentation speed on BN performance We have also investigated the effect of indenter speed on the BN sheet resistance to indentation by conducting nanoindentation simulations of 100 Å × 100 Å BN sheet with indentation speeds ranging from 0.001 Å/ps to 1 Å/ps using a pyramid shaped diamond indenter. All MD

192

Computational Materials Science 153 (2018) 183–199

A.R. Alian, S.A. Meguid

D= 1.90 Å

R= 1.90 Å R= 2.0 Å R= 2.05 Å R= 2.10 Å R= 2.30 Å

70

D= 2.05 Å 60

D= 2.15 Å

40

Tensile stress (GPa)

Indentation Force (nN)

50

D= 2.20 Å D= 2.30 Å

30 20 10

50 40 30 20 10 0

0 0

10

20 30 Indentation time (ps)

40

0

50

0.025

0.05

0.1

0.125

Tensile strain

(a)

(a) 72.5

55 FMax = 52.7 nN

70

50

Fracture strength (GPa)

Maximum Indentation Force (nN)

0.075

45

40

35

67.5 65 62.5 60 57.5 55

FMin = 31.2 nN

52.5

30 1.9

1.95

2 2.05 2.1 2.15 2.2 Cut-off minimum radius R(Å)

2.25

1.9

2.3

1.95

2

2.05

2.1

2.15

2.2

2.25

2.3

Cut-off distance (Å)

(b)

(b)

Fig. 13. Tensile MD simulations along the armchair direction of a 100 Å × 100 Å BN sheet using different R values in the cut-off function in Tersoff potential: (a) Stress-strain curves and (b) Variation of the fracture stress with R.

Fig. 12. Indentation MD simulations of 100 Å × 100 Å graphene sheet with pyramid-shaped diamond tip using different R1 values in the cut-off function in AIREBO potential: (a) Stress-strain curves and (b) Variation of the maximum indentation force with R1.

r (τ + δτ ) = 2r (τ )−r (τ )−r (τ −δτ ) + a (τ ) δτ 2 + O (δt 4 ) The relationship between the indentation force and indenter speed can be described with a logarithmic function, see Fig. 18(b).

(5)

where δτ is the time step, a (τ ) is the acceleration, and O (δt 4 ) is the truncation error from the Taylor series expansion [56]. It is very clear that the accuracy of the above integration increases significantly with the decrease in the time step. Using a very small time step will result in excessive computations [76]. On the other hand, large time step will result in an unstable molecular system of overlapping atoms and consequently a rapid increase in the potential energy of the system due to the repulsive force [76], as shown in Fig. 19. In this investigation, uniaxial tensile simulations were performed with time steps ranging from 0.1 fs to 15 fs to determine the effect of time step on the predicted mechanical properties of CNTs, graphene, and BN sheet. The output of

3.5. Effect of time step on mechanical properties In MD simulations of large atomic systems and structures, it is definitely better to use the largest possible step size in order to increase the trajectory length and the entire simulation time with the available computational resources [19,76]. The general form of the Verlet algorithm used in MD to solve Newton's equations of motion of the interacting atoms is as follows:

193

Computational Materials Science 153 (2018) 183–199

7.4

7.5

7.35

7.45

Potential energy/atom (ev/atom)

7.3 7.25

7.2 7.15 7.1 7.05 7

6.95

7.4 7.35 7.3 7.25 7.2

7.15

6.9 25

100 175 250 325 400 graphene sheet side size (Å)

7.1

475

25

100

(a)

80

785

75

765 Young's Modulus 745

Fracture strength

70

Young's modulus (GPa)

85 805

Fracture strength (GPa)

Young's Modulus (GPa)

90

825

475

(a)

95

845

175 250 325 400 graphene sheet side size (Å)

700

85

650

80 75

600

70

Young's modulus

550

65

Fracture strength 500

60

450

55

65

400

Fracture strength (GPa)

Potential energy/atom (ev/atom)

A.R. Alian, S.A. Meguid

50

60

725 0

50

100 150 200 250 300 Graphene sheet side size (Å)

350

45 0

(b)

100 200 300 400 Graphene sheet side size (Å)

500

(b)

Fig. 14. Tensile MD simulations of the effect of sheet size on the mechanical properties of graphene (armchair direction): (a) variation of potential energy/ atom with the sheet length and (b) variation of Young’s modulus and the fracture stress with the sheet length.

Fig. 15. Tensile MD simulation concerned with the effect of sheet size on the mechanical properties of BN sheet along the armchair direction: (a) variation of potential energy/atom with the sheet length and (b) variation of Young’s modulus and the fracture stress with the sheet length.

this study will help finding the maximum allowable time step beyond it the system becomes unstable.

Table 4 Mechanical properties of a BN sheet obtained from the literature and the present work. The thickness of the BN sheet was assumed to be 3.33 Å which is equivalent to the interlayer distance between stacking BN layers [74].

3.5.1. Effect of time step on mechanical properties of SWCNTs Uniaxial tension simulations of an armchair SWCNT of (6,6) chirality and length of 66 Å were conducted with an inner cut-off radius 2.0 Å in the AIREBO potential using a time step ranging from 0.1 to 6 fs. One end of the CNT was displaced with a constant speed of 0.01 A/ps Å, which will result in a strain rate of ∼0.16 1/ns and all simulations were conducted with the NVT ensemble at room temperature. Fig. 20 shows that the stress-strain curves of the CNT are significantly influenced by the time step used in the integration process. The obtained elastic moduli and fracture stress and strains for all cases are summarized in Table 6. All properties increased with increasing the time step value

194

Source

Young’s modulus (GPa)

Fracture stress (GPa)

Fracture strain (%)

25 Å × 25 Å BN sheet 500 Å × 500 Å BN sheet MD simulations [15] Ab initio [49] Experimental [50] Experimental [51]

396.0 660.0 692.7 835.0 500–600 811.0

82.4 51.4 114.1 – – –

19.5 8.6 27.0 – – –

Computational Materials Science 153 (2018) 183–199

A.R. Alian, S.A. Meguid

250

v=0.25 A/ps v=0.75 A/ps v=2.5 A/ps v=5.0 A/ps v=10.0 A/ps

Indentation force (nN)

200

150

100

50

0 0

5

10 15 20 Indentation time (ps)

25

30

(a)

(a)

Max. indentation force (nN)

250 Fmax = 11.282 v (A/ps) + 137.02 nN

200 150 100

50 0 0.25

2.5 Indentation speed (A/ps)

(b)

(b)

Fig. 17. Nanoindentation MD simulations of 100 Å × 100 Å graphene sheet. (a) Force-depth curves at different indentation speeds and (b) variation of indentation force at 30 Å depth with the indentation speed.

Fig. 16. Tensile MD simulations of (6,6) armchair CNT with different strain rates. (a) Stress-strain curves for different strain rates and (b) variation of fracture stress with the strain rate.

from 0.1 fs to 1 fs. Specifically, Young’s modulus increased by 14.4% from 804.7 GPa for τ = 0.10 fs to 920.2 GPa for τ = 1.0 fs and the fracture stress and the fracture strain increased by 12.9% and 7.3%, respectively. The mechanical properties then started to fluctuate with increasing time step over 1 fs. Finally, at time steps higher than 5 fs, the system becomes unstable and the CNT did not respond to the applied load. These results conclude that the maximum allowable time step for CNTs MD simulations is around 5 fs.

Table 5 The fracture properties of 100 Å × 100 Å graphene obtained with uniaxial tensile MD simulation using different interatomic potentials. Loading rate

3.5.2. Effect of time step on mechanical properties of graphene Uniaxial tension simulations of a 50 Å × 50 Å graphene sheet along the armchair direction were conducted with an inner cut-off radius

195

Strain rate (ns−1)

Young’s modulus (GPa)

Fracture stress (GPa)

Fracture strain (%)

Quasi-static

0.0002

847.3

94.4

16.5

Dynamic

0.016 0.16 1.6 16.0

871.8 860.9 855.25 844.68

95.8 97.9 99.6 101.7

17.3 18.1 18.7 19.4

Computational Materials Science 153 (2018) 183–199

A.R. Alian, S.A. Meguid

Nanoindentation force (nN)

25

V = 1 Å/ps V = 0.1 Å/ps V = 0.01 Å/ps V = 0.001 Å/ps

20

15

10

5

0

0

5

10 15 Indentation depth (A)

20

(a)

Fig. 20. Stress-strain curves from the conducted tensile MD simulations of (6,6) armchair CNT with different time steps.

Max. indentation force (nN)

30 2.0 Å in the AIREBO potential using a time step ranging from 0.1 to 10 fs. One end of the CNT was displaced with a constant speed of 0.01 A/ps which will result in a strain rate of ∼0.12 1/ns and all simulations were conducted with the NVT ensemble at room temperature. Fig. 21 shows that the stress-strain curves of the graphene sheet are significantly influenced by the time step. We summarized the determined elastic moduli and fracture stress and strains for all time steps in Table 7. The results show that the fracture stress and strain at time steps above than 1 fs are larger than those obtained for lower time steps. However, Young’s modulus is generally fluctuating over the entire range around 710 GPa. In addition, Young’s modulus value started to drop at time steps higher than 7 fs as the system did not responded to the applied deformation and the MD simulations become unstable due to the large truncation errors. This investigation concludes that elastic MD simulations of graphene can be conducted with time steps as high as 7 fs, while fracture simulations can be conducted with time steps up to 1 fs.

Fmax = 1.7505 ln (v) + 23.785 nN 20

10

0 0.001

0.01 0.1 Indentation speed (Å/ps)

1

3.5.3. Effect of time step on mechanical properties of BN sheets Uniaxial tension simulations of a 50 Å × 50 Å BN sheet along the armchair direction were conducted with an inner cut-off radius 2.0 Å in the Tersoff potential using a time step ranging from 0.1 to 20 fs. One end of the CNT was displaced with a constant speed of 0.01 A/ps Å which will result in a strain rate of ∼0.12 1/ns and all simulations were

(b) Fig. 18. Nanoindentation MD simulation of a 100 Å × 100 Å BN sheet. (a) Indentation force-depth relationship at different indentation speeds and (b) Variation of indentation force with the indentation speed at 30 Å depth.

Fig. 19. Trajectories evolution of atoms with time: (a) initial locations and velocities at time t, (b) locations of atoms when using appropriate time step (δτ ≅ 0.033–0.01 bond vibration period), and (c) overlapping of atoms when using large time step, known as molecular structure explosion. 196

Computational Materials Science 153 (2018) 183–199

A.R. Alian, S.A. Meguid

Table 6 Elastic and fracture properties of an armchair SWCNT of (6,6) chirality and length of 66 Å obtained with uniaxial tensile MD simulation using different time steps. Time step (fs)

Young's modulus (GPa)

Fracture stress (GPa)

Fracture strain (%)

0.1 0.25 0.5 0.75 1 2 3 4 5

804.7 835.4 861.2 893.6 920.2 843 844.4 846.7 957.7

91.7 95.7 97.5 101.2 103.5 100.3 100.6 101.2 106.5

17.7 17.9 17.8 18.8 19 18.8 18.9 19.2 18.5

Fig. 22. Stress-strain curves from the conducted tensile MD simulations of 50 Å × 50 Å BN sheet along the armchair direction with different time steps.

Table 8 The elastic and fracture properties of a 50 Å × 50 Å BN sheet along the armchair direction with different time steps.

Fig. 21. The Stress-strain curves from the conducted tensile MD simulations of a 50 Å × 50 Å graphene sheet along the armchair direction with different time steps.

Fracture stress (GPa)

Fracture strain (%)

0.1 0.25 0.5 0.75 1 2 4 5 7 8 10

678.3 755.85 606.64 715.82 749.98 722.1 666.76 726.71 723.23 622.66 448

80.4 91.6 85.1 83.3 88.7 104.9 106.2 99.7 101.4 101.6 100.6

13.3 15.5 14.7 13.6 14.7 19.2 19.1 20 17.1 17.6 17.1

Fracture stress (GPa)

Fracture strain (%)

0.1 0.25 0.5 0.75 1 2 3 4 5 10 15

568.9 552.96 509.34 539 477.31 555.9 598 562.2 579.7 613.9 579.5

55.9 61.4 60.6 68.7 73.5 75.2 73.7 86 83.5 86.5 80.9

10.4 11.7 12.2 13.6 15.4 15.3 14.4 18.7 17.4 18.8 17.5

4. Conclusions Currently, there exists anomalies and inconsistencies in the literature among different atomistic simulations. The purpose of this extensive study is to shed light on the sources of these anomalies and inconsistencies and help guide future atomistic research. Specifically, we conducted comprehensive molecular dynamics studies to investigate the effect of the type of the interatomic potential (AIREBO, Tersoff, CVFF, and ReaxFF), cut-off function parameters in AIREBO and Tersoff potentials, domain size (sheet widths ranging from 25 Å up to 500 Å), strain rate (ε ̇ ranging from 2 × 10−4to ∼ 16 ns−1) , indentation speed (v ranging from 0.001 to 10Å/ps, and time step (τ ranging from 0.1 to 20 fs) on the elastic and fracture properties of CNTs, graphene and boron nitride monolayers. Both uniaxial tensile and nanoindentation simulations were performed to explore the influence of all these parameters on the deformation and failure mechanisms in each material under different loading conditions. The following is a summary of our findings:

Table 7 The elastic and fracture properties of a 50 Å × 50 Å graphene sheet along the armchair direction with different time steps. Young's modulus (GPa)

Young's modulus (GPa)

than 10 fs. We summarized the obtained elastic moduli and fracture stress and strains for all time steps in Table 8. This investigation concludes that elastic MD simulations of BN sheet can be conducted with time steps as high as 10 fs, while fracture simulations can be conducted with time steps up to 1 fs.

conducted with the NVT ensemble at room temperature. Fig. 22 shows that the stress-strain curves of the graphene sheet are significantly influenced by the time step with Young’s modulus, fracture stress, and strain ranging from 477 to 614 GPa, 55.9 to 86.5 GPa, and 10.8% to 18.8%, respectively. In detailed, Young’s modulus value fluctuates with the time step around an average value of 555 GPa. The fracture stress and strain increase with the time step used in the simulation. The results also show that the BN sheet become unstable for time steps higher

Time step (fs)

Time step (fs)

197

Computational Materials Science 153 (2018) 183–199

A.R. Alian, S.A. Meguid

• The elastic moduli obtained by the Tersoff potential for CNT and •

•

•

•

• •

References

graphene are 20% higher than those obtained by AIREBO potential. CVFF and ReaxFF force fields can successfully predict the elastic behavior of CNT and graphene for tensile strains up to ∼5% and 9%, respectively. The harmonic representation of the covalent bonds between carbon atoms in CVFF force field leads to a linear stress-strain behavior without any possibility of bond breakage. ReaxFF force field overestimates the fracture stress of CNT by 50% due to the noticeable strain hardening at tensile strains over 10%. However, ReaxFF predicts lower fracture stress and fracture strain of graphene. The cut-off function parameters in AIREBO and Tersoff potentials only influence the fracture behavior of the simulated structures. A cut-off radius R1 in AIREBO potentials lower than 1.97 Å leads to strain hardening behavior in CNT and graphene structures. Cut-off inner and outer radii of 2.0 Å in Tersoff potential offer the best match with reported DFT studies. Young’s modulus and the potential energy/atom of graphene and BN sheets increase with increasing the sheet size due to the reduction of edge atoms. The fracture stress and fracture strain of graphene and BN sheets decrease with increasing the sheet size and all mechanical properties reach the bulk level above a sheet length of 200 Å. The elastic properties of the considered nanostructures are independent of strain rates below 2 × 10−9 s1. On the other hand, the fracture stress and fracture strain increase logarithmically with the strain rate. The atomic system becomes unstable during tensile MD simulations when deformed at strain rates higher than 2 × 10−9 s−1. The maximum indentation force and fracture depth of graphene and BN sheets decrease with the decrease in the indenter speed. The elastic and fracture properties of SWCNTs increase linearly with increasing the time step within the 0.1–1.0 fs range. The mechanical properties fluctuate with the increase in the time step before and the system become unstable at 5 fs. The elastic behavior of graphene and BN sheets can be successfully captured with the aid of tensile MD simulations conducted using time steps as high as 7 fs and 10 fs, respectively, while the fracture properties show a significant increase for time steps over 1 fs.

[1] A.R. Alian, S.A. Meguid, Multiscale modeling of nanoreinforced composites, in: S.A. Meguid (Ed.), Adv. Nanocomposites, Springer International Publishing, Cham, 2016, pp. 1–39, , http://dx.doi.org/10.1007/978-3-319-31662-8_1. [2] W.F. van Gunsteren, H.J.C. Berendsen, Computer simulation of molecular dynamics: methodology, applications, and perspectives in chemistry, Angew. Chem. Int. Ed. Engl. 29 (1990) 992–1023, http://dx.doi.org/10.1002/anie.199009921. [3] B. Javvaji, P.R. Budarapu, V.K. Sutrakar, D.R. Mahapatra, M. Paggi, G. Zi, T. Rabczuk, Mechanical properties of graphene: molecular dynamics simulations correlated to continuum based scaling laws, Comput. Mater. Sci. 125 (2016) 319–327, http://dx.doi.org/10.1016/j.commatsci.2016.08.016. [4] J. Zhang, F. Xu, Y. Hong, Q. Xiong, J. Pan, A comprehensive review on the molecular dynamics simulation of the novel thermal properties of graphene, RSC Adv. 5 (2015) 89415–89426, http://dx.doi.org/10.1039/C5RA18579C. [5] R. Rafiee, R.M. Moghadam, On the modeling of carbon nanotubes: A critical review, Compos. Part B Eng. 56 (2014) 435–449, http://dx.doi.org/10.1016/j.compositesb. 2013.08.037. [6] A. Verma, A. Parashar, M. Packirisamy, Atomistic modeling of graphene/hexagonal boron nitride polymer nanocomposites: a review: atomistic modeling of graphene/ h-BN polymer nanocomposites, Wiley Interdiscip. Rev. Comput. Mol. Sci. (2017) e1346, http://dx.doi.org/10.1002/wcms.1346. [7] A.R. Alian, S.A. Meguid, Large-scale atomistic simulations of CNT-reinforced thermoplastic polymers, Compos. Struct. 191 (2018) 221–230, http://dx.doi.org/10. 1016/j.compstruct.2018.02.056. [8] J. Zhang, X. He, L. Yang, G. Wu, J. Sha, C. Hou, C. Yin, A. Pan, Z. Li, Y. Liu, Effect of tensile strain on thermal conductivity in monolayer graphene nanoribbons: a molecular dynamics study, Sensors 13 (2013) 9388–9395, http://dx.doi.org/10.3390/ s130709388. [9] A.R. Alian, S.I. Kundalwal, S.A. Meguid, Multiscale modeling of carbon nanotube epoxy composites, Polymer 70 (2015) 149–160, http://dx.doi.org/10.1016/j. polymer.2015.06.004. [10] A.R. Alian, S.I. Kundalwal, S.A. Meguid, Interfacial and mechanical properties of epoxy nanocomposites using different multiscale modeling schemes, Compos. Struct. 131 (2015) 545–555, http://dx.doi.org/10.1016/j.compstruct.2015.06.014. [11] Y. Li, M. Kröger, A theoretical evaluation of the effects of carbon nanotube entanglement and bundling on the structural and mechanical properties of buckypaper, Carbon 50 (2012) 1793–1806, http://dx.doi.org/10.1016/j.carbon.2011. 12.027. [12] S.S. Han, J.K. Kang, H.M. Lee, A.C.T. van Duin, W.A. Goddard, The theoretical study on interaction of hydrogen with single-walled boron nitride nanotubes. I. The reactive force field ReaxFFHBN development, J. Chem. Phys. 123 (2005) 114703, http://dx.doi.org/10.1063/1.1999628. [13] B.D. Jensen, K.E. Wise, G.M. Odegard, Simulation of the elastic and ultimate tensile properties of diamond, graphene, carbon nanotubes, and amorphous carbon using a revised ReaxFF parametrization, J. Phys. Chem. A 119 (2015) 9710–9721, http:// dx.doi.org/10.1021/acs.jpca.5b05889. [14] W.H. Duan, Q. Wang, K.M. Liew, X.Q. He, Molecular mechanics modeling of carbon nanotube fracture, Carbon 45 (2007) 1769–1776, http://dx.doi.org/10.1016/j. carbon.2007.05.009. [15] S. Zhao, J. Xue, Mechanical properties of hybrid graphene and hexagonal boron nitride sheets as revealed by molecular dynamic simulations, J. Phys. Appl. Phys. 46 (2013) 135303, http://dx.doi.org/10.1088/0022-3727/46/13/135303. [16] M.-Q. Le, Y. Umeno, Fracture of monolayer boronitrene and its interface with graphene, Int. J. Fract. 205 (2017) 151–168, http://dx.doi.org/10.1007/s10704017-0188-0. [17] N. Dewapriya, Molecular dynamics study of effects of geometric defects on the mechanical properties of graphene, Univ. British Columbia (2012), http://dx.doi. org/10.14288/1.0072708. [18] A.R. Alian, S.A. Meguid, Hybrid molecular dynamics–finite element simulations of the elastic behavior of polycrystalline graphene, Int. J. Mech. Mater. Des. (2017), http://dx.doi.org/10.1007/s10999-017-9389-y. [19] B.D. Jensen, K.E. Wise, G.M. Odegard, The effect of time step, thermostat, and strain rate on ReaxFF simulations of mechanical failure in diamond, graphene, and carbon nanotube, J. Comput. Chem. 36 (2015) 1587–1596, http://dx.doi.org/10.1002/jcc. 23970. [20] K. Mylvaganam, L. Zhang, Important issues in a molecular dynamics simulation for characterising the mechanical properties of carbon nanotubes, Carbon 42 (2004) 2025–2032, http://dx.doi.org/10.1016/j.carbon.2004.04.004. [21] Z.-Y. Ong, E. Pop, Molecular dynamics simulation of thermal boundary conductance between carbon nanotubes and SiO2, Phys. Rev. B 81 (2010), http://dx.doi.org/10. 1103/PhysRevB.81.155408. [22] H. Bu, Y. Chen, M. Zou, H. Yi, K. Bi, Z. Ni, Atomistic simulations of mechanical properties of graphene nanoribbons, Phys. Lett. A 373 (2009) 3359–3362, http:// dx.doi.org/10.1016/j.physleta.2009.07.048. [23] P.R. Budarapu, B. Javvaji, V.K. Sutrakar, D. Roy Mahapatra, G. Zi, T. Rabczuk, Crack propagation in graphene, J. Appl. Phys. 118 (2015) 064307, http://dx.doi. org/10.1063/1.4928316. [24] M.Q. Chen, S.S. Quek, Z.D. Sha, C.H. Chiu, Q.X. Pei, Y.W. Zhang, Effects of grain size, temperature and strain rate on the mechanical properties of polycrystalline graphene – A molecular dynamics study, Carbon 85 (2015) 135–146, http://dx.doi. org/10.1016/j.carbon.2014.12.092. [25] L. Yi, Z. Yin, Y. Zhang, T. Chang, A theoretical evaluation of the temperature and strain-rate dependent fracture strength of tilt grain boundaries in graphene, Carbon 51 (2013) 373–380, http://dx.doi.org/10.1016/j.carbon.2012.08.069.

Finally, we believe that the findings of this work would help in understanding the role played by certain atomistic parameters which ultimately lead to inconsistent and contradictory predictions among the reported data in the literature and also in finding the appropriate parameters for each specific loading and material type in future studies. However, future investigations, such as the one reported here, should be considered to examine the effect of the parameters which govern MD simulations of the mechanical properties of larger systems such as advanced polymeric and metallic based nanocomposites.

Acknowledgements The authors wish to thank the Natural Sciences and Engineering Research Council of Canada (NSERC) for their financial support of this research (Project No. RGPIN-2018-03804). The first author wishes to thank Cairo University for giving him leave of absence to complete this work at the University of Toronto.

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.commatsci.2018.06. 028.

198

Computational Materials Science 153 (2018) 183–199

A.R. Alian, S.A. Meguid

[52] J.M. Haile, Molecular Dynamics Simulation: Elementary Methods, Wiley, New York, 1992. [53] S. Plimpton, Fast parallel algorithms for short-range molecular dynamics, J. Comput. Phys. 117 (1995) 1–19, http://dx.doi.org/10.1006/jcph.1995.1039. [54] C. Sevik, A. Kinaci, J.B. Haskins, T. Çağın, Characterization of thermal transport in low-dimensional boron nitride nanostructures, Phys. Rev. B 84 (2011), http://dx. doi.org/10.1103/PhysRevB.84.085409. [55] A.T. Hagler, E. Huler, S. Lifson, Energy functions for peptides and proteins. I. Derivation of a consistent force field including the hydrogen bond from amide crystals, J. Am. Chem. Soc. 96 (1974) 5319–5327, http://dx.doi.org/10.1021/ ja00824a004. [56] L. Verlet, Computer “Experiments” on classical fluids. I. Thermodynamical properties of Lennard-Jones molecules, Phys. Rev. 159 (1967) 98–103, http://dx.doi. org/10.1103/PhysRev. 159.98. [57] A.R. Alian, S. El-Borgi, S.A. Meguid, Multiscale modeling of the effect of waviness and agglomeration of CNTs on the elastic properties of nanocomposites, Comput. Mater. Sci. 117 (2016) 195–204, http://dx.doi.org/10.1016/j.commatsci.2016.01. 029. [58] A.R. Alian, S.A. Meguid, Molecular dynamics simulations of the effect of waviness and agglomeration of CNTs on interface strength of thermoset nanocomposites, Phys. Chem. Chem. Phys. 19 (2017) 4426–4434, http://dx.doi.org/10.1039/ C6CP07464B. [59] A.R. Alian, M.A.N. Dewapriya, S.A. Meguid, Molecular dynamics study of the reinforcement effect of graphene in multilayered polymer nanocomposites, Mater. Des. 124 (2017) 47–57, http://dx.doi.org/10.1016/j.matdes.2017.03.052. [60] A.R. Alian, S.A. Meguid, S.I. Kundalwal, Unraveling the influence of grain boundaries on the mechanical properties of polycrystalline carbon nanotubes, Carbon 125 (2017) 180–188, http://dx.doi.org/10.1016/j.carbon.2017.09.056. [61] R. LeSar, Introduction to computational materials science: fundamentals to applications, Cambridge University Press, Cambridge , New York, 2013 (ISBN: 9780521845878). [62] M.A. González, Force fields and molecular dynamics simulations, Éc. Thématique Société Fr. Neutron. 12 (2011) 169–200, http://dx.doi.org/10.1051/sfn/ 201112009. [63] D. Frenkel, B. Smit, Understanding Molecular Simulation: From Algorithms to Applications, second ed., Academic Press, San Diego, 2002 ISBN: 9780080519982. [64] D.W. Brenner, Empirical potential for hydrocarbons for use in simulating the chemical vapor deposition of diamond films, Phys. Rev. B 42 (1990) 9458–9471, http://dx.doi.org/10.1103/PhysRevB.42.9458. [65] S.A. Meguid, A.R. Alian, M.A.N. Dewapriya, Atomistic modelling of nanoindentation of multilayered graphene-reinforced nanocomposites, in: S.A. Meguid, G.J. Weng (Eds.), Micromechanics Nanomechanics Compos. Solids, Springer International Publishing, Cham, 2018, pp. 39–70, , http://dx.doi.org/10.1007/9783-319-52794-9_2. [66] S. Thamaraikannan, S.C. Pradhan, atomistic study of carbon nanotubes: effect of cut-off distance, in: The Minerals, Metals and Materials Society (Eds.), TMS 2016 145th Annu. Meet. Exhib. Springer International Publishing, Cham, 2016, pp. 293–300, , http://dx.doi.org/10.1007/978-3-319-48254-5_36. [67] M. Izadifar, R. Abadi, A.N. Jam, T. Rabczuk, Investigation into the effect of doping of boron and nitrogen atoms in the mechanical properties of single-layer polycrystalline graphene, Comput. Mater. Sci. 138 (2017) 435–447, http://dx.doi.org/ 10.1016/j.commatsci.2017.06.038. [68] V. Verma, V.K. Jindal, K. Dharamvir, Elastic moduli of a boron nitride nanotube, Nanotechnology 18 (2007) 435711, http://dx.doi.org/10.1088/0957-4484/18/43/ 435711. [69] Y. Wei, J. Wu, H. Yin, X. Shi, R. Yang, M. Dresselhaus, The nature of strength enhancement and weakening by pentagon–heptagon defects in graphene, Nat. Mater. 11 (2012) 759–763, http://dx.doi.org/10.1038/nmat3370. [70] B. Zhang, L. Mei, H. Xiao, Nanofracture in graphene under complex mechanical stresses, Appl. Phys. Lett. 101 (2012) 121915, http://dx.doi.org/10.1063/1. 4754115. [71] T. Zhang, X. Li, S. Kadkhodaei, H. Gao, Flaw insensitive fracture in nanocrystalline graphene, Nano Lett. 12 (2012) 4605–4610, http://dx.doi.org/10.1021/ nl301908b. [72] M.A.N. Dewapriya, R.K.N.D. Rajapakse, A.S. Phani, Atomistic and continuum modelling of temperature-dependent fracture of graphene, Int. J. Fract. 187 (2014) 199–212, http://dx.doi.org/10.1007/s10704-014-9931-y. [73] K.G.S. Dilrukshi, M.A.N. Dewapriya, U.G.A. Puswewala, Size dependency and potential field influence on deriving mechanical properties of carbon nanotubes using molecular dynamics, Theor. Appl. Mech. Lett. 5 (2015) 167–172, http://dx.doi.org/ 10.1016/j.taml.2015.05.005. [74] E.-S. Oh, Elastic properties of boron-nitride nanotubes through the continuum lattice approach, Mater. Lett. 64 (2010) 859–862, http://dx.doi.org/10.1016/j.matlet. 2010.01.041. [75] A.J. Kulkarni, M. Zhou, F.J. Ke, Orientation and size dependence of the elastic properties of zinc oxide nanobelts, Nanotechnology 16 (2005) 2749–2756, http:// dx.doi.org/10.1088/0957-4484/16/12/001. [76] S. Kim, Issues on the choice of a proper time step in molecular dynamics, Phys. Procedia. 53 (2014) 60–62, http://dx.doi.org/10.1016/j.phpro.2014.06.027.

[26] X. Liu, Q.-S. Yang, Plastic deformation and failure mechanisms of collapsed-carbon nanotube fibers by coarse-grained molecular dynamic simulations, Int. J. Plast. 64 (2015) 104–112, http://dx.doi.org/10.1016/j.ijplas.2014.08.004. [27] H. Li, H. Zhang, X. Cheng, The effect of temperature, defect and strain rate on the mechanical property of multi-layer graphene: coarse-grained molecular dynamics study, Phys. E Low-Dimens. Syst. Nanostruct. 85 (2017) 97–102, http://dx.doi.org/ 10.1016/j.physe.2016.07.003. [28] C. Wei, K. Cho, D. Srivastava, Tensile strength of carbon nanotubes under realistic temperature and strain rate, Phys. Rev. B 67 (2003), http://dx.doi.org/10.1103/ PhysRevB.67.115407. [29] Y.-Y. Zhang, Q.-X. Pei, Y.-W. Mai, Y.-T. Gu, Temperature and strain-rate dependent fracture strength of graphynes, J. Phys. Appl. Phys. 47 (2014) 425301, http://dx. doi.org/10.1088/0022-3727/47/42/425301. [30] H. Zhao, N.R. Aluru, Temperature and strain-rate dependent fracture strength of graphene, J. Appl. Phys. 108 (2010) 064321, http://dx.doi.org/10.1063/1. 3488620. [31] T. Han, Y. Luo, C. Wang, Effects of temperature and strain rate on the mechanical properties of hexagonal boron nitride nanosheets, J. Phys. Appl. Phys. 47 (2014) 025303, http://dx.doi.org/10.1088/0022-3727/47/2/025303. [32] B. Mortazavi, Y. Rémond, Investigation of tensile response and thermal conductivity of boron-nitride nanosheets using molecular dynamics simulations, Phys. E LowDimens. Syst. Nanostruct. 44 (2012) 1846–1852, http://dx.doi.org/10.1016/j. physe.2012.05.007. [33] T. Zhang, X. Li, H. Gao, Fracture of graphene: a review, Int. J. Fract. 196 (2015) 1–31, http://dx.doi.org/10.1007/s10704-015-0039-9. [34] A. Agius Anastasi, K. Ritos, G. Cassar, M.K. Borg, Mechanical properties of pristine and nanoporous graphene, Mol. Simul. 42 (2016) 1502–1511, http://dx.doi.org/ 10.1080/08927022.2016.1209753. [35] Z. Ni, H. Bu, M. Zou, H. Yi, K. Bi, Y. Chen, Anisotropic mechanical properties of graphene sheets from molecular dynamics, Phys. B Condens. Matter. 405 (2010) 1301–1306, http://dx.doi.org/10.1016/j.physb.2009.11.071. [36] Y.J. Sun, F. Ma, K.W. Xu, Size dependent mechanical properties of graphene nanoribbons: molecular dynamics simulation, Mater. Sci. Forum. 749 (2013) 456–460, http://dx.doi.org/10.4028/www.scientific.net/MSF.749.456. [37] S. Thomas, K.M. Ajith, M.C. Valsakumar, Directional anisotropy, finite size effect and elastic properties of hexagonal boron nitride, J. Phys. Condens. Matter. 28 (2016) 295302, http://dx.doi.org/10.1088/0953-8984/28/29/295302. [38] Q. Lu, W. Gao, R. Huang, Atomistic simulation and continuum modeling of graphene nanoribbons under uniaxial tension, Model. Simul. Mater. Sci. Eng. 19 (2011) 054006, http://dx.doi.org/10.1088/0965-0393/19/5/054006. [39] H. Zhao, K. Min, N.R. Aluru, Size and chirality dependent elastic properties of graphene nanoribbons under uniaxial tension, Nano Lett. 9 (2009) 3012–3015, http://dx.doi.org/10.1021/nl901448z. [40] Q. Lu, R. Huang, Excess energy and deformation along free edges of graphene nanoribbons, Phys. Rev. B 81 (2010), http://dx.doi.org/10.1103/PhysRevB.81. 155410. [41] M.-F. Yu, B.S. Files, S. Arepalli, R.S. Ruoff, Tensile loading of ropes of single wall carbon nanotubes and their mechanical properties, Phys. Rev. Lett. 84 (2000) 5552–5555, http://dx.doi.org/10.1103/PhysRevLett. 84.5552. [42] M. Yu, Strength and breaking mechanism of multiwalled carbon nanotubes under tensile load, Science 287 (2000) 637–640, http://dx.doi.org/10.1126/science.287. 5453.637. [43] G. Van Lier, C. Van Alsenoy, V. Van Doren, P. Geerlings, Ab initio study of the elastic properties of single-walled carbon nanotubes and graphene, Chem. Phys. Lett. 326 (2000) 181–185, http://dx.doi.org/10.1016/S0009-2614(00)00764-8. [44] M. Meo, M. Rossi, Tensile failure prediction of single wall carbon nanotube, Eng. Fract. Mech. 73 (2006) 2589–2599, http://dx.doi.org/10.1016/j.engfracmech. 2006.05.005. [45] Q.X. Pei, Y.W. Zhang, V.B. Shenoy, A molecular dynamics study of the mechanical properties of hydrogen functionalized graphene, Carbon 48 (2010) 898–904, http://dx.doi.org/10.1016/j.carbon.2009.11.014. [46] Z.G. Fthenakis, N.N. Lathiotakis, Graphene allotropes under extreme uniaxial strain: an ab initio theoretical study, Phys. Chem. Chem. Phys. 17 (2015) 16418–16427, http://dx.doi.org/10.1039/C5CP02412A. [47] F. Liu, P. Ming, J. Li, Ab initio calculation of ideal strength and phonon instability of graphene under tension, Phys. Rev. B 76 (2007), http://dx.doi.org/10.1103/ PhysRevB.76.064120. [48] C. Lee, X. Wei, J.W. Kysar, J. Hone, Measurement of the elastic properties and intrinsic strength of monolayer graphene, Science 321 (2008) 385–388, http://dx. doi.org/10.1126/science.1157996. [49] Q. Peng, W. Ji, S. De, Mechanical properties of the hexagonal boron nitride monolayer: Ab initio study, Comput. Mater. Sci. 56 (2012) 11–17, http://dx.doi. org/10.1016/j.commatsci.2011.12.029. [50] D. Golberg, P.M.F.J. Costa, O. Lourie, M. Mitome, X. Bai, K. Kurashima, C. Zhi, C. Tang, Y. Bando, Direct force measurements and kinking under elastic deformation of individual multiwalled boron nitride nanotubes, Nano Lett. 7 (2007) 2146–2151, http://dx.doi.org/10.1021/nl070863r. [51] A. Bosak, J. Serrano, M. Krisch, K. Watanabe, T. Taniguchi, H. Kanda, Elasticity of hexagonal boron nitride: inelastic x-ray scattering measurements, Phys. Rev. B 73 (2006), http://dx.doi.org/10.1103/PhysRevB.73.041402.

199