Molecular Dynamics-Based Cohesive Zone Model for Mg / Mg 17 Al 12 Interface

: The fracture of the Mg / Mg 17 Al 12 interface was investigated by molecular dynamics simulations. The interface crack extends in a brittle manner without noticeable plasticity. The distributions of normal stress and separation along the interface were examined to render a quantitative picture of the fracture process. A normal traction–separation curve was extracted from simulation and compared with three cohesive zone models, i.e., cubic polynomial cohesive zone model, exponential cohesive zone model, and bilinear cohesive zone model. The exponential cohesive zone model exhibits the best agreement with simulation results, followed by the bilinear cohesive zone model.


Introduction
As one of the lightest structural materials, magnesium (Mg) alloys have attracted substantial interest from both the scientific community and various industries such as transportation and aerospace, due to their high specific strength, excellent castability, and superior damping capacity [1][2][3]. However, compared to structural materials such as steel, titanium and aluminum alloys, their fracture toughness is lower, which restricts their applications. It has been reported that the matrix/precipitate interface in Mg alloys is the potential site for crack nucleation [4,5]. Therefore, it is of enormous significance to investigate the fracture properties of the matrix/precipitate interface.
Cohesive zone models (CZMs), which were first proposed by Dugdale [6] and Barenblatt [7], have been widely employed to study interfacial fracture. In the framework of CZM, fracture is described by a traction-separation relation, which is typically obtained based on macroscale fracture properties such as fracture strength and toughness. Such traction-separation relations represent the collective response of all microstructure constituents within the specimen, rather than the unique response of the interface where fracture takes place [8]. To establish a CZM that can capture the local response in the vicinity of a crack, one needs to have the knowledge of local nanoscale properties, which can be obtained through molecular dynamics (MD) simulations. After the pioneering work of Gall et al. [9] that parameterized CZMs for Al/Si bimaterial by using MD simulation results as input, substantial efforts have been devoted to developing MD-based CZMs for various materials, including Al grain boundary [8], SiC/Mg nanocomposite [10], Fe single crystal [11], Al/Si bimaterials [12], and FCC single crystals [13]. However, despite numerous efforts, no attention has been paid to the fracture of the matrix/precipitate interface in Mg alloys.
In this work, we investigated by MD simulations the fracture of the Mg/Mg 17 Al 12 interface, which is commonly found in commercial Mg-Al alloys. Tensile loading was applied to the bimaterial with an interface crack to initiate crack propagation. The evolution of atomic configuration during tensile loading was visualized to reveal the fracture mechanism. The traction-separation relation for the Mg/Mg 17 Al 12 interface was extracted from MD simulations and then compared with three CZMs.

Materials and Methods
As shown in Figure 1, the simulation model has a dimension of 38.7 (w) × 29.6 (h) × 14.6 (t) nm, with the Mg 17 [14]. The precipitate, which has a BCC structure with 58 atoms per unit cell, was constructed according to the structure given by Gharghouri et al. [15]. An interface crack of length c was inserted into the simulation system by blocking the interactions between atoms above and below crack surfaces, followed by a relaxation process by the conjugate gradient method. Then, boundary atoms (shaded by light orange) were frozen, whereas active atoms (all atoms except boundary atoms) were equilibrated in the NVT ensemble at 1 K for 20 ps. Subsequently, a mode I loading with a strain rate of 6.8 × 10 8 s −1 was applied to the bimaterial by moving the two boundaries oppositely along the z-axis. We have checked the possible effects arising from model size and determined that the current model is large enough to yield steady state propagation.
Metals 2020, 10, x FOR PEER REVIEW 2 of 6 an interface crack to initiate crack propagation. The evolution of atomic configuration during tensile loading was visualized to reveal the fracture mechanism. The traction-separation relation for the Mg/Mg17Al12 interface was extracted from MD simulations and then compared with three CZMs.

Materials and Methods
As shown in Figure [14]. The precipitate, which has a BCC structure with 58 atoms per unit cell, was constructed according to the structure given by Gharghouri et al. [15]. An interface crack of length c was inserted into the simulation system by blocking the interactions between atoms above and below crack surfaces, followed by a relaxation process by the conjugate gradient method. Then, boundary atoms (shaded by light orange) were frozen, whereas active atoms (all atoms except boundary atoms) were equilibrated in the NVT ensemble at 1 K for 20 ps. Subsequently, a mode I loading with a strain rate of 6.8 × 10 8 s −1 was applied to the bimaterial by moving the two boundaries oppositely along the z-axis. We have checked the possible effects arising from model size and determined that the current model is large enough to yield steady state propagation. MD simulations were carried out using the large-scale atomic/molecular massively parallel simulator (LAMMPS) [16] with a timestep of 1 fs. Atomeye, which is an atomistic configuration viewer developed by Li [17], was used to visualize atomic configurations. The Mg-Mg, Al-Al, Mg-Al atomic interactions were described by the embedded atom method (EAM) potential [18]. The stress tensor of the simulation system was calculated according to the definition of virial stress. A horizontal layer (shaded by light yellow in Figure 1) consisting of half Mg and half Mg17Al12 was divided into n (45) small regions, each of which is a cohesive zone volume element (CZVE) [8], with a dimension of 0.9 (lx) × 2.1 (lz) × 14.6 (t) nm. Two variables were used to monitor the state of each CZVE, i.e., Tn and un. Tn is the average atomic stress σzz of all atoms in each CZVE, while un is the average atomic displacement Δz of all atoms in the Mg17Al12 half with reference to the Mg half.

Results
The atomic configuration evolution of the bimaterial as the increase of strain is shown Figure 2, with atoms colored according to the zz component of the atomic stress. Due to its hexagonal closepacked (HCP) structure, the available deformation mode in Mg is limited, leading to its relatively low ductility [19]. In this study, tensile loading along the [0001] direction, which is the <c> axis of Mg, was applied to the simulation system. At room temperature, the <c + a> slip system is difficult to activate because of its high critical resolved shear stress, and therefore, twinning is the primary MD simulations were carried out using the large-scale atomic/molecular massively parallel simulator (LAMMPS) [16] with a timestep of 1 fs. Atomeye, which is an atomistic configuration viewer developed by Li [17], was used to visualize atomic configurations. The Mg-Mg, Al-Al, Mg-Al atomic interactions were described by the embedded atom method (EAM) potential [18]. The stress tensor of the simulation system was calculated according to the definition of virial stress. A horizontal layer (shaded by light yellow in Figure 1) consisting of half Mg and half Mg 17 Al 12 was divided into n (45) small regions, each of which is a cohesive zone volume element (CZVE) [8], with a dimension of 0.9 (l x ) × 2.1 (l z ) × 14.6 (t) nm. Two variables were used to monitor the state of each CZVE, i.e., T n and u n . T n is the average atomic stress σ zz of all atoms in each CZVE, while u n is the average atomic displacement ∆z of all atoms in the Mg 17 Al 12 half with reference to the Mg half.

Results
The atomic configuration evolution of the bimaterial as the increase of strain is shown Figure 2, with atoms colored according to the zz component of the atomic stress. Due to its hexagonal close-packed (HCP) structure, the available deformation mode in Mg is limited, leading to its relatively low ductility [19]. In this study, tensile loading along the [0001] direction, which is the <c> axis of Mg, was applied to the simulation system. At room temperature, the <c + a> slip system is difficult to activate because of its high critical resolved shear stress, and therefore, twinning is the primary mechanism to accommodate plastic deformation in the <c> axis [19]. However, neither twinning nor any dislocation activity was observed in our simulations. Instead, the interface crack was observed to propagate in a brittle manner without any noticeable plasticity, despite crack surfaces not being very clean, as shown in Figure 2. This finding is consistent with the previous experimental observation that decohesion occurs in the Mg/precipitate interface [5]. Due to the existence of the inserted interface crack, stress concentration was induced near the crack tips, which may make crack propagation along the interface more favorable than twinning or dislocation activities. Stress concentration at the crack tips can be observed from Figure 2. In addition, simulations were carried out at a low temperature of 1 K, which might suppress plastic deformation to some extent.
Metals 2020, 10, x FOR PEER REVIEW 3 of 6 σzz mechanism to accommodate plastic deformation in the <c> axis [19]. However, neither twinning nor any dislocation activity was observed in our simulations. Instead, the interface crack was observed to propagate in a brittle manner without any noticeable plasticity, despite crack surfaces not being very clean, as shown in Figure 2. This finding is consistent with the previous experimental observation that decohesion occurs in the Mg/precipitate interface [5]. Due to the existence of the inserted interface crack, stress concentration was induced near the crack tips, which may make crack propagation along the interface more favorable than twinning or dislocation activities. Stress concentration at the crack tips can be observed from Figure 2. In addition, simulations were carried out at a low temperature of 1 K, which might suppress plastic deformation to some extent. To render a general quantitative picture of the fracture process, the distribution of Tn, un, and Δz along the interface for two strain values corresponding to Figure 2a, c is shown in Figure 3. Both Tn and un exhibit a near-symmetrical distribution with respect to an axis parallel to the z-axis. un remains constant as x varies in the bonded region and increases gradually from the crack tips to the axis of symmetry, where the maximum value is achieved. The displacement of Mg is larger than that of Mg17Al12, which is reasonable since Mg17Al12 is stiffer than Mg [15]. Tn reaches its maximum value at the crack tips and decreases rapidly to a nearly constant value for the bonded interface away from the crack tips. Tn is zero for the cracked interface. The overall distribution of Tn and un is generally consistent with the prediction of linear elastic fracture mechanics [20]. A similar observation was reported in previous studies on the fracture of Al/Si interfaces [12] and interfaces composed of BCC materials [21]. It can be noted from Figure 3 that the interface crack has propagated for a certain distance when ε increases from 0.031 to 0.046, in agreement with the observation in Figure 2. To render a general quantitative picture of the fracture process, the distribution of T n , u n , and ∆z along the interface for two strain values corresponding to Figure 2a,c is shown in Figure 3. Both T n and u n exhibit a near-symmetrical distribution with respect to an axis parallel to the z-axis. u n remains constant as x varies in the bonded region and increases gradually from the crack tips to the axis of symmetry, where the maximum value is achieved. The displacement of Mg is larger than that of Mg 17 Al 12 , which is reasonable since Mg 17 Al 12 is stiffer than Mg [15]. T n reaches its maximum value at the crack tips and decreases rapidly to a nearly constant value for the bonded interface away from the crack tips. T n is zero for the cracked interface. The overall distribution of T n and u n is generally consistent with the prediction of linear elastic fracture mechanics [20]. A similar observation was reported in previous studies on the fracture of Al/Si interfaces [12] and interfaces composed of BCC materials [21]. It can be noted from Figure 3 that the interface crack has propagated for a certain distance when ε increases from 0.031 to 0.046, in agreement with the observation in Figure 2.
The T n − u n curve was extracted from simulation, as shown in Figure 4. T n increases nearly linearly as u n increases, until reaching the maximum value which is known as cohesive strength σ c (corresponding u n is denoted by δ c ), followed by a descending trend to zero. The u n at which T n reduces to zero is δ s , indicating the complete separation of the interface. The T n − u n curve extracted from MD simulation was compared with that predicted by three widely used CZMs, namely cubic polynomial CZM [22], exponential CZM [23], and bilinear CZM [24]. In the cubic polynomial CZM, T n is related to u n by T n for the exponential CZM is expressed as σ c e 2 u n δ s exp − 16 9 e u n δ s In the bilinear CZM, is given as where u* = max (δ c , u n ). Note from Equations (1)-(3) that cubic polynomial and exponential CZMs are associated with two parameters (σ c and δ s ), while bilinear CZM is related to one more parameter δ c . These three parameters were obtained from MD simulation as σ c = 5.63 GPa, δ s = 7.80 Å, δ c = 1.53 Å, which were then input to Equations (1)-(3) to plot T n −u n curves, as shown in Figure 4. Exponential CZM exhibits the best agreement with MD results since it can capture the nonlinear processes associated with bond breaking. The Tn−un curve was extracted from simulation, as shown in Figure 4. Tn increases nearly linearly as un increases, until reaching the maximum value which is known as cohesive strength σc (corresponding un is denoted by δc), followed by a descending trend to zero. The un at which Tn reduces to zero is δs, indicating the complete separation of the interface. The Tn−un curve extracted from MD simulation was compared with that predicted by three widely used CZMs, namely cubic polynomial CZM [22], exponential CZM [23], and bilinear CZM [24]. In the cubic polynomial CZM, Tn is related to un by  The Tn−un curve was extracted from simulation, as shown in Figure 4. Tn increases nearly linearly as un increases, until reaching the maximum value which is known as cohesive strength σc (corresponding un is denoted by δc), followed by a descending trend to zero. The un at which Tn reduces to zero is δs, indicating the complete separation of the interface. The Tn−un curve extracted from MD simulation was compared with that predicted by three widely used CZMs, namely cubic polynomial CZM [22], exponential CZM [23], and bilinear CZM [24]. In the cubic polynomial CZM, Tn is related to un by

Conclusions
MD simulations were performed to investigate the fracture of Mg/Mg 17 Al 12 interface. The interface crack extended in a brittle manner without noticeable plasticity. The expressions for three CZMs were established by using parameters obtained from MD simulations. The normal traction-separation curves predicted by the three CZMs were compared with that extracted from MD simulations. Exponential CZM shows the closest match with MD results, followed by bilinear CZM.