- Email: [email protected]
Atomistic study of lattice trapping behavior for brittle fracture in bcc-iron
Computational Materials Science 40 (2007) 376–381 www.elsevier.com/locate/commatsci
Atomistic study of lattice trapping behavior for brittle fracture in bcc-iron Ya-Fang Guo a
a,*
, Chong-Yu Wang
b
Institute of Engineering Mechanics, Beijing Jiaotong University, Beijing 100044, China b Department of Physics, Tsinghua University, Beijing 100084, China
Received 31 December 2005; received in revised form 27 November 2006; accepted 29 January 2007 Available online 16 April 2007
Abstract The lattice trapping behavior for brittle fracture in bcc-iron has been studied by atomistic simulations. A pronounced anisotropy for brittle cleavage fracture of a mode I crack is observed in the discrete atomistic scale, and the preferred direction for cleavage is along the h1 1 0i direction on both {1 0 0} and {0 1 1} planes. The analysis of the atomic structure indicates that, due to the shear effect at the crack tip, the stacking faults or partial dislocations are formed before crack cleavage occurs. For the crack with a h1 1 0i crack front, shear occurs easily along the slip direction of bcc crystals, whereas for the crack with a h1 0 0i front, shear can only occur along the non-slip direction, and a strong lattice trapping is exhibited. We conclude that the anisotropy for cleavage fracture and the lattice trapping behavior are closely related with the slip systems of bcc crystals. 2007 Elsevier B.V. All rights reserved. PACS: 71.15.Pd; 62.20.Mk Keywords: Atomistic simulations; Lattice trapping; Anisotropy; Brittle fracture; Slip systems; Bcc-iron
1. Introduction The macroscopic failure of materials is ultimately determined by events occurring on the atomic scale, it is therefore a typical multiscale phenomenon in physics. Traditionally, engineering fracture mechanics, which is based on the continuum mechanics, is used to study material fracture properties, and many useful criterions are obtained. However, it cannot give a fracture criterion from the physical point of view. For example, based on the Griffith criterion for materials brittle fracture [1], it can be predicted that a perfectly brittle crack in a crystal is expected to choose a cleavage plane with low surface energy and to propagate on this plane with equal ease in all directions. But this theoretical prediction shows a clear deviation from
*
Corresponding author. Fax: +86 10 51682094. E-mail address: [email protected] (Y.-F. Guo).
0927-0256/$ - see front matter 2007 Elsevier B.V. All rights reserved. doi:10.1016/j.commatsci.2007.01.019
the experimental observations [2,3] that a crack may not cleavage along the low energy surface preferentially and sometimes a very pronounced anisotropy is observed for cleavage fracture. Moreover, results from atomistic simulations [4] indicate that the anisotropy for cleavage fracture and the lattice trapping behavior are related with atoms bonding, which can not be fully explained by the continuum theory. Lattice trapping effect was proposed about 30 years ago by Thomson et al. [5] when the first atomistic study of fracture was applied using a simple pair-wise potential, and the results indicated that crack propagation was ultimately determined by atoms bonding at the crack tip on the discrete atomic scale. Owing to the lattice trapping, a crack remains stable and not to advance or heal until loads somewhat larger or smaller than the Griffith load are reached. The critical loads for a crack to advance or heal are called the upper and lower trapping limits. Further studies [6,7] showed that the magnitude of the lattice trapping effect
Y.-F. Guo, C.-Y. Wang / Computational Materials Science 40 (2007) 376–381
strongly depends on atoms bonding at the crack tip and the forces needed to separate bonds successively. Subsequently, studies of brittle fracture in Silicon by experiments and simulations [8,9] revealed that the anisotropy for cleavage fracture was strongly related with the lattice trapping effect. Then the ab-initio study of the cleavage fracture [8] indicated that the magnitude of the trapping range decreases with the increasing of the system size. Moreover, the effect of boundary conditions on the lattice trapping behavior for crack propagation in bcc-iron was studied by atomistic simulations [10], and many types of interatomic potentials were used [10–12]. Furthermore, a combined finite element and atomistic simulation method (FEAt) [12] and the molecular static simulation [13] were both used to study the anisotropy for cleavage fracture and the lattice trapping behavior in bcc-iron. From above results, we find that the anisotropy of cleavage fracture and lattice trapping behavior have been widely investigated by atomic scale simulations, and the atoms bonding at the crack tip in particular has been studied. In this paper, the analysis of the slip systems for bcc crystals is introduced. Based on the observation of the atomic structure of the crack tip at the critical load for different types of crack in bcc-iron, we aim to find the relationship between the anisotropy behavior for cleavage fracture and the slip systems for bcc crystals. 2. Simulation model and procedures As shown in Fig. 1, the model used in our study includes two regions: the interior atomistic region I containing the crack tip and the boundary region II. The outer most two layers of atoms in the atomistic region are identified as the boundary region, and its thickness is larger than the cut-off distance of the Finnis–Sinclair N-body potential [14], which has been proved to be efficient for bcc crystal calculations. We consider the plane strain problem and apply the periodic boundary condition in Z direction. The lengths of the atomistic region in X and Y directions ˚ , and the initial crack length is about 37 A ˚. are both 350 A In atomistic simulations we first get the pre-existing crack by shifting the atoms from their positions in a perfect
377
crystal to positions specified by the anisotropic elastic continuum mechanics equations [14] for a certain load. According to anisotropic linear elastic continuum theory, the displacement fields of an open mode I crack for the plane strain condition [15] is given by fug ¼ K I ðrÞ
1=2
fgðhÞg;
ð1Þ
where T
fug ¼ fux ; uy g ;
fgg ¼ fgx ; gy g
T
ð2Þ
and K I ¼ rðapÞ
1=2
ð3Þ
In Eqs. (1)–(3), KI is the stress intensity factor; r and h are the polar coordinates of the atom for which the displacement is calculated; g is a function of the angle h between ~ r and the crack plane; r is the applied load; and a is the length of the crack. In our study we use the stress intensity factor KI to express the corresponding external load r, which can be defined as ‘nKIC’, where KIC is the critical stress intensity factor. In the present work, {0 1 0}h1 0 1i{1 1 0}h1 0 1i,{1 1 0} h0 0 1i and {0 1 0}h0 0 1i mode I (opening mode) cracks in a bcc-iron crystal are studied, where {0 1 0} and {1 1 0} represent the cleavage planes of the crack, h1 0 1i and h0 0 1i are the directions of the crack front. After getting the initial configurations of the crack at different loads and different temperatures, the atoms at the crack tip are fully relaxed for 2000 steps of magnitude 5 · 1015 s with the fixed-displacement boundary condition. The temperature of the system is invariant throughout the simulation, which is achieved by scaling all atoms instantaneous velocities with the appropriate Maxwell–Boltzmann distribution at a specified temperature. By observing atoms bonding at the crack tip, the critical load for crack extension or heal is obtained. 3. Results and discussion 3.1. Calculations of the surface energy Using the F–S N-body potential we firstly calculated the surface energies of the {1 0 0} and {1 1 0} planes for bcciron. The calculated values are: cð100Þ ¼ 1:90 J=m2 cð110Þ ¼ 1:70 J=m
2
and thus cð100Þ > cð110Þ :
Fig. 1. Scheme of the calculation model: region I is the atomistic region, and region II is the border region.
We notice that the calculated surface energy of the {1 0 0} plane is higher than that of the {1 1 0} plane, which is also in good agreement with the previous calculation results by using the same interatomic potential [12]. However, the experimentally determined surface energy of 2.32 J/m2 [16] is rather higher than our calculated values. Moreover, first principle calculations give values of the
378
Y.-F. Guo, C.-Y. Wang / Computational Materials Science 40 (2007) 376–381
surface energy that are close to the experimental value [17– 19], and the {1 0 0} energy value of 2.29 J/m2 is almost identical to the {1 1 0} value of 2.29 J/m2 [17]. For the consistency of the results, we use our calculated values in the present simulation. 3.2. Lattice trapping behavior for crack propagation After getting the initial configurations of the crack at different loads according to anisotropic linear elastic continuum mechanics equations, atoms at the crack tip are fully relaxed with the fixed-displacement boundary condition. Then the critical loads for crack extension at different temperatures are obtained as shown in Table 1. In the previous studies, it has been found that the upper and lower trapping limits are sensitive to the simulation conditions [10,12,20] such as the type of atomistic potential, the boundary condition, the crystallography and the temperature. Table 1 lists our simulation results. It is shown that the calculated upper trapping limits K þ I and the lower trapping limits K I for the {0 1 0}h1 0 1i crack in bcc-iron by using the F–S potential are very close to the theoretical critical stress intensity factor KIC according to the Griffith criterion, while cleavage occurs along the {1 0 0} planes. Comparing with the results obtained by fixed stressed-boundary condition and flexible-boundary condition [10], we find that the displacement-border condition will be valid for simulating crack propagation if the system size is selected sufficiently large. From the results for different types of cracks in Table 1, a strong directional anisotropy for the cleavage of a model I crack is observed. A crack with the h0 0 1i crack front produces pronounced lattice trapping, while a crack with the h0 1 1i front propagates at about the Griffith load. This means that cracks with the h0 0 1i crack front present significantly more lattice trapping than cracks with the h1 1 0i crack front. Moreover, the strongest trapping effects are observed for cleavage along the {1 1 0} planes and it is found that lattice trapping strongly favors cleavage along the {1 0 0} planes. In Table 2 we have compared our results with previous studies. Although the different interatomic potential and atomistic simulation method are used, similar trends are observed that the {1 0 0} plane is the preferred cleavage plane in bcc-iron, and the critical load for a {0 1 0}h1 0 1i crack starting propagation is the lowest compared with that of other cracks. However, there still exist differences between the above calculated values. For examTable 1 The critical load (KI/KIC) for crack propagation at different temperatures Temperature {0 1 0}h1 0 1i {0 1 0}h1 0 1i {1 1 0}h1 1 0i {1 1 0}h0 0 1i {0 1 0}h1 0 0i
K+ K K+ K+ K+
5K
100 K
200 K
300 K
400 K
600 K
1.11 0.89 1.28 1.5 2.00
1.09 0.90 1.25 1.48 1.80
1.07 0.94 1.18 1.46 1.72
1.04 0.98 1.15 1.45 1.71
1.04 0.99 1.13 1.44 1.65
1.05 1.01 1.13 1.40 1.64
Table 2 A comparison with previous calculations Crack type {0 1 0}h1 0 1i {1 1 0}h1 1 0i {1 1 0}h0 0 1i {0 1 0}h1 0 0i
K+ K+ K+ K+
Our work
Ref. [11]
Ref. [12]a
Ref. [12]b
1.11 1.28 1.5 2.00
1.0 1.8 2.4 2.0
1.11 1.72 1.19 1.41
1.07 1.62 1.17 1.75
The critical load for crack propagation is expressed as KI/KIC. a Using the Finnis–Sinclair potential. b Using the Johnson potential.
ple, our results agree well with the results in Ref. [11] that the {1 1 0}h1 1 0i crack can propagate more easily than the {1 1 0}h0 0 1i and {0 1 0}h1 0 0i cracks. Whereas in Ref. [12] the {1 1 0}h0 0 1i crack is predicted to propagate more easily. According to the Griffith theory for brittle fracture, a plane with the lowest surface energy should be the preferred cleavage plane, and a crack will propagate on this plane with equal ease in all directions. In our calculation, a strong anisotropy for crack cleavage exists, and the critical load for a {0 1 0}h1 0 1i crack starting propagation is the lowest compared with that of other cracks, which is also consistent with the experimental results that {1 0 0} is the most likely cleavage plane of iron [21]. However, our calculated surface energy of the {1 0 0} plane is about 12% higher than that of the {1 1 0} plane, and first principle calculations show the {1 0 0} and {1 1 0} surfaces have similar stabilities because the {1 0 0} energy value is almost identical to the {1 1 0} value. This means that the lattice trapping in the discrete atomistic scale causes the anisotropy for cleavage, and the preferred plane for cleavage may not be with the lowest surface energy. Furthermore, we notice that the differences between K þ I , K I and KIC at high temperatures are smaller than those at lower temperatures, which indicates that the trapping effect becomes weak at high temperatures and the trapping barrier can be overcome by thermal activation. It is also shown that the upper and lower trapping limits are almost invariant with the temperature increasing when the temperature is higher than 200 K. As we have known, crack propagation is ultimately determined by atoms bonding at the crack tip, thus we give the relationship between the length of the bond at a crack tip and the load for different types of cracks in Fig. 2. We observe that, when the load is smaller than the upper trapping limits, the lengths of the bond at the crack tip almost have no change with the load increasing for all types of cracks. But if the load is greater than the upper trapping limit, the length of the bond increases rapidly due to the bond breaking, and the crack extends subsequently.
3.3. Lattice trapping effect and slip systems of the crystal Fig. 3 is the critical atomic configuration of a crack tip at 5 K just before the crack tip bond breaking. The solid
Y.-F. Guo, C.-Y. Wang / Computational Materials Science 40 (2007) 376–381
Fig. 2. The relationship between the length of the bond at the crack tip and the load; , indicates the {0 1 0}h1 0 1i crack, indicates the {1 1 0}h1 0 1i crack, n indicates the {1 1 0}h0 0 1i crack, h indicates the {0 1 0}h0 0 1i crack.
and open circles represent atoms in the layer A and B, respectively along the periodic direction. In Fig. 3 we observe that shear occurs at the crack tip at the critical load. For {0 1 0}h1 0 1i and {1 1 0}h1 0 1i cracks in Fig. 3a and b, shear occurs along the h1 1 1i direction, while the stacking fault (SF) is formed at the crack tip. Whereas for {1 1 0}h0 0 1i and {0 1 0}h0 0 1i cracks in Fig. 3c and d,
379
shear occurs along the h0 1 0i direction, which causes the generation of partial dislocations. For bcc crystals, it is well known that the slip planes are {1 1 0}, {2 1 1} or {321}, where the slip direction is always in the h1 1 1i direction. But at low temperatures, the slip system is mainly as the {2 1 1} plane with the h1 1 1i slip direction. As shown in Fig. 4a, the crack fronts for {0 1 0}h1 0 1i and {1 1 0}h1 0 1i cracks are both perpendicular to the {1 0 1} plane (ABCD). Thus the atoms in the {1 0 1} plane can move along the h1 1 1i slip direction, which causes the formation of the SF at the crack tip. The angles a and b between the slip direction and the crack cleavage plane are 35.26 and 54.74, respectively. The slip plane is {2 1 1}. But for {1 1 0} h0 0 1i and {0 1 0}h0 0 1i cracks, the crack fronts are perpendicular to the {0 0 1} plane (ABEF) as shown in Fig. 4b, and the atoms in the {0 0 1} plane cannot move along the h1 1 1i slip direction due to the constraint of the plane strain condition. Because the resistance for atoms movement along the non-slip direction is much greater than that along the slip system in a crystal, thus for the crack with a h0 0 1i front, shear can only occur along the non-slip directions h1 1 0i and h0 1 0i to create partial dislocations when the load is sufficiently large. Therefore, the lattice trapping for the crack with a h0 0 1i front is much stronger than that of the crack with a h1 1 0i front, and the anisotropy for brittle cleavage fracture
Fig. 3. Atomic configuration of the crack tip at the critical load: (a) {0 1 0}h1 0 1i crack, KI = 1.11KIC; (b) {1 1 0}h1 0 1i crack, KI = 1.27KIC; (c) {1 1 0}h0 0 1i crack, KI = 1.49KIC; (d) {0 1 0}h0 0 1i crack, KI = 1.9KIC.
380
Y.-F. Guo, C.-Y. Wang / Computational Materials Science 40 (2007) 376–381
Fig. 4. Lattice structure of a bcc crystal: (a) for the {0 1 0}h1 0 1i and {1 1 0}h1 0 1i cracks; (b) for the {1 1 0}h0 0 1i and {0 1 0}h0 0 1i cracks.
is exhibited, which is closely related with the slip systems of the crystal. Moreover, we notice in Fig. 3 that shear occurs for all types of cracks before crack propagation, which causes the SF formation or partial dislocations creation. We have to mention that the shear would change atoms bonding at the crack tip, and it is closely correlated with lattice trapping behavior. But the precise analysis requires further the first principle calculations.
Fig. 5 gives atoms displacements at a crack tip after relaxation. Open circles represent atoms positions before relaxation, and solid circles represent atoms positions after relaxation. We see clearly that there are more atoms moving away from their initial positions in Fig. 5c and d than those in Fig. 5a and b, and that the magnitudes of atoms displacements in Fig. 5c and d are also greater. Table 3 lists the energy of atomistic region at the critical load for crack propagation. Although the average energy of atoms are
Fig. 5. Atoms displacements at the crack tip after relaxation: (a) {0 1 0}h1 0 1i crack, KI = 1.11KIC; (b) {1 1 0}h1 0 1i crack, KI = 1.27KIC; (c) {1 1 0}h0 0 1i crack, KI = 1.49KIC; and (d) {0 1 0}h0 0 1i crack, KI = 1.9KIC; hollow circles indicate atoms positions before relaxation, and solid circles indicate atoms positions after relaxation.
Y.-F. Guo, C.-Y. Wang / Computational Materials Science 40 (2007) 376–381 Table 3 The energy of atomistic region at the critical load for crack propagation at 5K Crack type
Number of atoms
Eeatom (eV)
Etot (eV)
{0 1 0}h1 0 1i {1 1 0}h1 1 0i {1 1 0}h0 0 1i {0 1 0}h1 0 0i
138,240 138,240 97,200 98,304
4.2606 4.2603 4.2603 4.2545
588,985 588,943 414,101 418,234
almost the same for four types of cracks, the total energy has greater difference between cracks with the h1 0 0i front and h1 1 0i front. For the crack with a h1 0 0i front, higher strain energy is stored at the crack tip before the crack cleavage occurring. Therefore, a higher critical load for the crack propagation is needed for bond breaking at the crack tip, and a strong lattice trapping is found for the crack with a h0 0 1i crack front. 4. Conclusions The lattice trapping behavior and the anisotropy for cleavage fracture are closely related with the slip systems of bcc crystals. For a mode I crack in bcc-iron, the preferred direction for cleavage is along the h1 1 0i direction on both {1 0 0} and {0 1 1} planes, and a stronger lattice trapping is found for the crack with a h1 0 0i crack front than for the crack with a h1 1 0i front. Due to the shear effect at the crack tip, the SF is formed at the crack tip along the slip direction of the bcc crystal for the crack with a h1 1 0i front, just before crack cleavage occurs. While for the crack with a h1 0 0i front, shear can only occur along the non-slip directions and partial dislocations are generated. Meanwhile, high strain energy is found to be stored at the crack tip for the crack with a h0 0 1i front, and a higher critical load for crack propagation is needed for atoms bond breaking, thus a strong lattice trapping is exhibited.
381
Acknowledgements The research was supported by ‘‘973’’ Project from the Ministry of Science and Technology of China (Grant No. 2006CD605102), Chinese Nature Science Foundation (Grant No. 10672016) and NJTU Science Foundation (Grant No. 2005SM0035). References [1] A.A. Griffith, Philos. Trans. Roy. Soc. A 221 (1920) 163. [2] A. George, G. Michot, Mater. Sci. Eng. A 164 (1993) 118. [3] J. Riedle, P. Gumbsch, H.F. Fischmeister, Phys. Rev. Lett. 76 (1996) 3594. [4] P. Gumbsch, Mater. Sci. Eng. A 319–321 (2001) 1. [5] R. Thomson, C. Hsieh, V. Rana, J. Appl. Phys. 42 (1971) 3154. [6] J.E. Sinclair, Philos. Mag. 31 (1975) 647. [7] W.A. Curtin, J. Mater. Res. 5 (1990) 1549. [8] R. Perez, P. Gumbsch, Acta Mater. 48 (2000) 4517. [9] J.C.H. Spence, Y.M. Huang, O. Sankey, Acta Metall. Mater. 41 (1993) 2815. [10] B. DeCelis, A.S. Argon, S. Yip, J. Appl. Phys. 54 (1983) 4864. [11] D. Farkas, M.J. Mehl, D.A. Papaconstantopoulos, Mater. Res. Soc. Symp. Proc. Z. 653 (2000) 37860. [12] S. Kohlhoff, P. Gumbsch, H.F. Fischmeister, Philos. Mag. A 64 (1991) 851. [13] V. Shastry, D. Farkas, Modelling Simul. Mater. Sci. Eng. 4 (1996) 473. [14] M.W. Finnis, J.E. Sinclair, Philos. Mag. A 50 (1984) 45; M.W. Finnis, J.E. Sinclair, Erratum 53 (1986) 161. [15] G.C. Sih, H. Liebowitz, Mathematical theories of brittle facture, In: Fracture: An Advanced Treatise, vol. 2, Academic Press, New York, 1968 (Chapter 2). [16] A.R. Miedema, Z. Metallkd. 69 (1978) 287. [17] M.J.S. Spencer, A. Hung, I.K. Snook, I. Yarovsky, Surf. Sci. 513 (2002) 389. [18] L. Vitos, A.V. Ruban, H.L. Skriver, J. Kolla´r, Surf. Sci. 411 (1998) 186. [19] P. Błon´ski, A. Kiejna, Vacuum 74 (2004) 179. [20] P. Gumbsch, S.J. Zhou, B.L. Holian, Phys. Rev. B 55 (1997) 3445. [21] W.R. Tyson, R.A. Ayres, D.F. Stein, Acta Metall. 21 (1973) 621.
Atomistic study of lattice trapping behavior for brittle fracture in bcc-iron Ya-Fang Guo a
a,*
, Chong-Yu Wang
b
Institute of Engineering Mechanics, Beijing Jiaotong University, Beijing 100044, China b Department of Physics, Tsinghua University, Beijing 100084, China
Received 31 December 2005; received in revised form 27 November 2006; accepted 29 January 2007 Available online 16 April 2007
Abstract The lattice trapping behavior for brittle fracture in bcc-iron has been studied by atomistic simulations. A pronounced anisotropy for brittle cleavage fracture of a mode I crack is observed in the discrete atomistic scale, and the preferred direction for cleavage is along the h1 1 0i direction on both {1 0 0} and {0 1 1} planes. The analysis of the atomic structure indicates that, due to the shear effect at the crack tip, the stacking faults or partial dislocations are formed before crack cleavage occurs. For the crack with a h1 1 0i crack front, shear occurs easily along the slip direction of bcc crystals, whereas for the crack with a h1 0 0i front, shear can only occur along the non-slip direction, and a strong lattice trapping is exhibited. We conclude that the anisotropy for cleavage fracture and the lattice trapping behavior are closely related with the slip systems of bcc crystals. 2007 Elsevier B.V. All rights reserved. PACS: 71.15.Pd; 62.20.Mk Keywords: Atomistic simulations; Lattice trapping; Anisotropy; Brittle fracture; Slip systems; Bcc-iron
1. Introduction The macroscopic failure of materials is ultimately determined by events occurring on the atomic scale, it is therefore a typical multiscale phenomenon in physics. Traditionally, engineering fracture mechanics, which is based on the continuum mechanics, is used to study material fracture properties, and many useful criterions are obtained. However, it cannot give a fracture criterion from the physical point of view. For example, based on the Griffith criterion for materials brittle fracture [1], it can be predicted that a perfectly brittle crack in a crystal is expected to choose a cleavage plane with low surface energy and to propagate on this plane with equal ease in all directions. But this theoretical prediction shows a clear deviation from
*
Corresponding author. Fax: +86 10 51682094. E-mail address: [email protected] (Y.-F. Guo).
0927-0256/$ - see front matter 2007 Elsevier B.V. All rights reserved. doi:10.1016/j.commatsci.2007.01.019
the experimental observations [2,3] that a crack may not cleavage along the low energy surface preferentially and sometimes a very pronounced anisotropy is observed for cleavage fracture. Moreover, results from atomistic simulations [4] indicate that the anisotropy for cleavage fracture and the lattice trapping behavior are related with atoms bonding, which can not be fully explained by the continuum theory. Lattice trapping effect was proposed about 30 years ago by Thomson et al. [5] when the first atomistic study of fracture was applied using a simple pair-wise potential, and the results indicated that crack propagation was ultimately determined by atoms bonding at the crack tip on the discrete atomic scale. Owing to the lattice trapping, a crack remains stable and not to advance or heal until loads somewhat larger or smaller than the Griffith load are reached. The critical loads for a crack to advance or heal are called the upper and lower trapping limits. Further studies [6,7] showed that the magnitude of the lattice trapping effect
Y.-F. Guo, C.-Y. Wang / Computational Materials Science 40 (2007) 376–381
strongly depends on atoms bonding at the crack tip and the forces needed to separate bonds successively. Subsequently, studies of brittle fracture in Silicon by experiments and simulations [8,9] revealed that the anisotropy for cleavage fracture was strongly related with the lattice trapping effect. Then the ab-initio study of the cleavage fracture [8] indicated that the magnitude of the trapping range decreases with the increasing of the system size. Moreover, the effect of boundary conditions on the lattice trapping behavior for crack propagation in bcc-iron was studied by atomistic simulations [10], and many types of interatomic potentials were used [10–12]. Furthermore, a combined finite element and atomistic simulation method (FEAt) [12] and the molecular static simulation [13] were both used to study the anisotropy for cleavage fracture and the lattice trapping behavior in bcc-iron. From above results, we find that the anisotropy of cleavage fracture and lattice trapping behavior have been widely investigated by atomic scale simulations, and the atoms bonding at the crack tip in particular has been studied. In this paper, the analysis of the slip systems for bcc crystals is introduced. Based on the observation of the atomic structure of the crack tip at the critical load for different types of crack in bcc-iron, we aim to find the relationship between the anisotropy behavior for cleavage fracture and the slip systems for bcc crystals. 2. Simulation model and procedures As shown in Fig. 1, the model used in our study includes two regions: the interior atomistic region I containing the crack tip and the boundary region II. The outer most two layers of atoms in the atomistic region are identified as the boundary region, and its thickness is larger than the cut-off distance of the Finnis–Sinclair N-body potential [14], which has been proved to be efficient for bcc crystal calculations. We consider the plane strain problem and apply the periodic boundary condition in Z direction. The lengths of the atomistic region in X and Y directions ˚ , and the initial crack length is about 37 A ˚. are both 350 A In atomistic simulations we first get the pre-existing crack by shifting the atoms from their positions in a perfect
377
crystal to positions specified by the anisotropic elastic continuum mechanics equations [14] for a certain load. According to anisotropic linear elastic continuum theory, the displacement fields of an open mode I crack for the plane strain condition [15] is given by fug ¼ K I ðrÞ
1=2
fgðhÞg;
ð1Þ
where T
fug ¼ fux ; uy g ;
fgg ¼ fgx ; gy g
T
ð2Þ
and K I ¼ rðapÞ
1=2
ð3Þ
In Eqs. (1)–(3), KI is the stress intensity factor; r and h are the polar coordinates of the atom for which the displacement is calculated; g is a function of the angle h between ~ r and the crack plane; r is the applied load; and a is the length of the crack. In our study we use the stress intensity factor KI to express the corresponding external load r, which can be defined as ‘nKIC’, where KIC is the critical stress intensity factor. In the present work, {0 1 0}h1 0 1i{1 1 0}h1 0 1i,{1 1 0} h0 0 1i and {0 1 0}h0 0 1i mode I (opening mode) cracks in a bcc-iron crystal are studied, where {0 1 0} and {1 1 0} represent the cleavage planes of the crack, h1 0 1i and h0 0 1i are the directions of the crack front. After getting the initial configurations of the crack at different loads and different temperatures, the atoms at the crack tip are fully relaxed for 2000 steps of magnitude 5 · 1015 s with the fixed-displacement boundary condition. The temperature of the system is invariant throughout the simulation, which is achieved by scaling all atoms instantaneous velocities with the appropriate Maxwell–Boltzmann distribution at a specified temperature. By observing atoms bonding at the crack tip, the critical load for crack extension or heal is obtained. 3. Results and discussion 3.1. Calculations of the surface energy Using the F–S N-body potential we firstly calculated the surface energies of the {1 0 0} and {1 1 0} planes for bcciron. The calculated values are: cð100Þ ¼ 1:90 J=m2 cð110Þ ¼ 1:70 J=m
2
and thus cð100Þ > cð110Þ :
Fig. 1. Scheme of the calculation model: region I is the atomistic region, and region II is the border region.
We notice that the calculated surface energy of the {1 0 0} plane is higher than that of the {1 1 0} plane, which is also in good agreement with the previous calculation results by using the same interatomic potential [12]. However, the experimentally determined surface energy of 2.32 J/m2 [16] is rather higher than our calculated values. Moreover, first principle calculations give values of the
378
Y.-F. Guo, C.-Y. Wang / Computational Materials Science 40 (2007) 376–381
surface energy that are close to the experimental value [17– 19], and the {1 0 0} energy value of 2.29 J/m2 is almost identical to the {1 1 0} value of 2.29 J/m2 [17]. For the consistency of the results, we use our calculated values in the present simulation. 3.2. Lattice trapping behavior for crack propagation After getting the initial configurations of the crack at different loads according to anisotropic linear elastic continuum mechanics equations, atoms at the crack tip are fully relaxed with the fixed-displacement boundary condition. Then the critical loads for crack extension at different temperatures are obtained as shown in Table 1. In the previous studies, it has been found that the upper and lower trapping limits are sensitive to the simulation conditions [10,12,20] such as the type of atomistic potential, the boundary condition, the crystallography and the temperature. Table 1 lists our simulation results. It is shown that the calculated upper trapping limits K þ I and the lower trapping limits K I for the {0 1 0}h1 0 1i crack in bcc-iron by using the F–S potential are very close to the theoretical critical stress intensity factor KIC according to the Griffith criterion, while cleavage occurs along the {1 0 0} planes. Comparing with the results obtained by fixed stressed-boundary condition and flexible-boundary condition [10], we find that the displacement-border condition will be valid for simulating crack propagation if the system size is selected sufficiently large. From the results for different types of cracks in Table 1, a strong directional anisotropy for the cleavage of a model I crack is observed. A crack with the h0 0 1i crack front produces pronounced lattice trapping, while a crack with the h0 1 1i front propagates at about the Griffith load. This means that cracks with the h0 0 1i crack front present significantly more lattice trapping than cracks with the h1 1 0i crack front. Moreover, the strongest trapping effects are observed for cleavage along the {1 1 0} planes and it is found that lattice trapping strongly favors cleavage along the {1 0 0} planes. In Table 2 we have compared our results with previous studies. Although the different interatomic potential and atomistic simulation method are used, similar trends are observed that the {1 0 0} plane is the preferred cleavage plane in bcc-iron, and the critical load for a {0 1 0}h1 0 1i crack starting propagation is the lowest compared with that of other cracks. However, there still exist differences between the above calculated values. For examTable 1 The critical load (KI/KIC) for crack propagation at different temperatures Temperature {0 1 0}h1 0 1i {0 1 0}h1 0 1i {1 1 0}h1 1 0i {1 1 0}h0 0 1i {0 1 0}h1 0 0i
K+ K K+ K+ K+
5K
100 K
200 K
300 K
400 K
600 K
1.11 0.89 1.28 1.5 2.00
1.09 0.90 1.25 1.48 1.80
1.07 0.94 1.18 1.46 1.72
1.04 0.98 1.15 1.45 1.71
1.04 0.99 1.13 1.44 1.65
1.05 1.01 1.13 1.40 1.64
Table 2 A comparison with previous calculations Crack type {0 1 0}h1 0 1i {1 1 0}h1 1 0i {1 1 0}h0 0 1i {0 1 0}h1 0 0i
K+ K+ K+ K+
Our work
Ref. [11]
Ref. [12]a
Ref. [12]b
1.11 1.28 1.5 2.00
1.0 1.8 2.4 2.0
1.11 1.72 1.19 1.41
1.07 1.62 1.17 1.75
The critical load for crack propagation is expressed as KI/KIC. a Using the Finnis–Sinclair potential. b Using the Johnson potential.
ple, our results agree well with the results in Ref. [11] that the {1 1 0}h1 1 0i crack can propagate more easily than the {1 1 0}h0 0 1i and {0 1 0}h1 0 0i cracks. Whereas in Ref. [12] the {1 1 0}h0 0 1i crack is predicted to propagate more easily. According to the Griffith theory for brittle fracture, a plane with the lowest surface energy should be the preferred cleavage plane, and a crack will propagate on this plane with equal ease in all directions. In our calculation, a strong anisotropy for crack cleavage exists, and the critical load for a {0 1 0}h1 0 1i crack starting propagation is the lowest compared with that of other cracks, which is also consistent with the experimental results that {1 0 0} is the most likely cleavage plane of iron [21]. However, our calculated surface energy of the {1 0 0} plane is about 12% higher than that of the {1 1 0} plane, and first principle calculations show the {1 0 0} and {1 1 0} surfaces have similar stabilities because the {1 0 0} energy value is almost identical to the {1 1 0} value. This means that the lattice trapping in the discrete atomistic scale causes the anisotropy for cleavage, and the preferred plane for cleavage may not be with the lowest surface energy. Furthermore, we notice that the differences between K þ I , K I and KIC at high temperatures are smaller than those at lower temperatures, which indicates that the trapping effect becomes weak at high temperatures and the trapping barrier can be overcome by thermal activation. It is also shown that the upper and lower trapping limits are almost invariant with the temperature increasing when the temperature is higher than 200 K. As we have known, crack propagation is ultimately determined by atoms bonding at the crack tip, thus we give the relationship between the length of the bond at a crack tip and the load for different types of cracks in Fig. 2. We observe that, when the load is smaller than the upper trapping limits, the lengths of the bond at the crack tip almost have no change with the load increasing for all types of cracks. But if the load is greater than the upper trapping limit, the length of the bond increases rapidly due to the bond breaking, and the crack extends subsequently.
3.3. Lattice trapping effect and slip systems of the crystal Fig. 3 is the critical atomic configuration of a crack tip at 5 K just before the crack tip bond breaking. The solid
Y.-F. Guo, C.-Y. Wang / Computational Materials Science 40 (2007) 376–381
Fig. 2. The relationship between the length of the bond at the crack tip and the load; , indicates the {0 1 0}h1 0 1i crack, indicates the {1 1 0}h1 0 1i crack, n indicates the {1 1 0}h0 0 1i crack, h indicates the {0 1 0}h0 0 1i crack.
and open circles represent atoms in the layer A and B, respectively along the periodic direction. In Fig. 3 we observe that shear occurs at the crack tip at the critical load. For {0 1 0}h1 0 1i and {1 1 0}h1 0 1i cracks in Fig. 3a and b, shear occurs along the h1 1 1i direction, while the stacking fault (SF) is formed at the crack tip. Whereas for {1 1 0}h0 0 1i and {0 1 0}h0 0 1i cracks in Fig. 3c and d,
379
shear occurs along the h0 1 0i direction, which causes the generation of partial dislocations. For bcc crystals, it is well known that the slip planes are {1 1 0}, {2 1 1} or {321}, where the slip direction is always in the h1 1 1i direction. But at low temperatures, the slip system is mainly as the {2 1 1} plane with the h1 1 1i slip direction. As shown in Fig. 4a, the crack fronts for {0 1 0}h1 0 1i and {1 1 0}h1 0 1i cracks are both perpendicular to the {1 0 1} plane (ABCD). Thus the atoms in the {1 0 1} plane can move along the h1 1 1i slip direction, which causes the formation of the SF at the crack tip. The angles a and b between the slip direction and the crack cleavage plane are 35.26 and 54.74, respectively. The slip plane is {2 1 1}. But for {1 1 0} h0 0 1i and {0 1 0}h0 0 1i cracks, the crack fronts are perpendicular to the {0 0 1} plane (ABEF) as shown in Fig. 4b, and the atoms in the {0 0 1} plane cannot move along the h1 1 1i slip direction due to the constraint of the plane strain condition. Because the resistance for atoms movement along the non-slip direction is much greater than that along the slip system in a crystal, thus for the crack with a h0 0 1i front, shear can only occur along the non-slip directions h1 1 0i and h0 1 0i to create partial dislocations when the load is sufficiently large. Therefore, the lattice trapping for the crack with a h0 0 1i front is much stronger than that of the crack with a h1 1 0i front, and the anisotropy for brittle cleavage fracture
Fig. 3. Atomic configuration of the crack tip at the critical load: (a) {0 1 0}h1 0 1i crack, KI = 1.11KIC; (b) {1 1 0}h1 0 1i crack, KI = 1.27KIC; (c) {1 1 0}h0 0 1i crack, KI = 1.49KIC; (d) {0 1 0}h0 0 1i crack, KI = 1.9KIC.
380
Y.-F. Guo, C.-Y. Wang / Computational Materials Science 40 (2007) 376–381
Fig. 4. Lattice structure of a bcc crystal: (a) for the {0 1 0}h1 0 1i and {1 1 0}h1 0 1i cracks; (b) for the {1 1 0}h0 0 1i and {0 1 0}h0 0 1i cracks.
is exhibited, which is closely related with the slip systems of the crystal. Moreover, we notice in Fig. 3 that shear occurs for all types of cracks before crack propagation, which causes the SF formation or partial dislocations creation. We have to mention that the shear would change atoms bonding at the crack tip, and it is closely correlated with lattice trapping behavior. But the precise analysis requires further the first principle calculations.
Fig. 5 gives atoms displacements at a crack tip after relaxation. Open circles represent atoms positions before relaxation, and solid circles represent atoms positions after relaxation. We see clearly that there are more atoms moving away from their initial positions in Fig. 5c and d than those in Fig. 5a and b, and that the magnitudes of atoms displacements in Fig. 5c and d are also greater. Table 3 lists the energy of atomistic region at the critical load for crack propagation. Although the average energy of atoms are
Fig. 5. Atoms displacements at the crack tip after relaxation: (a) {0 1 0}h1 0 1i crack, KI = 1.11KIC; (b) {1 1 0}h1 0 1i crack, KI = 1.27KIC; (c) {1 1 0}h0 0 1i crack, KI = 1.49KIC; and (d) {0 1 0}h0 0 1i crack, KI = 1.9KIC; hollow circles indicate atoms positions before relaxation, and solid circles indicate atoms positions after relaxation.
Y.-F. Guo, C.-Y. Wang / Computational Materials Science 40 (2007) 376–381 Table 3 The energy of atomistic region at the critical load for crack propagation at 5K Crack type
Number of atoms
Eeatom (eV)
Etot (eV)
{0 1 0}h1 0 1i {1 1 0}h1 1 0i {1 1 0}h0 0 1i {0 1 0}h1 0 0i
138,240 138,240 97,200 98,304
4.2606 4.2603 4.2603 4.2545
588,985 588,943 414,101 418,234
almost the same for four types of cracks, the total energy has greater difference between cracks with the h1 0 0i front and h1 1 0i front. For the crack with a h1 0 0i front, higher strain energy is stored at the crack tip before the crack cleavage occurring. Therefore, a higher critical load for the crack propagation is needed for bond breaking at the crack tip, and a strong lattice trapping is found for the crack with a h0 0 1i crack front. 4. Conclusions The lattice trapping behavior and the anisotropy for cleavage fracture are closely related with the slip systems of bcc crystals. For a mode I crack in bcc-iron, the preferred direction for cleavage is along the h1 1 0i direction on both {1 0 0} and {0 1 1} planes, and a stronger lattice trapping is found for the crack with a h1 0 0i crack front than for the crack with a h1 1 0i front. Due to the shear effect at the crack tip, the SF is formed at the crack tip along the slip direction of the bcc crystal for the crack with a h1 1 0i front, just before crack cleavage occurs. While for the crack with a h1 0 0i front, shear can only occur along the non-slip directions and partial dislocations are generated. Meanwhile, high strain energy is found to be stored at the crack tip for the crack with a h0 0 1i front, and a higher critical load for crack propagation is needed for atoms bond breaking, thus a strong lattice trapping is exhibited.
381
Acknowledgements The research was supported by ‘‘973’’ Project from the Ministry of Science and Technology of China (Grant No. 2006CD605102), Chinese Nature Science Foundation (Grant No. 10672016) and NJTU Science Foundation (Grant No. 2005SM0035). References [1] A.A. Griffith, Philos. Trans. Roy. Soc. A 221 (1920) 163. [2] A. George, G. Michot, Mater. Sci. Eng. A 164 (1993) 118. [3] J. Riedle, P. Gumbsch, H.F. Fischmeister, Phys. Rev. Lett. 76 (1996) 3594. [4] P. Gumbsch, Mater. Sci. Eng. A 319–321 (2001) 1. [5] R. Thomson, C. Hsieh, V. Rana, J. Appl. Phys. 42 (1971) 3154. [6] J.E. Sinclair, Philos. Mag. 31 (1975) 647. [7] W.A. Curtin, J. Mater. Res. 5 (1990) 1549. [8] R. Perez, P. Gumbsch, Acta Mater. 48 (2000) 4517. [9] J.C.H. Spence, Y.M. Huang, O. Sankey, Acta Metall. Mater. 41 (1993) 2815. [10] B. DeCelis, A.S. Argon, S. Yip, J. Appl. Phys. 54 (1983) 4864. [11] D. Farkas, M.J. Mehl, D.A. Papaconstantopoulos, Mater. Res. Soc. Symp. Proc. Z. 653 (2000) 37860. [12] S. Kohlhoff, P. Gumbsch, H.F. Fischmeister, Philos. Mag. A 64 (1991) 851. [13] V. Shastry, D. Farkas, Modelling Simul. Mater. Sci. Eng. 4 (1996) 473. [14] M.W. Finnis, J.E. Sinclair, Philos. Mag. A 50 (1984) 45; M.W. Finnis, J.E. Sinclair, Erratum 53 (1986) 161. [15] G.C. Sih, H. Liebowitz, Mathematical theories of brittle facture, In: Fracture: An Advanced Treatise, vol. 2, Academic Press, New York, 1968 (Chapter 2). [16] A.R. Miedema, Z. Metallkd. 69 (1978) 287. [17] M.J.S. Spencer, A. Hung, I.K. Snook, I. Yarovsky, Surf. Sci. 513 (2002) 389. [18] L. Vitos, A.V. Ruban, H.L. Skriver, J. Kolla´r, Surf. Sci. 411 (1998) 186. [19] P. Błon´ski, A. Kiejna, Vacuum 74 (2004) 179. [20] P. Gumbsch, S.J. Zhou, B.L. Holian, Phys. Rev. B 55 (1997) 3445. [21] W.R. Tyson, R.A. Ayres, D.F. Stein, Acta Metall. 21 (1973) 621.