- Email: [email protected]
Application of CTOD in atomistic modeling of fracture
Engineering Fracture Mechanics 150 (2015) 153–160
Contents lists available at ScienceDirect
Engineering Fracture Mechanics journal homepage: www.elsevier.com/locate/engfracmech
Application of CTOD in atomistic modeling of fracture Jørn Skogsrud ⇑, Christian Thaulow Department of Engineering Design and Materials, Norwegian University of Science and Technology (NTNU), NO-7491 Trondheim, Norway
a r t i c l e
i n f o
Article history: Received 17 December 2014 Received in revised form 24 August 2015 Accepted 29 August 2015 Available online 1 October 2015 Keywords: Molecular dynamics Nanobeams CTOD Fracture
a b s t r a c t In this presentation we first review a method to use crack tip displacements and an estimate of the plastic zone to estimate fracture toughness. The same procedure was then applied on nanosized bcc-Fe cantilever beams with an initial crack. The results were then compared to literature. It was found that the fracture mechanical calculations should reflect the deformation processes at the crack tip in order to give accurate results. Further investigation is needed to determine whether the flange extension method can be used as a robust measure of fracture toughness. Ó 2015 Elsevier Ltd. All rights reserved.
1. Introduction In atomistic modeling of fracture, most examinations have been 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 used as a meaningful physical parameter in atomistic modeling of material failure. We begin with a short review of previous results from atomistic modeling of fracture of single crystal silicon [1], 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. We then expect our investigation to include iron and present results from atomistic modeling of iron cantilever beams 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. Earlier work on silicon The system for the earlier work on silicon was described in detail in earlier studies [1,2]. These earlier simulations of a crack in a defect-free silicon crystal, a first principles based ReaxFF reactive force field was applied, which have an accuracy close to that of quantum chemistry based techniques, even for bond breaking events [3]. The silicon lattice was oriented with 1 as the x-y-z-directions, with the fracture plane (1 0 0) and initial fracture direction [0 1 1]. The model geom½1 0 0½0 1 1 0 1 etry was approximately 200 Å by 180 Å, with a thickness of 15 Å and a crack length of 40 Å which gives about 27,500 atoms and (a/h) = 0.2. The loading was done with pure Mode I, at 500 K, 800 K and 1200 K. ⇑ Corresponding author. Tel.: +47 73593873. E-mail address: [email protected] (J. Skogsrud). http://dx.doi.org/10.1016/j.engfracmech.2015.08.043 0013-7944/Ó 2015 Elsevier Ltd. All rights reserved.
154
J. Skogsrud, C. Thaulow / Engineering Fracture Mechanics 150 (2015) 153–160
Nomenclature a B bij C ij CMOD CTOD d E F GG
cs h I KG L ry
rf w
crack depth elastic properties of material plane strain moduli elastic stiffness values Crack Mouth Opening Displacement Crack Tip Opening Displacement CTOD elastic modulus loading force critical Griffith stress for brittle cleavage surface energy cantilever height moment of inertia theoretical Griffith fracture toughness distance from cantilever notch to center of load point radius of plastic zone stress at the cantilever top surface cantilever width
In silicon the dominating mechanism changed from completely brittle bond breaking to ductile behavior with generation of dislocation as the temperature was raised from 500 K to 800 K and 1200 K. The 500 K simulation experienced crack propagation directly through the specimen along the initial crack plane, while the 1200 K simulation had two minor brittle jumps before dislocation emission along the {1 1 1} slip plane as the crack tip came to rest. The CTOD measurements and crack flange extrapolation method gave reasonably good estimates for the fracture toughpffiffiffiffiffi ness of silicon, with the calculated K at fracture initiation within 0.25 MPa m of the applied K at 800 K and 1200 K. Howpffiffiffiffiffi ever, at 500 K the calculated K was found to be 0.4 MPa m higher. One of the explanations for the increase was lattice trapping [4], which is dependent on temperature since thermal fluctuations are necessary to pull atoms out of the trapped state. Amorphization mechanisms [5] were also considered to be influential in the fracture initiation process. The same methodology used in silicon is now used, in this work, to calculate the K from the geometrical factors in iron, to see whether CTOD and the flange extrapolation method can be useful measurements in atomistic simulations of iron. 3. Fracture model The iron atoms were described by the Mendelev II potential, for which the physical properties relevant for this work are listed in Table 1. The geometry used was a cantilever beam with either a square or pentagonal cross-section, connected to a fixed sidewall. The pentagonal cross section was created by adding a right isosceles triangle to the bottom of the square cantilever with the base side facing the square. Periodic conditions were applied in the wall plane to allow for the cantilever to be attached to a bulk wall of atoms, as well as allowing the cantilever to flex outside the original dimensions of the simulation box. Non-periodic conditions were applied in the cantilever direction. An atomistically sharp initial crack with (a/h) = 0.3 and (a/h) = 0.5 was created by deactivating interactions between the atoms on each side of the crack plane surface. The crack was placed one beam width away from the sidewall, and the load-controlled loading was performed by adding force directly to a half spherical group of atoms centered on the beam 0.5 beam widths from the free end of the cantilever. The displacement was measured by tracking the average position of this group of atoms. The loading rate was set to 0.193 eV/Å ps, or 309 N/s. The system is initially relaxed for 0.3 ns, over 200,000 time steps in an isobaric-isothermal ensemble, and was subsequently kept at a temperature of 5 K through a Nose–Hoover thermostat in the canonical ensemble. The system contained approximately 40 million atoms, and the side width of the cantilever is 400 Å. An overview over the cantilever dimensions can be seen in Fig. 1. Table 2 has the parameters for the four cantilevers, named Cantilever A to D. The crack plane was (0 1 0) with the crack front oriented along the [0 0 1] direction for all the cantilevers. The simulations were run using the MD software LAMMPS [6] on the supercomputer at the The Norwegian metacenter for computational science (Notur) facilities. Compared to the silicon work, the iron model was three dimensional, larger, and has a different lattice orientation and crack depth. The simulations were also performed at considerably lower temperature, and with a simpler EAM potential compared to the ReaxFF reactive force field. The position of the atoms close to the crack tip can be measured directly in atomistic simulations, and the output can be used to calculate Crack Tip Opening Displacement (CTOD) and Crack Mouth Opening Displacement (CMOD) with very high accuracy. Similar to the simulations in Si by Thaulow [1], the crack flanges away from the crack tip are fairly straight, which makes it possible to do extrapolation from these. Fig. 2 visualizes the extrapolation from the crack flanges. The position of
J. Skogsrud, C. Thaulow / Engineering Fracture Mechanics 150 (2015) 153–160
155
Table 1 Relevant physical properties of the Mendelev II EAM Fe potential. Lattice constant (a) at T = 0 K Lattice constant (a) at T = 1820 K C 11 C 12 C 44 Eð100Þ (calculated) Eð110Þ (calculated)
2.8553 Å 2.926 Å 243.4 GPa 145.0 GPa 116.0 GPa 135.0 GPa 224.0 GPa 1.79 N/m 1.65 N/m
c100 c110
Fig. 1. Schematic overview over the system. The values for relevant system lengths are marked in the figure. At the left we see a 3D rendered image of the cantilever with pentagonal cross-section, and on the right is a 2D-slice of the cantilever with square cross-section.
Table 2 Simulation model parameters. Width (w) Simulation timestep and total time Height (h) Crack depth (a), (a/h) Wall width Wall to crack length Temperature Crack to load length (L) Orientations, crack plane
400 Å 1.5 fs timestep, 1.5 ns in total. 400 Å 120 Å (a/h) = 0.3 and 200 Å (a/h) = 0.5 150 Å 400 Å 5K 1000 Å {1 0 0}, (1 0 0) [0 0 1]
Cantilever Cantilever Cantilever Cantilever
a = 120 Å, a = 200 Å, a = 120 Å, a = 200 Å,
A B C D
pentagonal pentagonal square square
atoms close to the crack tip is tracked, seen as point A and B in Fig. 2, and is used to calculate two planes on each side of the crack. The planes are then moved along their normal until they are at the crack surface, and the extrapolation is done to find the intersection between them, marked m in Fig. 2. K has been calculated for square cantilever beams by a number of groups, among them Armstrong et al. [7] and Takashima et al. [8,9], who used the following form for KI
K I ¼ rf
pffiffiffiffiffiffi a paf h
where the dimensionless shape factor f
a
ð1Þ a h
is given by
a a2 a 3 a 4 ¼ 1:22 1:40 þ 7:33 f 13:08 þ 14:0 h h h h h
ð2Þ
For cantilevers with a pentagonal cross-section, the literature fitting the cantilevers simulated in this work the best was done by Zhao et al. [10] who used the same formula as Armstrong et al. for the KI, but approximated the stress as
rf ¼
FLy I
ð3Þ
156
J. Skogsrud, C. Thaulow / Engineering Fracture Mechanics 150 (2015) 153–160
Fig. 2. The figure illustrates how the measurements of CTOD, CMOD and the extrapolation from the crack flanges are done. Four evenly spaced rows of atoms in the (0 0 1)-direction (parallel to the crack front) on each side of the crack tip are tracked, and the CTOD and CMOD are tracked from the movement of these atoms. Necessary compensations are done in order to compensate for the fact that the tracked atoms are not exactly at the crack surface, as shown with the extrapolation line. The extrapolation is done between points A and B, marked in the figure, and a parallel line is drawn originating from the crack flange surface, indicated by the dotted line. The radius ry is also shown in the figure.
where the moment of inertia of the cantilever cross section was calculated by
I¼
2 2 3 wh h w4 w2 h þ ðh yÞ hw þ þ y þ 2 12 288 4 6
ð4Þ
and y is the vertical distance between the upper surface and the neutral plane, defined by
y¼
h2 w 2
2 þ w4 h þ w6 hw þ w4
2
ð5Þ
They calculated the dimensionless shape factor by
f
a
a a2 a3 a 4 a5 ¼ 1:122 1:121 þ 3:74 þ 3:873 19:05 þ 22:55 h h h h h h
ð6Þ
4. Results From the earlier work [1], we see that he values calculated according to the crack flange extrapolation procedure were found to be in good agreement with observations of the extension of the damaged structures as viewed on the atomistic models. The iron model experienced cleavage fracture along the (1 1 0) plane in the middle of the cantilever beam, but to varying degrees. The most brittle behavior was found in the pentagonal cross-section cantilever with the deepest crack (Cantilever B), while both shallower crack (Cantilever A) and square cross-section (Cantilever D) resulted in more ductile behavior. We observed some blunting at the crack tip before fracture initiation for all the cantilevers, as well as heavy blunting close to the sidewall for all simulations. The critical stress intensity factor K Ic at initiation of the first brittle fracture was calculated using both the classical equations and the CTOD model, Table 3. The classical equations were taken from Eqs. (1) and (2) for the square cross-section beams, and Eqs. (3), (4) and (6) for the pentagonal cross-section beams. The initial deformation mechanism observed in the cantilever beams were quite similar for the two different crosssections and crack lengths. After the elastic part of the loading curve, all the simulations showed bond-breaking events at the crack tip, where some atoms on the loading side of the crack moved a small distance relative to the rest of the atoms, compressing in the (1 1 0) direction downwards from the crack front. This corresponded to clearly visible movements in the CTOD-Displacement curve, Fig. 3, and happened just before dislocation emission in the pentagonal beams, while the time from bond breaking to dislocation emission was longer in the square cross-section beams. A visualization of the process can be seen in Fig. 4. At the onset of plastic deformation, the dislocations always originated from where the crack front met the
157
J. Skogsrud, C. Thaulow / Engineering Fracture Mechanics 150 (2015) 153–160 Table 3 Summary of the fracture mechanical results from the cantilever simulations. Cantilever
(a/h)
CTOD
2r y;Flange
A B
0.3 0.5
5.7 Å 5.4 Å
40 Å 42 Å
C D
0.3 0.5
5.7 Å 5.3 Å
34 Å 32 Å
pffiffiffiffiffi m) pffiffiffiffiffi 1.04 MPa m pffiffiffiffiffi 1.42 MPa m pffiffiffiffiffi 1.29 MPa m pffiffiffiffiffi 1.49 MPa m K I (MPa
[10] [10] [7] [7]
Fig. 3. CTOD as a function of displacement of the group of atoms used to load the cantilever beam. Bumps before the onset of plastic deformation are visible for all the cantilevers, and are marked as bond breaking events. The rapid upswing of CTOD corresponds to the first dislocation emission from the pillar sidewall.
Fig. 4. The figure shows the crack tip at the first bond-breaking events, and the initial dislocation emission. The BCC atoms are colored white, while nonBCC are gray. The picture to the top left shows the crack tip just before the first bond breaking, with a gray line drawn from one of the atoms at the crack tip 0)-plane. through two different directions. The lines are drawn trough the same atoms in the next pictures, showing sliding along the (1 1 0)-plane and (1 1 The top right picture shows the initial dislocation, which is marked in the CTOD-displacement curve Fig. 3, where the first dislocations nucleates from the sidewall, and start traveling along the crack front into the cantilever Fig. 5. The middle series of pictures shows the same timesteps, zoomed in, and the measurement of the CTOD is visualized.
158
J. Skogsrud, C. Thaulow / Engineering Fracture Mechanics 150 (2015) 153–160
Fig. 5. The figure shows a small cut from the cantilever, where BCC atoms are removed. The initial dislocation is colored red, and the magnified pictures showing how the dislocation is first created near the sidewall, and then travels into the beam along the crack front. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
direction), and the surface Fig. 6. The figure shows the formation of twinning close to the sides of the pillar. The view is slightly above (from the [0 1 3] atoms at the sidewall are removed to show the twinning plane marked in red. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
surface of the cantilever sidewall. These dislocations then connected the beam sidewall and the crack front, moving inwards toward the center of the crack, as shown in Fig. 5. The dislocations tended to be mixed mode dislocations, starting as pure edge dislocations at the cantilever sidewall, transforming to pure screw dislocations near the center of the cantilever. As the dislocations moved, the crack grew in the center of the beam, initially along the (1 1 0)-plane, while at the edges the crack blunted through emission of more dislocations. At higher strain, twinning was observed on the beam edges, traveling on {1 1 2}-planes in (1 1 1)-directions, as seen in Fig. 6. This is consistent with findings in similar systems [11,12]. Increased crack depth also resulted in fewer dislocations and more crack extension, as the shallower crack had significantly more dislocations at similar CTOD. 5. Discussion According to the Griffith concept, the theoretical fracture toughness is given by
KG ¼
qffiffiffiffiffiffiffiffiffiffiffiffiffi GG B1
ð7Þ
where GG ¼ 2cs and
vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 0sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi1ffi u u ub11 b22 @ b22 2b12 þ b66 A B¼t þ 2 b11 2b11
ð8Þ
and bij are the plane strain moduli obtained from the anisotropic elastic compliance [13]. Values for the fracture surface energy and the elastic modulus for the (1 0 0) plane in bcc iron from the potential are pffiffiffiffiffi c ¼ 1:79 N/m, E(1 0 0) = 1.35 1011 Pa and Poisson’s ratio of 0.3. This yields a Griffith fracture toughness K G of 0.92 MPa m [14] for the (1 0 0) [0 0 1] crack system. The equation from [15] shows the use of Eq. (9) to estimate the K from the plastic zone:
159
J. Skogsrud, C. Thaulow / Engineering Fracture Mechanics 150 (2015) 153–160
Fig. 7. The figure shows the extrapolation procedure from Fig. 2 on Cantilever A. The left part shows a 20 Å thick slice from the middle of the cantilever, with the atoms colored after von Mises stress. At the right we have the CTOD-displacement curve, with the arrow showing where on the plot the image is taken. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
Table 4 Comparison between the stress intensity for the Griffith criterion and the flange extrapolation method when using surface energy values for the (1 0 0)-plane and elastic modulus for the [1 0 0]-direction, and alternative calculations using values for the (1 1 0)-plane and [1 1 0]direction. Cantilever
GG (MPa
A B C D
0.90 0.90 0.90 0.90
KI ¼
pffiffiffiffiffi m)
Flange100 (MPa 0.54 0.50 0.58 0.56
pffiffiffiffiffiffiffi E 2p d pffiffiffiffi : ry 8
pffiffiffiffiffi m)
Flangealt (MPa
pffiffiffiffiffi m)
0.89 0.83 0.97 0.93
ð9Þ
The first step in calculating K I is to estimate 2ry by extrapolating the crack flanges. Fig. 7 shows an overview over how the extrapolation is done in Cantilever A, with a line being drawn through atom A and B in the figure, and a line parallel to this line is drawn at the crack surface. This is then done for both sides of the crack surface, and the distance between the crack tip and the intersection is taken as an estimate for 2ry . For the case visualized in Fig. 7, 2ry is estimated to be 40 Å and pffiffiffiffiffi CTOD = 5.7 Å. Applying Eq. (9) gives K Ic =0.54 MPa m, while calculations using the approach by Zhao et al. yield pffiffiffiffiffi 1.04 MPa m. When calculating CTOD in atomistic simulations, care must be taken to obtain accurate measurements that are consistent across orientations and crack widths. The approach we used was to track atoms close to the crack tip, so that we could get an expression for the surface plane, indicated by green1 lines in the middle series of Fig. 4. The CTOD is then measured from the center of one of the atoms near the crack tip, straight over the crack, to where it hits the surface (green line) on the other side. The calculated stress intensities are proportional to the CTOD squared as can be seen in Eq. (9), and this can reduce the accuracy when performing measurements on such small systems. When reviewing the silicon work, we observed that the crack extension happened on the same (0 1 1)-plane as the initial 0)-planes which were different from crack. In iron, all dislocation activity and crack extension happened on (1 1 0) and (1 1 the initial (1 0 0)-plane. The (1 0 0) direction in bcc iron had the lowest elastic modulus of any direction, and thus having the crack plane in the (1 0 0) plane and using the elastic modulus for this direction yielded the lowest estimate for the fracture toughness. Calculating the stress intensity using the values for the (1 1 0)-plane and [1 1 0] direction for surface energy and elastic modulus, which are 1.65 N/m and 224 GPa respectively, changes the results obtained from the flange extrapolation method, summarized in Table 4. The results obtained from using physical properties from the crystallographic direction and plane where the actual crack extension and dislocation movements happened, yielded values for the stress intensity closer to those found in earlier studies [12]. The stress intensity was calculated at the point when the crack front started to move, and since we then observed activity only on {1 1 0} planes, we used the physical properties tied to the active planes and directions in our calculations. We 1
For interpretation of color in Fig. 4, the reader is referred to the web version of this article.
160
J. Skogsrud, C. Thaulow / Engineering Fracture Mechanics 150 (2015) 153–160
also observe that the calculations from Armstrong and Zhao yields a larger variation in stress intensity across the different geometries than for the crack flange extrapolation method, and given that the initial crack across all geometries has the same geometry and similar size, the stress intensity was not expected to vary much. 6. Conclusion While the crack flange extension method has given reasonable results in earlier work in silicon, the results from applying this method to iron using the physical parameters for the nominal directions h1 0 0i and crack surfaces {1 0 0} were not as good. The crack flange extrapolation method gave K crit values far lower than the calculations from the continuum derived formulas, and while the results for K crit were closer to the Griffith criterion than the calculations from continuum mechanics, it was significantly lower than calculations from similar systems with the same potential [16,12]. However, observing the crack at failure, the deformation and dislocation activity was exclusively located in {1 1 0} planes, and using elastic modulus and surface energy from these planes and directions, the K I produced by the flange extrapolation method yielded results more in line with existing literature. This indicates that knowledge of the specific deformation mechanism may be necessary for the flange extrapolation method to be accurate, and further investigation is necessary before a robust procedure can be established. Acknowledgements The authors acknowledge the master student Marie Jørum for running the atomistic simulations, and NOTUR (Norwegian High Performance Computing) for providing the computational resources used in this work. We also thank the support from the EU project MultiHy (EU-FP7-NMP Project No. 263335). References [1] Thaulow C, Schieffer SV, Vatne IR, Sen D, Østby E. Crack tip opening displacement in atomistic modeling of fracture of silicon. Comput Mater Sci 2011;50(9):2621–7. [2] Sen D, Thaulow C, Schieffer SV, Cohen A, Buehler MJ. Atomistic study of crack-tip cleavage to dislocation emission transition in silicon single crystals. Phys Rev Lett 2010;104(23):235502. [3] Buehler MJ, Tang H, van Duin AC, Goddard III WA. Threshold crack speed controls dynamical fracture of silicon single crystals. Phys Rev Lett 2007;99 (16):165502. [4] Gumbsch P. Modelling brittle and semi-brittle fracture processes. Mater Sci Engng: A 2001;319321(0):1–7. http://dx.doi.org/10.1016/S0921-5093(01) 01062-0.. [5] Huang S, Zhang S, Belytschko T, Terdalkar SS, Zhu T. Mechanics of nanocrack: fracture, dislocation emission, and amorphization. J Mech Phys Solids 2009;57(5):840–50. http://dx.doi.org/10.1016/j.jmps.2009.01.006. . [6] Plimpton S. Fast parallel algorithms for short-range molecular dynamics. J Comput Phys 1995;117(1):1–19. [7] Armstrong D, Haseeb A, Roberts S, Wilkinson A, Bade K. Nanoindentation and micro-mechanical fracture toughness of electrodeposited nanocrystalline Ni–W alloy films. Thin Solid Films 2012;520(13):4369–72. [8] Takashima K, Higo Y, Sugiura S, Shimojo M. Fatigue crack growth behavior of micro-sized specimens prepared from an electroless plated Ni–P amorphous alloy thin film. Mater Trans-JIM 2001;42(1):68–73. [9] Takashima K, Higo Y, Swain M. Fatigue crack growth behaviour of micro-sized specimens prepared from amorphous alloy thin films. In: ICF10, Honolulu (USA) 2001; 2013. [10] Zhao X, Langford R, Tan J, Xiao P. Mechanical properties of sic coatings on spherical particles measured using the micro-beam method. Scripta Mater 2008;59(1):39–42. [11] Hora P, Pelikán V, Machová A, Spielmannová A, Prahl J, Landa M, et al. Crack induced slip processes in 3D. Engng Fract Mech 2008;75(12):3612–23. http://dx.doi.org/10.1016/j.engfracmech.2007.05.013. . [12] Ersland CH, Vatne IR, Thaulow C. Atomistic modeling of penny-shaped and through-thickness cracks in bcc iron. Model Simul Mater Sci Engng 2012;20 (7):075004. . [13] Lieberman DS, Zirinsky S. A simplified calculation for the elastic constants of arbitrarily oriented single crystals. Acta Crystallogr 1956;9(5):431–6. http://dx.doi.org/10.1107/S0365110X56003144. [14] Möller JJ, Bitzek E. Comparative study of embedded atom potentials for atomistic simulations of fracture in a-iron. Model Simul Mater Sci Engng 2014;22(4):045002. [15] Anderson TL, Anderson T. Fracture mechanics: fundamentals and applications. CRC Press; 2005. [16] Vatne IR, Østby E, Thaulow C, Farkas D. Quasicontinuum simulation of crack propagation in bcc-fe. Mater Sci Engng: A 2011;528(15):5122–34.
Contents lists available at ScienceDirect
Engineering Fracture Mechanics journal homepage: www.elsevier.com/locate/engfracmech
Application of CTOD in atomistic modeling of fracture Jørn Skogsrud ⇑, Christian Thaulow Department of Engineering Design and Materials, Norwegian University of Science and Technology (NTNU), NO-7491 Trondheim, Norway
a r t i c l e
i n f o
Article history: Received 17 December 2014 Received in revised form 24 August 2015 Accepted 29 August 2015 Available online 1 October 2015 Keywords: Molecular dynamics Nanobeams CTOD Fracture
a b s t r a c t In this presentation we first review a method to use crack tip displacements and an estimate of the plastic zone to estimate fracture toughness. The same procedure was then applied on nanosized bcc-Fe cantilever beams with an initial crack. The results were then compared to literature. It was found that the fracture mechanical calculations should reflect the deformation processes at the crack tip in order to give accurate results. Further investigation is needed to determine whether the flange extension method can be used as a robust measure of fracture toughness. Ó 2015 Elsevier Ltd. All rights reserved.
1. Introduction In atomistic modeling of fracture, most examinations have been 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 used as a meaningful physical parameter in atomistic modeling of material failure. We begin with a short review of previous results from atomistic modeling of fracture of single crystal silicon [1], 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. We then expect our investigation to include iron and present results from atomistic modeling of iron cantilever beams 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. Earlier work on silicon The system for the earlier work on silicon was described in detail in earlier studies [1,2]. These earlier simulations of a crack in a defect-free silicon crystal, a first principles based ReaxFF reactive force field was applied, which have an accuracy close to that of quantum chemistry based techniques, even for bond breaking events [3]. The silicon lattice was oriented with 1 as the x-y-z-directions, with the fracture plane (1 0 0) and initial fracture direction [0 1 1]. The model geom½1 0 0½0 1 1 0 1 etry was approximately 200 Å by 180 Å, with a thickness of 15 Å and a crack length of 40 Å which gives about 27,500 atoms and (a/h) = 0.2. The loading was done with pure Mode I, at 500 K, 800 K and 1200 K. ⇑ Corresponding author. Tel.: +47 73593873. E-mail address: [email protected] (J. Skogsrud). http://dx.doi.org/10.1016/j.engfracmech.2015.08.043 0013-7944/Ó 2015 Elsevier Ltd. All rights reserved.
154
J. Skogsrud, C. Thaulow / Engineering Fracture Mechanics 150 (2015) 153–160
Nomenclature a B bij C ij CMOD CTOD d E F GG
cs h I KG L ry
rf w
crack depth elastic properties of material plane strain moduli elastic stiffness values Crack Mouth Opening Displacement Crack Tip Opening Displacement CTOD elastic modulus loading force critical Griffith stress for brittle cleavage surface energy cantilever height moment of inertia theoretical Griffith fracture toughness distance from cantilever notch to center of load point radius of plastic zone stress at the cantilever top surface cantilever width
In silicon the dominating mechanism changed from completely brittle bond breaking to ductile behavior with generation of dislocation as the temperature was raised from 500 K to 800 K and 1200 K. The 500 K simulation experienced crack propagation directly through the specimen along the initial crack plane, while the 1200 K simulation had two minor brittle jumps before dislocation emission along the {1 1 1} slip plane as the crack tip came to rest. The CTOD measurements and crack flange extrapolation method gave reasonably good estimates for the fracture toughpffiffiffiffiffi ness of silicon, with the calculated K at fracture initiation within 0.25 MPa m of the applied K at 800 K and 1200 K. Howpffiffiffiffiffi ever, at 500 K the calculated K was found to be 0.4 MPa m higher. One of the explanations for the increase was lattice trapping [4], which is dependent on temperature since thermal fluctuations are necessary to pull atoms out of the trapped state. Amorphization mechanisms [5] were also considered to be influential in the fracture initiation process. The same methodology used in silicon is now used, in this work, to calculate the K from the geometrical factors in iron, to see whether CTOD and the flange extrapolation method can be useful measurements in atomistic simulations of iron. 3. Fracture model The iron atoms were described by the Mendelev II potential, for which the physical properties relevant for this work are listed in Table 1. The geometry used was a cantilever beam with either a square or pentagonal cross-section, connected to a fixed sidewall. The pentagonal cross section was created by adding a right isosceles triangle to the bottom of the square cantilever with the base side facing the square. Periodic conditions were applied in the wall plane to allow for the cantilever to be attached to a bulk wall of atoms, as well as allowing the cantilever to flex outside the original dimensions of the simulation box. Non-periodic conditions were applied in the cantilever direction. An atomistically sharp initial crack with (a/h) = 0.3 and (a/h) = 0.5 was created by deactivating interactions between the atoms on each side of the crack plane surface. The crack was placed one beam width away from the sidewall, and the load-controlled loading was performed by adding force directly to a half spherical group of atoms centered on the beam 0.5 beam widths from the free end of the cantilever. The displacement was measured by tracking the average position of this group of atoms. The loading rate was set to 0.193 eV/Å ps, or 309 N/s. The system is initially relaxed for 0.3 ns, over 200,000 time steps in an isobaric-isothermal ensemble, and was subsequently kept at a temperature of 5 K through a Nose–Hoover thermostat in the canonical ensemble. The system contained approximately 40 million atoms, and the side width of the cantilever is 400 Å. An overview over the cantilever dimensions can be seen in Fig. 1. Table 2 has the parameters for the four cantilevers, named Cantilever A to D. The crack plane was (0 1 0) with the crack front oriented along the [0 0 1] direction for all the cantilevers. The simulations were run using the MD software LAMMPS [6] on the supercomputer at the The Norwegian metacenter for computational science (Notur) facilities. Compared to the silicon work, the iron model was three dimensional, larger, and has a different lattice orientation and crack depth. The simulations were also performed at considerably lower temperature, and with a simpler EAM potential compared to the ReaxFF reactive force field. The position of the atoms close to the crack tip can be measured directly in atomistic simulations, and the output can be used to calculate Crack Tip Opening Displacement (CTOD) and Crack Mouth Opening Displacement (CMOD) with very high accuracy. Similar to the simulations in Si by Thaulow [1], the crack flanges away from the crack tip are fairly straight, which makes it possible to do extrapolation from these. Fig. 2 visualizes the extrapolation from the crack flanges. The position of
J. Skogsrud, C. Thaulow / Engineering Fracture Mechanics 150 (2015) 153–160
155
Table 1 Relevant physical properties of the Mendelev II EAM Fe potential. Lattice constant (a) at T = 0 K Lattice constant (a) at T = 1820 K C 11 C 12 C 44 Eð100Þ (calculated) Eð110Þ (calculated)
2.8553 Å 2.926 Å 243.4 GPa 145.0 GPa 116.0 GPa 135.0 GPa 224.0 GPa 1.79 N/m 1.65 N/m
c100 c110
Fig. 1. Schematic overview over the system. The values for relevant system lengths are marked in the figure. At the left we see a 3D rendered image of the cantilever with pentagonal cross-section, and on the right is a 2D-slice of the cantilever with square cross-section.
Table 2 Simulation model parameters. Width (w) Simulation timestep and total time Height (h) Crack depth (a), (a/h) Wall width Wall to crack length Temperature Crack to load length (L) Orientations, crack plane
400 Å 1.5 fs timestep, 1.5 ns in total. 400 Å 120 Å (a/h) = 0.3 and 200 Å (a/h) = 0.5 150 Å 400 Å 5K 1000 Å {1 0 0}, (1 0 0) [0 0 1]
Cantilever Cantilever Cantilever Cantilever
a = 120 Å, a = 200 Å, a = 120 Å, a = 200 Å,
A B C D
pentagonal pentagonal square square
atoms close to the crack tip is tracked, seen as point A and B in Fig. 2, and is used to calculate two planes on each side of the crack. The planes are then moved along their normal until they are at the crack surface, and the extrapolation is done to find the intersection between them, marked m in Fig. 2. K has been calculated for square cantilever beams by a number of groups, among them Armstrong et al. [7] and Takashima et al. [8,9], who used the following form for KI
K I ¼ rf
pffiffiffiffiffiffi a paf h
where the dimensionless shape factor f
a
ð1Þ a h
is given by
a a2 a 3 a 4 ¼ 1:22 1:40 þ 7:33 f 13:08 þ 14:0 h h h h h
ð2Þ
For cantilevers with a pentagonal cross-section, the literature fitting the cantilevers simulated in this work the best was done by Zhao et al. [10] who used the same formula as Armstrong et al. for the KI, but approximated the stress as
rf ¼
FLy I
ð3Þ
156
J. Skogsrud, C. Thaulow / Engineering Fracture Mechanics 150 (2015) 153–160
Fig. 2. The figure illustrates how the measurements of CTOD, CMOD and the extrapolation from the crack flanges are done. Four evenly spaced rows of atoms in the (0 0 1)-direction (parallel to the crack front) on each side of the crack tip are tracked, and the CTOD and CMOD are tracked from the movement of these atoms. Necessary compensations are done in order to compensate for the fact that the tracked atoms are not exactly at the crack surface, as shown with the extrapolation line. The extrapolation is done between points A and B, marked in the figure, and a parallel line is drawn originating from the crack flange surface, indicated by the dotted line. The radius ry is also shown in the figure.
where the moment of inertia of the cantilever cross section was calculated by
I¼
2 2 3 wh h w4 w2 h þ ðh yÞ hw þ þ y þ 2 12 288 4 6
ð4Þ
and y is the vertical distance between the upper surface and the neutral plane, defined by
y¼
h2 w 2
2 þ w4 h þ w6 hw þ w4
2
ð5Þ
They calculated the dimensionless shape factor by
f
a
a a2 a3 a 4 a5 ¼ 1:122 1:121 þ 3:74 þ 3:873 19:05 þ 22:55 h h h h h h
ð6Þ
4. Results From the earlier work [1], we see that he values calculated according to the crack flange extrapolation procedure were found to be in good agreement with observations of the extension of the damaged structures as viewed on the atomistic models. The iron model experienced cleavage fracture along the (1 1 0) plane in the middle of the cantilever beam, but to varying degrees. The most brittle behavior was found in the pentagonal cross-section cantilever with the deepest crack (Cantilever B), while both shallower crack (Cantilever A) and square cross-section (Cantilever D) resulted in more ductile behavior. We observed some blunting at the crack tip before fracture initiation for all the cantilevers, as well as heavy blunting close to the sidewall for all simulations. The critical stress intensity factor K Ic at initiation of the first brittle fracture was calculated using both the classical equations and the CTOD model, Table 3. The classical equations were taken from Eqs. (1) and (2) for the square cross-section beams, and Eqs. (3), (4) and (6) for the pentagonal cross-section beams. The initial deformation mechanism observed in the cantilever beams were quite similar for the two different crosssections and crack lengths. After the elastic part of the loading curve, all the simulations showed bond-breaking events at the crack tip, where some atoms on the loading side of the crack moved a small distance relative to the rest of the atoms, compressing in the (1 1 0) direction downwards from the crack front. This corresponded to clearly visible movements in the CTOD-Displacement curve, Fig. 3, and happened just before dislocation emission in the pentagonal beams, while the time from bond breaking to dislocation emission was longer in the square cross-section beams. A visualization of the process can be seen in Fig. 4. At the onset of plastic deformation, the dislocations always originated from where the crack front met the
157
J. Skogsrud, C. Thaulow / Engineering Fracture Mechanics 150 (2015) 153–160 Table 3 Summary of the fracture mechanical results from the cantilever simulations. Cantilever
(a/h)
CTOD
2r y;Flange
A B
0.3 0.5
5.7 Å 5.4 Å
40 Å 42 Å
C D
0.3 0.5
5.7 Å 5.3 Å
34 Å 32 Å
pffiffiffiffiffi m) pffiffiffiffiffi 1.04 MPa m pffiffiffiffiffi 1.42 MPa m pffiffiffiffiffi 1.29 MPa m pffiffiffiffiffi 1.49 MPa m K I (MPa
[10] [10] [7] [7]
Fig. 3. CTOD as a function of displacement of the group of atoms used to load the cantilever beam. Bumps before the onset of plastic deformation are visible for all the cantilevers, and are marked as bond breaking events. The rapid upswing of CTOD corresponds to the first dislocation emission from the pillar sidewall.
Fig. 4. The figure shows the crack tip at the first bond-breaking events, and the initial dislocation emission. The BCC atoms are colored white, while nonBCC are gray. The picture to the top left shows the crack tip just before the first bond breaking, with a gray line drawn from one of the atoms at the crack tip 0)-plane. through two different directions. The lines are drawn trough the same atoms in the next pictures, showing sliding along the (1 1 0)-plane and (1 1 The top right picture shows the initial dislocation, which is marked in the CTOD-displacement curve Fig. 3, where the first dislocations nucleates from the sidewall, and start traveling along the crack front into the cantilever Fig. 5. The middle series of pictures shows the same timesteps, zoomed in, and the measurement of the CTOD is visualized.
158
J. Skogsrud, C. Thaulow / Engineering Fracture Mechanics 150 (2015) 153–160
Fig. 5. The figure shows a small cut from the cantilever, where BCC atoms are removed. The initial dislocation is colored red, and the magnified pictures showing how the dislocation is first created near the sidewall, and then travels into the beam along the crack front. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
direction), and the surface Fig. 6. The figure shows the formation of twinning close to the sides of the pillar. The view is slightly above (from the [0 1 3] atoms at the sidewall are removed to show the twinning plane marked in red. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
surface of the cantilever sidewall. These dislocations then connected the beam sidewall and the crack front, moving inwards toward the center of the crack, as shown in Fig. 5. The dislocations tended to be mixed mode dislocations, starting as pure edge dislocations at the cantilever sidewall, transforming to pure screw dislocations near the center of the cantilever. As the dislocations moved, the crack grew in the center of the beam, initially along the (1 1 0)-plane, while at the edges the crack blunted through emission of more dislocations. At higher strain, twinning was observed on the beam edges, traveling on {1 1 2}-planes in (1 1 1)-directions, as seen in Fig. 6. This is consistent with findings in similar systems [11,12]. Increased crack depth also resulted in fewer dislocations and more crack extension, as the shallower crack had significantly more dislocations at similar CTOD. 5. Discussion According to the Griffith concept, the theoretical fracture toughness is given by
KG ¼
qffiffiffiffiffiffiffiffiffiffiffiffiffi GG B1
ð7Þ
where GG ¼ 2cs and
vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 0sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi1ffi u u ub11 b22 @ b22 2b12 þ b66 A B¼t þ 2 b11 2b11
ð8Þ
and bij are the plane strain moduli obtained from the anisotropic elastic compliance [13]. Values for the fracture surface energy and the elastic modulus for the (1 0 0) plane in bcc iron from the potential are pffiffiffiffiffi c ¼ 1:79 N/m, E(1 0 0) = 1.35 1011 Pa and Poisson’s ratio of 0.3. This yields a Griffith fracture toughness K G of 0.92 MPa m [14] for the (1 0 0) [0 0 1] crack system. The equation from [15] shows the use of Eq. (9) to estimate the K from the plastic zone:
159
J. Skogsrud, C. Thaulow / Engineering Fracture Mechanics 150 (2015) 153–160
Fig. 7. The figure shows the extrapolation procedure from Fig. 2 on Cantilever A. The left part shows a 20 Å thick slice from the middle of the cantilever, with the atoms colored after von Mises stress. At the right we have the CTOD-displacement curve, with the arrow showing where on the plot the image is taken. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
Table 4 Comparison between the stress intensity for the Griffith criterion and the flange extrapolation method when using surface energy values for the (1 0 0)-plane and elastic modulus for the [1 0 0]-direction, and alternative calculations using values for the (1 1 0)-plane and [1 1 0]direction. Cantilever
GG (MPa
A B C D
0.90 0.90 0.90 0.90
KI ¼
pffiffiffiffiffi m)
Flange100 (MPa 0.54 0.50 0.58 0.56
pffiffiffiffiffiffiffi E 2p d pffiffiffiffi : ry 8
pffiffiffiffiffi m)
Flangealt (MPa
pffiffiffiffiffi m)
0.89 0.83 0.97 0.93
ð9Þ
The first step in calculating K I is to estimate 2ry by extrapolating the crack flanges. Fig. 7 shows an overview over how the extrapolation is done in Cantilever A, with a line being drawn through atom A and B in the figure, and a line parallel to this line is drawn at the crack surface. This is then done for both sides of the crack surface, and the distance between the crack tip and the intersection is taken as an estimate for 2ry . For the case visualized in Fig. 7, 2ry is estimated to be 40 Å and pffiffiffiffiffi CTOD = 5.7 Å. Applying Eq. (9) gives K Ic =0.54 MPa m, while calculations using the approach by Zhao et al. yield pffiffiffiffiffi 1.04 MPa m. When calculating CTOD in atomistic simulations, care must be taken to obtain accurate measurements that are consistent across orientations and crack widths. The approach we used was to track atoms close to the crack tip, so that we could get an expression for the surface plane, indicated by green1 lines in the middle series of Fig. 4. The CTOD is then measured from the center of one of the atoms near the crack tip, straight over the crack, to where it hits the surface (green line) on the other side. The calculated stress intensities are proportional to the CTOD squared as can be seen in Eq. (9), and this can reduce the accuracy when performing measurements on such small systems. When reviewing the silicon work, we observed that the crack extension happened on the same (0 1 1)-plane as the initial 0)-planes which were different from crack. In iron, all dislocation activity and crack extension happened on (1 1 0) and (1 1 the initial (1 0 0)-plane. The (1 0 0) direction in bcc iron had the lowest elastic modulus of any direction, and thus having the crack plane in the (1 0 0) plane and using the elastic modulus for this direction yielded the lowest estimate for the fracture toughness. Calculating the stress intensity using the values for the (1 1 0)-plane and [1 1 0] direction for surface energy and elastic modulus, which are 1.65 N/m and 224 GPa respectively, changes the results obtained from the flange extrapolation method, summarized in Table 4. The results obtained from using physical properties from the crystallographic direction and plane where the actual crack extension and dislocation movements happened, yielded values for the stress intensity closer to those found in earlier studies [12]. The stress intensity was calculated at the point when the crack front started to move, and since we then observed activity only on {1 1 0} planes, we used the physical properties tied to the active planes and directions in our calculations. We 1
For interpretation of color in Fig. 4, the reader is referred to the web version of this article.
160
J. Skogsrud, C. Thaulow / Engineering Fracture Mechanics 150 (2015) 153–160
also observe that the calculations from Armstrong and Zhao yields a larger variation in stress intensity across the different geometries than for the crack flange extrapolation method, and given that the initial crack across all geometries has the same geometry and similar size, the stress intensity was not expected to vary much. 6. Conclusion While the crack flange extension method has given reasonable results in earlier work in silicon, the results from applying this method to iron using the physical parameters for the nominal directions h1 0 0i and crack surfaces {1 0 0} were not as good. The crack flange extrapolation method gave K crit values far lower than the calculations from the continuum derived formulas, and while the results for K crit were closer to the Griffith criterion than the calculations from continuum mechanics, it was significantly lower than calculations from similar systems with the same potential [16,12]. However, observing the crack at failure, the deformation and dislocation activity was exclusively located in {1 1 0} planes, and using elastic modulus and surface energy from these planes and directions, the K I produced by the flange extrapolation method yielded results more in line with existing literature. This indicates that knowledge of the specific deformation mechanism may be necessary for the flange extrapolation method to be accurate, and further investigation is necessary before a robust procedure can be established. Acknowledgements The authors acknowledge the master student Marie Jørum for running the atomistic simulations, and NOTUR (Norwegian High Performance Computing) for providing the computational resources used in this work. We also thank the support from the EU project MultiHy (EU-FP7-NMP Project No. 263335). References [1] Thaulow C, Schieffer SV, Vatne IR, Sen D, Østby E. Crack tip opening displacement in atomistic modeling of fracture of silicon. Comput Mater Sci 2011;50(9):2621–7. [2] Sen D, Thaulow C, Schieffer SV, Cohen A, Buehler MJ. Atomistic study of crack-tip cleavage to dislocation emission transition in silicon single crystals. Phys Rev Lett 2010;104(23):235502. [3] Buehler MJ, Tang H, van Duin AC, Goddard III WA. Threshold crack speed controls dynamical fracture of silicon single crystals. Phys Rev Lett 2007;99 (16):165502. [4] Gumbsch P. Modelling brittle and semi-brittle fracture processes. Mater Sci Engng: A 2001;319321(0):1–7. http://dx.doi.org/10.1016/S0921-5093(01) 01062-0.