Fracture Toughness Estimation of Single-Crystal Aluminum at Nanoscale

Abstract

In this publication, molecular dynamics simulations are used to investigate the fracture behavior of single-crystal aluminum. The stress intensity factor is estimated by means of four different methods, the accuracy is assessed for each approach and the fracture toughness is estimated. The proposed methodology is also applied to estimate the fracture toughness for graphene and diamond using published data from other scientific articles. The obtained fracture toughness for the single-crystal aluminum is compared with other nanomaterials that have similar microstructures. Dislocation emission during the fracture simulation of the cracked nano-crystal of aluminum is analyzed to study the fracture behavior. Brittle fracture behavior is the predominant failure mode for the nanomaterials studied in this research.

1. Introduction

Molecular dynamics (MD) simulations are successfully utilized to investigate tensile properties of nano-crystals such as ultimate tensile strength ($S_U$) and elastic modulus (E) [1,2,3,4,5,6]. Recently, cracked nano-crystals have been modeled by means of MD to investigate the fracture mechanics (FM) properties, which have been estimated for different nano-crystalline materials [7,8,9,10,11,12]. At the nanoscale, different methodologies have been proposed to evaluate fracture toughness by means of MD simulations of cracked nano-crystals [13,14,15,16,17,18,19,20,21,22,23,24]. Fracture toughness is the FM property that quantifies the material resistance to unstable crack propagation, which is usually obtained by means of the Energy Release Rate (G). This parameter (G) is estimated as the measure of the energy available for an increment of crack extension [25]. Using MD simulations, G is calculated from a cracked nano-crystal which is monotonically loaded until breaking [15,18,22,26]. However, the fracture toughness is established modeling a nano-crystal with a specific crack length and geometry that change from one investigation to other, therefore a fracture toughness that depends on the crack size is obtained, instead of a unique value that only depends on the material properties. The J-integral (J) is also an energy method that is extensively used to measure the fracture toughness [27,28,29]. Nevertheless, as other energy parameters, J lacks accuracy when irreversibility such as dislocations appear in the MD simulation [30]. In order to calculate J, Hardy’s approach is frequently implemented to estimate the local stress for atomistic systems as proposed by Zimmerman and Jones [31]. As in other atomistic investigations, Zimmerman and Jones estimate J modeling a nano-crystal with a unique initial crack length arbitrarily selected. Another method that is also widely utilized to estimate the fracture toughness is the crack tip opening displacement ($CTOD$). This FM parameter ($CTOD$) has been implemented by Thaulow et al. [32] to estimate the fracture toughness for silicon at the atomic scale. In addition, Thaulow et al. propose a method to estimate the plastic zone at the crack tip and the $CTOD$ from MD simulations. In other research, Skogsrud and Thaulow perform atomistic modeling for a bbc-Fe cantilever beam at 500 K, 800 K and 1200 K, and estimate the fracture toughness using $CTOD$, but also using a unique and arbitrary initial crack length [16]. Skogsrud and Thaulow obtain a fracture toughness two times larger than values reported in the literature for the same material. In other investigation, Skogsrud and Thaulow [33] show that different crystallographic orientations in the same cantilever beam affect the fracture behavior, varying from ductile fracture for $\left(10\overline{1}\right)\left[101\right]$ crack system to brittle for $\left(100\right)\left[0\overline{1}1\right]$. Additionally, Skogsrud and Thaulow report that $CTOD$ stops growing when dislocations are emitted from the crack tip during the MD simulation [33].

In the present work, based on MD simulations the fracture toughness $K_c$ is estimated for a single-crystal of Al by means of four different FM methods. The first method is the classical Irwing’s formulation for K. The second approach is a methodology proposed by Thaulow et al. [32]. The third method is based on a relationship between $CTOD$ and the stress of the system. Lastly, G is estimated to calculate $K_c$ using a relationship between G and K. In order to establish an independent $K_c$ on the crack length using the four aforementioned FM methods, a methodology is proposed based on the analysis of cracked single-crystals where the crack length is varied. The proposed methodology is validated using published experimental data for graphene and data from MD simulations for diamond. The obtained $K_c$ for Al single-crystal is compare with other nanomaterials that present similar microstructures. Finally, the results are discussed and a comparison of the FM methods is carried out.

2. Methodology

2.1. Molecular Dynamics Simulation

In this study, MD simulations are carried out using the code Large-scale Atomic/ Molecular Massively Parallel Simulator (LAMMPS) [34]. Single-crystals of Al with four different initial crack lengths are modeled. For each crystal, a uniaxial tensile test under controlled displacement is simulated. Initially, an isobaric-isothermal ensemble (NPT) is implemented to equilibrate the atomistic system using 20,000 steps. A deformation rate of 0.01% is applied to the single-crystal models in z-direction for 10,000 steps. The atomistic systems are relaxed using 20,000 steps (after every step deformation of 0.01%). A strain rate of $1\times {10}^{-4}$ /ps and time-step of 0.001 ps are selected. The described process (deforming-relaxing) is repeated until fracture of the single-crystal. The embedded atom method (EAM) proposed by Mendelev et al. [35] is implemented in order to study the fracture process. EAM potential has been used to investigate crack propagation in Al crystals by several authors [21,30,36,37,38,39,40,41,42]. All simulations are carried out at 300 K and 1.01 bar using a Nose/Hoover thermostat and barostat. During the deformation process, free surface condition is established in the x-direction and periodic conditions in the others. Approximately 200,000 atoms conformed the atomistic system with dimension of $24\times 8\times 16$ nm${}^3$ (Figure 1), with a lattice parameter $a=0.405$ nm for aluminum and initial crack lengths $l_0=5a,10a,15a,20a$. Virial stress tensor is utilized to compute the global stress as follow [43]:

$${\sigma }_{ij}=\frac{1}{V}\sum _{m=1}^N\left(\frac{1}{2}\sum _{m\ne n}^N{F}_{mn}^i{r}_{mn}^j-s_m{v}_m^i{v}_m^j\right)$$

(1)

where ${\sigma }_{ij}$ is the virial stress tensor, V is the atomistic system volume, $F_{mn}^i$ is the force vector between particle m and particle n, $r_{mn}^j$ is the distance vector between particle m and particle n, $s_m$ is the mass of the atom m and $v_m^i$ is the velocity vector of the atom m. The volume of the atoms is obtained using voro++ [44]. OVITO [45] is used to calculate the distance among atoms at the crack tip in order to estimate the $CTOD$. Using the dislocation extraction algorithm (DXA), the dislocations were analyzed [46].

2.2. Stress Intensity Factor Estimation

The stress intensity factor K is firstly estimated using the equation proposed by Irwin, as follows [47]:

$$K_{I,1}=f\sigma \sqrt{\pi l}$$

(2)

where $\sigma$ is the global stress, l is the crack length and f is the geometry factor for an edge crack, which is given by [47],

$$f=0.265{(1-\alpha )}^4+\frac{0.857+0.265\alpha }{{(1-\alpha )}^{3/2}}$$

(3)

where $\alpha =l/60a$. The same approach has been previously used to estimate $K_c$ from MD simulations in several publications [18,19,32,48].

Another approach considered in this research is an equation implemented by Thaulow [32] to estimate K for atomistic cracked elements. In this equation, K is estimated based on the $CTOD$ and the plastic zone size at the crack tip, as follows:

$$K_{I,2}=\frac{E\sqrt{2\pi }}{8}\frac{CTOD}{\sqrt{r_y}}$$

(4)

where the Young’s modulus is $E=60.58$ GPa for an Al single-crystal [30,49] and $r_y$ is the plastic zone radius, which is estimated as proposed by Thaulow in [32] (Figure 2). At the atomistic scale, some researchers relate the dislocation emission region with the plastic zone size [50]. These dislocation emissions are considered as a source of energy release before crack propagation [51], however dislocation emissions are not observed during the simulations of the cracked single-crystals (Figure 3).

Another equation to estimate K from $CTOD$ is given in [25], and it is defined by:

$$K_{I,3}=\sqrt{\frac{\pi ·E·CTOD·\sigma }{4·(1-\nu )}}$$

(5)

where the Poisson’s ratio is $\nu =0.36$ for an Al single-crystal [30]. This approach has been implemented previously in [18].

Finally, K was also computed based on G using the Equation [25]:

$$K_{I,4}=\sqrt{\frac{E}{1-{\nu }^2}·G}$$

(6)

where G is estimated as the dissipated energy during the fracture per unit of new crack area created. G is calculated as follow:

$$G=-\frac{\Delta (U-W)}{\Delta A}$$

(7)

where U is the potential energy of the system, $\Delta A$ is the new crack surface created due to fracture and W is the work applied during the tensile test before the fracture (${\epsilon }_f=6.7\%$, Figure 4). The initial work and potential energy are assumed to be negligible. Thus, W is estimated from the stress–strain curve (Figure 4) as follows [47]:

$$W=V{\int }_0^{{\epsilon }_f}\sigma d\epsilon$$

(8)

where V is the atomistic system volume. Jung et al. has used this method to estimate G from MD simulation [22].

3. Results and Discussion

3.1. Fracture Behavior

By means of dislocation extraction algorithm (DXA) [46], a dislocation analysis is carried out. The DXA identifies the dislocation lines in a nano-crystal in case that dislocations appear during the simulation. Based on this analysis, dislocation emissions are not evidenced during the fracture process (Figure 3), therefore brittle behavior is obtained in the MD simulations.

3.2. Stress Intensity Factor Assessment

The data used to calculate K by means of the different approaches, viz. Equations (2), (4), (5) and (6), are summarized in

Table 1. Data from the molecular dynamics (MD) simulations to estimate K.

${\mathit{l}}_0$ (m)	$\mathit{\alpha }$	f	$\mathit{\sigma }$ (GPa)	$\mathbf{CTOD}$ (m)	${\mathit{r}}_{\mathit{y}}$ (m)	W (fJ)	U (fJ)	$\mathbf{\Delta }\mathit{A}$ (m${}^2$)
2.03 $\times {10}^{-9}$	8.33 $\times {10}^{-2}$	1.19	2.96	0.137 $\times {10}^{-8}$	1.95 $\times {10}^{-9}$	0.62	0.21	1.61 $\times {10}^{-16}$
4.05 $\times {10}^{-9}$	1.67 $\times {10}^{-1}$	1.31	2.56	0.141 $\times {10}^{-8}$	3.81 $\times {10}^{-9}$	0.48	0.17	1.58 $\times {10}^{-16}$
6.08 $\times {10}^{-9}$	2.50 $\times {10}^{-1}$	1.51	2.38	0.140 $\times {10}^{-8}$	5.51 $\times {10}^{-9}$	0.44	0.13	1.34 $\times {10}^{-16}$
8.10 $\times {10}^{-9}$	3.33 $\times {10}^{-1}$	1.79	2.13	0.143 $\times {10}^{-8}$	7.79 $\times {10}^{-9}$	0.41	0.15	1.28 $\times {10}^{-16}$

. These values are computed just before the fracture. The K values are presented in

Table 2. Stress intensity factor just before fracture.

${\mathit{l}}_0$ (m)	${\mathit{K}}_{\mathit{I},1}$ (MPa$·\sqrt{\mathbf{m}}$)	${\mathit{K}}_{\mathit{I},2}$ (MPa$·\sqrt{\mathbf{m}}$)	${\mathit{K}}_{\mathit{I},3}$ (MPa$·\sqrt{\mathbf{m}}$)	${\mathit{K}}_{\mathit{I},4}$ (MPa$·\sqrt{\mathbf{m}}$)
2.03 $\times {10}^{-9}$	0.281	0.588	0.471	0.421
4.05 $\times {10}^{-9}$	0.378	0.433	0.443	0.369
6.08 $\times {10}^{-9}$	0.494	0.357	0.426	0.405
8.10 $\times {10}^{-9}$	0.606	0.307	0.407	0.375

.

3.3. Fracture Toughness for Al Single-Crystal

In order to establish a $K_c$ value that is independent on the crack length, an equivalent stress ${\sigma }_{eq}$ is computed based on the $K_I$ in Table 2. ${\sigma }_{eq}$ is calculated as:

$${\sigma }_{eq}=\frac{K_I}{f·\sqrt{\pi ·l_0}}$$

(9)

A dimensionless value is obtained dividing ${\sigma }_{eq}$ by $S_U$ for a single-crystal of Al without defects, viz., $S_U=6.182$ GPa [30].

Table 3. Dimensionless equivalent stress.

${\mathit{l}}_0$ (m)	${\mathit{\sigma }}_{\mathbf{eq}}/{\mathit{S}}_{\mathit{U}}\left({\mathit{K}}_{\mathit{I},1}\right)$	${\mathit{\sigma }}_{\mathbf{eq}}/{\mathit{S}}_{\mathit{U}}\left({\mathit{K}}_{\mathit{I},2}\right)$	${\mathit{\sigma }}_{\mathbf{eq}}/{\mathit{S}}_{\mathit{U}}\left({\mathit{K}}_{\mathit{I},3}\right)$	${\mathit{\sigma }}_{\mathbf{eq}}/{\mathit{S}}_{\mathit{U}}\left({\mathit{K}}_{\mathit{I},4}\right)$
2.03 $\times {10}^{-9}$	0.48	1.00	0.80	0.72
4.05 $\times {10}^{-9}$	0.41	0.47	0.48	0.40
6.08 $\times {10}^{-9}$	0.38	0.28	0.33	0.31
8.10 $\times {10}^{-9}$	0.34	0.17	0.23	0.21

summarizes ${\sigma }_{eq}/S_U$, which is computed using the data from Table 2. ${\sigma }_{eq}/S_U$ from Table 3 vs. $l_0/L$ (where $L=60a$) was plotted in Figure 5. Using the least squares method a curve based on Equation (9) is fitted to the data, thus $K_c$ is the $K_I$ that yields the curve with best fit (Figure 5). $K_c$ values estimated by means of the different approaches are presented in

Table 4. Fracture toughness for Al single-crystal.

${\mathit{K}}_{\mathit{c}}-\mathit{K}$ (MPa$·\sqrt{\mathbf{m}})$	${\mathit{K}}_{\mathit{c}}-\mathbf{CTOD}$-pz (MPa$·\sqrt{\mathbf{m}})$	${\mathit{K}}_{\mathit{c}}-\mathbf{CTOD}-\mathit{\sigma }$ (MPa$·\sqrt{\mathbf{m}})$	${\mathit{K}}_{\mathit{c}}-\mathit{G}$ (MPa$·\sqrt{\mathbf{m}})$
0.35	0.50	0.45	0.40

. The residual sum of squares (RSS) is implemented to determine the accuracy of the methods to predict $K_c$, and the results are summarized in

Table 5. Residual sum of squares for the different methods.

K	$\mathbf{CTOD}$-pz	$\mathbf{CTOD}-\mathit{\sigma }$	G
0.04871	0.05234	0.00260	0.00252

. The obtained lowest error is 0.00252, which corresponds to G method, however $CTOD-\sigma$ obtained almost the same error, viz., 0.00260. Additionally, $K_c$ is in the same order of magnitude reported by Chandra [36]. However, the fracture toughness reported by Chandra is determined modeling a single-crystal with a unique initial crack length that is arbitrarily chosen.

3.4. Fracture Toughness for Graphene and Diamond

The same methodology that is applied in Section 3.3 to establish $K_c$ for Al single-crystal is also used to estimate the $K_c$ for other two nanomaterials by analyzing data from published scientific articles.

Table 6. Fracture toughness and crack length for graphene and diamond.

Material	Crack Length (nm)	${\mathit{K}}_{\mathit{c}}$ (MPa$·\sqrt{\mathbf{m}})$	Reference
Graphene	33	3.1	[52]
	438	4.1
	518	3.7
	600	4.9
	1256	4.1
Diamond	1.785	8.360	[18]
	2.499	8.463
	3.213	8.405

presents the $K_c$ for cracked specimens of graphene and diamond with different crack lengths. The data for graphene are experimentally obtained in [52] and for the diamond by means of MD simulations in [18]. Using data from Table 6, ${\sigma }_{eq}$ is calculated by Equation (9) and normalized dividing ${\sigma }_{eq}$ by $S_U$, where $S_U=130$ GPa for graphene [53] and $S_U=241.2$ GPa for diamond [18]. Figure 6a,b show ${\sigma }_{eq}/S_U$ vs. $l_0/L$ and the respective adjusted curves to obtain the fracture toughness. The $K_c$ that is obtained for graphene and diamond based on the proposed methodology are shown in

Table 7. Fracture toughness for graphene and diamond and the respective error.

Material	${\mathit{K}}_{\mathit{c}}$ (MPa$·\sqrt{\mathbf{m}})$	E (GPa)	${\mathit{G}}_{\mathit{c}}$ (J/m${}^2$)	Error
Graphene	3.20	1000	10.3	1.03 $\times {10}^{-4}$
Diamond	8.40	1011.5	69.8	8.65 $\times {10}^{-6}$

with the respective errors which are calculated by means of the residual sum of squares.

3.5. Fracture Toughness Comparison

In order to compare the fracture toughness of the Al single-crystal to other nanomaterials with similar microstructures, $K_c$ vs. E is plotted in Figure 7. The different nanomaterials and their properties are summarized in

Table 8. FM properties for different nanomaterials.

Material	${\mathit{K}}_{\mathit{c}}$ (MPa$·\sqrt{\mathbf{m}})$	E (GPa)	${\mathit{G}}_{\mathit{c}}$ (J/m${}^2$)	Reference
Au	0.45	79	2.56	[54,55]
Cu	0.64	169.9	2.43	[56]
Ni	0.88	188	4.12	[57]

, where $G_c$ is estimated as follows:

$$G_c=K_c^2/E$$

(10)

The fracture toughness for graphene and diamond are also plotted in the same figure as reference values. As seen in Figure 7, materials with FCC microstructures show similar fracture behavior, viz.: Al, Au, Cu and Ni. It is important to point out that $K_c$ for Au, Cu and Ni are estimated in [54,55,56,57] based on modeling nano-crystals with a unique and arbitrary selected crack length.

4. Conclusions

In this research, the fracture toughness for Al single-crystals has been assessed successfully by means of four different approaches. The fracture toughness for graphene and diamond has been also estimated using published data by other researchers. Finally, $K_c$ for the Al single-crystal has been compared with nanomaterials that have microstructures similar to the Al. From the results, it is concluded:

• The fracture behavior observed in the MD simulations is in accordance with observations reported by other researchers in their investigations on similar Al crystals.
• Despite that Al is considered a ductile material at 300 K, a brittle fracture behavior is observed for single-crystals in the MD simulations.
• The methodology that is proposed in this research provides a suitable method to obtain a fracture toughness value that is independent on the crack length.
• The parameters G and $CTOD-\sigma$ yield a good accuracy to predict the fracture of single-crystals.
• Dislocations are not observed during the simulations. Therefore, it is not possible to compare the plastic zone estimations with the dislocation emission zone.