- Email: [email protected]
Application of CTOD in Atomistic Modeling of Fracture
Available online at www.sciencedirect.com
ScienceDirect Procedia Materials Science 3 (2014) 1542 – 1547
20th European Conference on Fracture (ECF20)
Application of CTOD in atomistic modeling of fracture Christian Thaulow* Department of Engineering Design and Materials, The Norwegian University of Science and Technology, NO-7491 Trondheim, Norway
Abstract In atomistic modeling of fracture, most examinations are based on the Griffith energy, calculated from the surface energy, and the stress intensity factor derived from continuum mechanical formulations. Both approaches, however, reflect continuum mechanics and lack the direct connection to the atomistic crack tip mechanisms, and it is thus appealing to examine if Crack Tip Opening Displacement, CTOD, can be a meaningful physical parameter in atomistic modeling of material failure. In this presentation we first review our previous examinations of the crack tip displacements from fracture of single crystal silicon. We there established a procedure to determine the deformation zone ahead of the crack tip, based on CTOD90degree and CTODblunting, and demonstrated that CTOD could give a realistic estimate of the fracture toughness. The same procedure has then been applied on iron, loaded by a Modified Boundary Layer (MBL) model. In this case the crack flanges were not straight, and we were not able to extrapolate the flanges and define the center of rotation. However, by directly measure the deformed zone ahead of the crack tip, we could obtain realistic estimates of the fracture toughness via the CTOD. © 2014 Elsevier Ltd. This is an open access article under the CC BY-NC-ND license © 2014 The Authors. Published by Elsevier Ltd. (http://creativecommons.org/licenses/by-nc-nd/3.0/). Selectionand andpeer-review peer-review under responsibility the Norwegian University of Science and Technology (NTNU), Department of Selection under responsibility of theofNorwegian University of Science and Technology (NTNU), Department Structural of StructuralEngineering. Engineering Keywords: : Fracture toughness; CTOD; Silicon; Iron; Atomistic modeling
* Corresponding author. Tel.: +47 93059460 E-mail address: [email protected]
2211-8128 © 2014 Elsevier Ltd. This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/licenses/by-nc-nd/3.0/). Selection and peer-review under responsibility of the Norwegian University of Science and Technology (NTNU), Department of Structural Engineering doi:10.1016/j.mspro.2014.06.249
Christian Thaulow / Procedia Materials Science 3 (2014) 1542 – 1547
1543
1. Introduction In atomistic modeling of fracture, most examinations are based on the Griffith energy, calculated directly from the surface energy, and the stress intensity factor derived from continuum mechanical formulations. Both approaches, however, reflect continuum mechanics and lack the direct connection to the atomistic crack tip mechanisms, and it is thus appealing to examine if CTOD can be a meaningful physical parameter in atomistic modeling of material failure. The plan for the paper is first to review our previous results from atomistic modeling of fracture of single crystal silicon, Thaulow (2011), where the linear-elastic Griffith analysis was extended by including the mechanisms responsible for the formation of the deformation zone ahead of the crack tip. In the case of silicon this included covalent bond breaking, lattice trapping and amorphization. Next, we have extended our research to include iron and we present results from atomistic modeling of iron crystals and a procedure to estimate CTOD based on direct measurements of geometry and the deformation zone ahead of the crack tip at initiation of fracture. 2. Fracture Models Silica.The system set up for the work reported here is the same as described in our earlier studies, Sen (2010). We apply the first principles based ReaxFF reactive force field, which retains nearly the accuracy of quantum mechanics, even for bond breaking events, Buehler (2007). A perfect crystal with an edge crack is loaded by displacing the boundaries as shown in Figure 1. The crystal is oriented in the x-y-z directions, [100]-[011]-[0 1], creating a fracture plane (100) with initial fracture direction [011]. The model geometry is approximately 200Å by 180Å, with a thickness of 15Å and a crack length of 40Å. The number of atoms is approximately 27.500. We use periodic boundary conditions in the x and z directions. We note that this constraint may influence the possibility of nucleating certain type of dislocations from the tip of the crack.
Figure 1. Silicon single crystal under Mode I loading with an edge crack of {100}<110> character (left). Schematic presentation of CTOD blunting, CTOD 90degree and the deformation zone 2ry (right).
Iron. The iron single crystal is described by the Mendelev II potential, with a lattice parameter of 2.8553 Å, Mendelev (2003). The surface energies of the {100} and {110} planes are 1.79 J/m2 and 1.65 J/m2. The model employs a modified boundary layer (MBL) approach, where the analytical solution for the stress field around the crack tip is used as boundary conditions. To account for the anisotropy of bcc-Fe an anisotropic formulation of the MBL was implemented, Vatne (2011).
1544
Christian Thaulow / Procedia Materials Science 3 (2014) 1542 – 1547
The geometry is a square model with an edge crack. The crack plane is (010) with the crack front oriented along the [101] direction. Three dimensions were investigated, from a width of 1472 Å to 10.000 Å, by applying the Quasi Continuum approach, Miller (2007).
3. Results By shifting the temperature from 500 K to 800 K and further to 1200 K, the dominating fracture mechanism in silicon changes from completely brittle bond breaking to ductile behavior with generation of dislocations. A recent publication by the authors revealed that the emergence of very small geometric irregularities formed along the crack front was the prerequisite for the transition from brittle to ductile appearance, Sen (2010). As the critical temperature is approached, crack tip blunting dominates and is accompanied by changes in the perfect structure at the crack tip with bond rotations that lead to the formation of new bond structures that differ from the hexagonal rings typically found in silicon. In computational fracture mechanics the CTOD is defined according to displacement at the original crack tip or alternatively at the intersection of a 90 degree vertex with the crack flanges. We refer to these alternatives as CTODblunting and CTOD90degree, Figure 1 (right). We identify the surface atoms along the crack numerically and calculate the crack profiles. These profiles could successfully be fitted to a parabolic function. Based on this function both CTOD90degree, CTODblunting and the radius of the deformation zone, ry could be estimated for each load step. The total deformation zone, 2ry, is defined as the distance from the crack tip (determined by CTOD90degree) and the hinge rotation center (determined by CTODblunting), Figure 2. The values calculated according to this procedure were found to be in good agreement with observations of the extension of the damaged structures as viewed on the atomistic models.
Figure 2. Comparison between atomistic model, left, and calculation, right, at the initiation of brittle fracture of silicon at temperature 500K.
The iron model experienced cleavage fracture along the (010) plane. We observe some blunting at the crack tip before fracture initiation, Figure 3.The critical stress intensity factor KIc was calculated according to the deformation in the MBL model at fracture initiation, K Ic=0.94 MPam1/2.
1545
Christian Thaulow / Procedia Materials Science 3 (2014) 1542 – 1547
Figure 3. Crack tip in iron just before initiation of cleavage fracture. The color contours refers to strain energy density.
4 Discussion In this section we compare the measured CTOD with the stress intensity factor K calculated for an edge crack, Figure 1, and the crack resistance from the Griffith surface energy concept. According to the Griffith criterion, Griffith (1920), cleavage is expected when
K cleave
2J S E 1X 2
(1)
Typical values for the fracture surface energy and the elastic modulus for the (100) plane in silicon are J 100 2.41N / m , E(100)=1.40∙1011 Pa, Poisons ratio is 0.28. This yields a Griffith fracture toughness of 0.86 MPam1/2. The stress intensity factor for the edge crack in the silicon model with crack size a=40Å is expressed as, Tada (1985)
KI
§
V Sa ¨ 0.265 1 D 4 ¨ ©
0.857 0.265D ·¸ 3 ¸ (1 D ) 2 ¹
(2)
where α is the ratio between the crack size and the total length of the model. In the present case, α=0.22. Fracture initiation is observed at an applied stress intensity of KI=0.83MPam1/2 at 800K and 1200K, and KI=1.25MPam1/2 at 500 K. Hence, for the 800K and 1200K the applied stress intensity is close to the Griffith resistance, KI=0.86MPam1/2 while at 500K the stress intensity exceeds the resistance with 0.4MPam 1/2. We now also calculate the stress intensity factor KI by converting the measured CTOD to KI, [10]
KI
E 2S G 8 ry
(3)
where ry is the radius of the plastic zone at the crack tip. The calculations at initiation of the brittle fracture are based on the direct measurements of the CTOD and the size of the deformed zone, Figure 2. At 500 K, ry=4.7 Å and fracture is predicted to initiate at KI=1.52 MPam1/2, somewhat higher than KI=1.25 MPam1/2 from the continuum mechanics calculations, Eq. (2). For the higher temperatures the KI is lower, about KI=1.05 MPam1/2 but still higher than KI=0.83MPam1/2 .
1546
Christian Thaulow / Procedia Materials Science 3 (2014) 1542 – 1547
To summarize, the calculations based on CTOD increases the K with about 0.25MPam1/2 for all temperatures. According to other reports amorphization will increase the resistance to failure, Buehler (2007) and Kermode (2008), and we therefore propose that the increase of toughness compared with the Griffith calculation is due to this amorphization mechanism. An additional increase in stress intensity of about 0.4MPam1/2 is observed at 500K, and it is proposed that this is reflecting the lattice trapping mechanism, Gumbsch (2001). The lattice trapping phenomenon is related to the different length scales involved in cutting the bonds at the crack tip and elastically relaxing the crack tip region. Even if enough energy is stored to make it possible for the crack to advance, atoms just ahead at the crack tip are not pulled far enough from each other for the bonds between them to break. This is possible because the Griffith stored energy is ahead of the crack tip, while lattice trapping refers to the forces right at the crack tip region. The effect is temperature dependent, because thermal fluctuations should provide energy needed to pull the atoms out of the trapped state. We now apply the same procedure on the iron. The first approach was to estimate 2r y by extrapolating the crack flanges, Fig 4. But the flanges are not straight and it is difficult to perform an extrapolation. For the case visualized in Figure 4, ry was estimated to be 30 Å and CTOD=7.6 Å. Applying equation (3) gives KIc=0.60 MPam1/2, while the value from the MBL model was 0.93 MPam1/2. The lack of straight flanges violates this procedure, so alternatively we have measured the deformation zone directly at the crack tip just before crack initiation, Figure 3. The deformation zone is about two times the CTOD, hence we estimate ry=15 Å, with a corresponding KIc=0.84 MPam1/2. We are still underestimating the calculated stress intensity, but the results are encouraging and open up for further investigations.
Figure 4. Definitions of CTOD (left) and 2ry estimation by straight extrapolation of the crack flanges
Conclusions We demonstrate that measuring the Crack Tip Opening Displacement (CTOD) of single crystal silicon and iron at atomic level provides realistic estimate of the fracture toughness.
Christian Thaulow / Procedia Materials Science 3 (2014) 1542 – 1547
1547
The CTOD is a pure geometry related parameter and yields an independent measure of the fracture toughness in addition to the Griffith approach. CTOD is also capable to extend the linear-elastic approach into more non-linear behavior. Different strategies had to be applied for silicon and iron, and further examinations are necessary in order to establish robust procedures.
Acknowledgements: The author acknowledges the cooperation with the authors of Thaulow (2011) and also the atomistic calculations of iron crystals performed by Heidi Gryteland Holm during her master thesis at NTNU. The support from NOTUR (Norwegian High Performance Computing) and the cooperation with the Arctic Materials project at SINTEF are also acknowledged. References Griffith, A.A., The phenomena of rupture and flow in solids. Philosophical Transactions, Series A, 1920. 221: p. 163-198. Thaulow, C ; Scieffer, S V; Vatne, I R; Sen D and Østby E (2011) “Crack Tip Opening Displacement in Atomistic Modeling of Fracture of Silicon”, Journ Computational Materials Science 50 (2011) 2621–2627 Mendelev, M.I., Han S.,Srolovitz, D. J.,Ackland, G. J., Sun, D. Y. and Asta M (2003). Development of new interatomic potentials appropriate for crystalline and liquid iron. Phil Mag, 83 (35), 3977-3994. Vatne, Inga R; Østby E; Thaulow,C and Farkas,D (2011), “Quasicontinuum simulation of crack propagation in bcc-Fe”. Materials Science and Engineering A Buehler, M., et al., Threshold crack speed controls dynamical fracture of silicon single crystals. Physical Review Letters, 2007. 99(16): p. 165502. Sen, D., et al., Atomistic Study of Crack-Tip Cleavage to Dislocation Emission Transition in Silicon Single Crystals. Physical Review Letters, 2010. 104(23). Tada, H., P. Paris, and G. Irwin, The stress analysis handbook. Paris Productions (and Del Research Corp.) St. Louis, MO, 1985. Kermode, J., et al., Low-speed fracture instabilities in a brittle crystal. Nature, 2008. 455(7217): p. 1224-1227. Gumbsch, P., Modelling brittle and semi-brittle fracture processes. Materials Science and Engineering A, 2001. 319: p. 1-7. Miller, R. E. and Tadmor, E. B. (2007) QC Reference Manual version 1.3. Downloaded from www.qcmethod.com
ScienceDirect Procedia Materials Science 3 (2014) 1542 – 1547
20th European Conference on Fracture (ECF20)
Application of CTOD in atomistic modeling of fracture Christian Thaulow* Department of Engineering Design and Materials, The Norwegian University of Science and Technology, NO-7491 Trondheim, Norway
Abstract In atomistic modeling of fracture, most examinations are based on the Griffith energy, calculated from the surface energy, and the stress intensity factor derived from continuum mechanical formulations. Both approaches, however, reflect continuum mechanics and lack the direct connection to the atomistic crack tip mechanisms, and it is thus appealing to examine if Crack Tip Opening Displacement, CTOD, can be a meaningful physical parameter in atomistic modeling of material failure. In this presentation we first review our previous examinations of the crack tip displacements from fracture of single crystal silicon. We there established a procedure to determine the deformation zone ahead of the crack tip, based on CTOD90degree and CTODblunting, and demonstrated that CTOD could give a realistic estimate of the fracture toughness. The same procedure has then been applied on iron, loaded by a Modified Boundary Layer (MBL) model. In this case the crack flanges were not straight, and we were not able to extrapolate the flanges and define the center of rotation. However, by directly measure the deformed zone ahead of the crack tip, we could obtain realistic estimates of the fracture toughness via the CTOD. © 2014 Elsevier Ltd. This is an open access article under the CC BY-NC-ND license © 2014 The Authors. Published by Elsevier Ltd. (http://creativecommons.org/licenses/by-nc-nd/3.0/). Selectionand andpeer-review peer-review under responsibility the Norwegian University of Science and Technology (NTNU), Department of Selection under responsibility of theofNorwegian University of Science and Technology (NTNU), Department Structural of StructuralEngineering. Engineering Keywords: : Fracture toughness; CTOD; Silicon; Iron; Atomistic modeling
* Corresponding author. Tel.: +47 93059460 E-mail address: [email protected]
2211-8128 © 2014 Elsevier Ltd. This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/licenses/by-nc-nd/3.0/). Selection and peer-review under responsibility of the Norwegian University of Science and Technology (NTNU), Department of Structural Engineering doi:10.1016/j.mspro.2014.06.249
Christian Thaulow / Procedia Materials Science 3 (2014) 1542 – 1547
1543
1. Introduction In atomistic modeling of fracture, most examinations are based on the Griffith energy, calculated directly from the surface energy, and the stress intensity factor derived from continuum mechanical formulations. Both approaches, however, reflect continuum mechanics and lack the direct connection to the atomistic crack tip mechanisms, and it is thus appealing to examine if CTOD can be a meaningful physical parameter in atomistic modeling of material failure. The plan for the paper is first to review our previous results from atomistic modeling of fracture of single crystal silicon, Thaulow (2011), where the linear-elastic Griffith analysis was extended by including the mechanisms responsible for the formation of the deformation zone ahead of the crack tip. In the case of silicon this included covalent bond breaking, lattice trapping and amorphization. Next, we have extended our research to include iron and we present results from atomistic modeling of iron crystals and a procedure to estimate CTOD based on direct measurements of geometry and the deformation zone ahead of the crack tip at initiation of fracture. 2. Fracture Models Silica.The system set up for the work reported here is the same as described in our earlier studies, Sen (2010). We apply the first principles based ReaxFF reactive force field, which retains nearly the accuracy of quantum mechanics, even for bond breaking events, Buehler (2007). A perfect crystal with an edge crack is loaded by displacing the boundaries as shown in Figure 1. The crystal is oriented in the x-y-z directions, [100]-[011]-[0 1], creating a fracture plane (100) with initial fracture direction [011]. The model geometry is approximately 200Å by 180Å, with a thickness of 15Å and a crack length of 40Å. The number of atoms is approximately 27.500. We use periodic boundary conditions in the x and z directions. We note that this constraint may influence the possibility of nucleating certain type of dislocations from the tip of the crack.
Figure 1. Silicon single crystal under Mode I loading with an edge crack of {100}<110> character (left). Schematic presentation of CTOD blunting, CTOD 90degree and the deformation zone 2ry (right).
Iron. The iron single crystal is described by the Mendelev II potential, with a lattice parameter of 2.8553 Å, Mendelev (2003). The surface energies of the {100} and {110} planes are 1.79 J/m2 and 1.65 J/m2. The model employs a modified boundary layer (MBL) approach, where the analytical solution for the stress field around the crack tip is used as boundary conditions. To account for the anisotropy of bcc-Fe an anisotropic formulation of the MBL was implemented, Vatne (2011).
1544
Christian Thaulow / Procedia Materials Science 3 (2014) 1542 – 1547
The geometry is a square model with an edge crack. The crack plane is (010) with the crack front oriented along the [101] direction. Three dimensions were investigated, from a width of 1472 Å to 10.000 Å, by applying the Quasi Continuum approach, Miller (2007).
3. Results By shifting the temperature from 500 K to 800 K and further to 1200 K, the dominating fracture mechanism in silicon changes from completely brittle bond breaking to ductile behavior with generation of dislocations. A recent publication by the authors revealed that the emergence of very small geometric irregularities formed along the crack front was the prerequisite for the transition from brittle to ductile appearance, Sen (2010). As the critical temperature is approached, crack tip blunting dominates and is accompanied by changes in the perfect structure at the crack tip with bond rotations that lead to the formation of new bond structures that differ from the hexagonal rings typically found in silicon. In computational fracture mechanics the CTOD is defined according to displacement at the original crack tip or alternatively at the intersection of a 90 degree vertex with the crack flanges. We refer to these alternatives as CTODblunting and CTOD90degree, Figure 1 (right). We identify the surface atoms along the crack numerically and calculate the crack profiles. These profiles could successfully be fitted to a parabolic function. Based on this function both CTOD90degree, CTODblunting and the radius of the deformation zone, ry could be estimated for each load step. The total deformation zone, 2ry, is defined as the distance from the crack tip (determined by CTOD90degree) and the hinge rotation center (determined by CTODblunting), Figure 2. The values calculated according to this procedure were found to be in good agreement with observations of the extension of the damaged structures as viewed on the atomistic models.
Figure 2. Comparison between atomistic model, left, and calculation, right, at the initiation of brittle fracture of silicon at temperature 500K.
The iron model experienced cleavage fracture along the (010) plane. We observe some blunting at the crack tip before fracture initiation, Figure 3.The critical stress intensity factor KIc was calculated according to the deformation in the MBL model at fracture initiation, K Ic=0.94 MPam1/2.
1545
Christian Thaulow / Procedia Materials Science 3 (2014) 1542 – 1547
Figure 3. Crack tip in iron just before initiation of cleavage fracture. The color contours refers to strain energy density.
4 Discussion In this section we compare the measured CTOD with the stress intensity factor K calculated for an edge crack, Figure 1, and the crack resistance from the Griffith surface energy concept. According to the Griffith criterion, Griffith (1920), cleavage is expected when
K cleave
2J S E 1X 2
(1)
Typical values for the fracture surface energy and the elastic modulus for the (100) plane in silicon are J 100 2.41N / m , E(100)=1.40∙1011 Pa, Poisons ratio is 0.28. This yields a Griffith fracture toughness of 0.86 MPam1/2. The stress intensity factor for the edge crack in the silicon model with crack size a=40Å is expressed as, Tada (1985)
KI
§
V Sa ¨ 0.265 1 D 4 ¨ ©
0.857 0.265D ·¸ 3 ¸ (1 D ) 2 ¹
(2)
where α is the ratio between the crack size and the total length of the model. In the present case, α=0.22. Fracture initiation is observed at an applied stress intensity of KI=0.83MPam1/2 at 800K and 1200K, and KI=1.25MPam1/2 at 500 K. Hence, for the 800K and 1200K the applied stress intensity is close to the Griffith resistance, KI=0.86MPam1/2 while at 500K the stress intensity exceeds the resistance with 0.4MPam 1/2. We now also calculate the stress intensity factor KI by converting the measured CTOD to KI, [10]
KI
E 2S G 8 ry
(3)
where ry is the radius of the plastic zone at the crack tip. The calculations at initiation of the brittle fracture are based on the direct measurements of the CTOD and the size of the deformed zone, Figure 2. At 500 K, ry=4.7 Å and fracture is predicted to initiate at KI=1.52 MPam1/2, somewhat higher than KI=1.25 MPam1/2 from the continuum mechanics calculations, Eq. (2). For the higher temperatures the KI is lower, about KI=1.05 MPam1/2 but still higher than KI=0.83MPam1/2 .
1546
Christian Thaulow / Procedia Materials Science 3 (2014) 1542 – 1547
To summarize, the calculations based on CTOD increases the K with about 0.25MPam1/2 for all temperatures. According to other reports amorphization will increase the resistance to failure, Buehler (2007) and Kermode (2008), and we therefore propose that the increase of toughness compared with the Griffith calculation is due to this amorphization mechanism. An additional increase in stress intensity of about 0.4MPam1/2 is observed at 500K, and it is proposed that this is reflecting the lattice trapping mechanism, Gumbsch (2001). The lattice trapping phenomenon is related to the different length scales involved in cutting the bonds at the crack tip and elastically relaxing the crack tip region. Even if enough energy is stored to make it possible for the crack to advance, atoms just ahead at the crack tip are not pulled far enough from each other for the bonds between them to break. This is possible because the Griffith stored energy is ahead of the crack tip, while lattice trapping refers to the forces right at the crack tip region. The effect is temperature dependent, because thermal fluctuations should provide energy needed to pull the atoms out of the trapped state. We now apply the same procedure on the iron. The first approach was to estimate 2r y by extrapolating the crack flanges, Fig 4. But the flanges are not straight and it is difficult to perform an extrapolation. For the case visualized in Figure 4, ry was estimated to be 30 Å and CTOD=7.6 Å. Applying equation (3) gives KIc=0.60 MPam1/2, while the value from the MBL model was 0.93 MPam1/2. The lack of straight flanges violates this procedure, so alternatively we have measured the deformation zone directly at the crack tip just before crack initiation, Figure 3. The deformation zone is about two times the CTOD, hence we estimate ry=15 Å, with a corresponding KIc=0.84 MPam1/2. We are still underestimating the calculated stress intensity, but the results are encouraging and open up for further investigations.
Figure 4. Definitions of CTOD (left) and 2ry estimation by straight extrapolation of the crack flanges
Conclusions We demonstrate that measuring the Crack Tip Opening Displacement (CTOD) of single crystal silicon and iron at atomic level provides realistic estimate of the fracture toughness.
Christian Thaulow / Procedia Materials Science 3 (2014) 1542 – 1547
1547
The CTOD is a pure geometry related parameter and yields an independent measure of the fracture toughness in addition to the Griffith approach. CTOD is also capable to extend the linear-elastic approach into more non-linear behavior. Different strategies had to be applied for silicon and iron, and further examinations are necessary in order to establish robust procedures.
Acknowledgements: The author acknowledges the cooperation with the authors of Thaulow (2011) and also the atomistic calculations of iron crystals performed by Heidi Gryteland Holm during her master thesis at NTNU. The support from NOTUR (Norwegian High Performance Computing) and the cooperation with the Arctic Materials project at SINTEF are also acknowledged. References Griffith, A.A., The phenomena of rupture and flow in solids. Philosophical Transactions, Series A, 1920. 221: p. 163-198. Thaulow, C ; Scieffer, S V; Vatne, I R; Sen D and Østby E (2011) “Crack Tip Opening Displacement in Atomistic Modeling of Fracture of Silicon”, Journ Computational Materials Science 50 (2011) 2621–2627 Mendelev, M.I., Han S.,Srolovitz, D. J.,Ackland, G. J., Sun, D. Y. and Asta M (2003). Development of new interatomic potentials appropriate for crystalline and liquid iron. Phil Mag, 83 (35), 3977-3994. Vatne, Inga R; Østby E; Thaulow,C and Farkas,D (2011), “Quasicontinuum simulation of crack propagation in bcc-Fe”. Materials Science and Engineering A Buehler, M., et al., Threshold crack speed controls dynamical fracture of silicon single crystals. Physical Review Letters, 2007. 99(16): p. 165502. Sen, D., et al., Atomistic Study of Crack-Tip Cleavage to Dislocation Emission Transition in Silicon Single Crystals. Physical Review Letters, 2010. 104(23). Tada, H., P. Paris, and G. Irwin, The stress analysis handbook. Paris Productions (and Del Research Corp.) St. Louis, MO, 1985. Kermode, J., et al., Low-speed fracture instabilities in a brittle crystal. Nature, 2008. 455(7217): p. 1224-1227. Gumbsch, P., Modelling brittle and semi-brittle fracture processes. Materials Science and Engineering A, 2001. 319: p. 1-7. Miller, R. E. and Tadmor, E. B. (2007) QC Reference Manual version 1.3. Downloaded from www.qcmethod.com