- Email: [email protected]
Combined atomistic simulation and continuum mechanics: Size-dependent behavior of atomistic simulation for brittle fracture in bcc-iron
Computational Materials Science 36 (2006) 432–439 www.elsevier.com/locate/commatsci
Combined atomistic simulation and continuum mechanics: Size-dependent behavior of atomistic simulation for brittle fracture in bcc-iron Ya-Fang Guo *, Yu-Chen Gao Institute of Mechanics, Beijing Jiaotong University, Beijing 100044, China Received 21 December 2004; accepted 18 June 2005
Abstract The size-dependent behavior of atomistic simulation for brittle fracture in bcc-iron is studied using combined continuum-atomistic method. The result of displacement distribution at the crack tip indicates the radius of discrete region at the crack tip is about ˚ , and the truncated distance of continuum elastic field is about 120 A ˚ away from the crack tip. Further investigations of energy 120 A and atomic structure show that the elastic field of continuum mechanics cannot affect the crack tip processes effectively through the ˚ . We assume that the range of the system size that can common region if the diameter of the atomistic region is smaller than 300 A ˚. bridge continuum with discrete region is at the scale of about 300 A 2005 Elsevier B.V. All rights reserved. PACS: 71.15.Pd; 62.20.Mk Keywords: Atomistic simulation; Size-dependent behavior; Brittle fracture; Bcc-iron
1. Introduction The materials fracture is a macroscopic phenomenon, but it is originated from the atomic-scale processes, such as atoms bonds breaking at a crack tip, it is therefore a typical multiscale phenomenon in physics. Traditionally engineering fracture mechanics, just as a methodology, is used to study material fracture properties, and many useful criterions are obtained. In recent years, atomistic simulation has become an effective method to probe the atomic-scale processes of fracture with the development of computer technique, and the combination of contin-
*
Corresponding author. Fax: +86 10 51 68 2094. E-mail address: [email protected] (Y.-F. Guo).
0927-0256/$ - see front matter 2005 Elsevier B.V. All rights reserved. doi:10.1016/j.commatsci.2005.06.014
uum mechanic and atomistic simulation technique has become a powerful tool for the research of material fracture properties [1–4]. Continuum mechanics describes the macroscopically observable behavior of materials by viewing components as ÔcontinuousÕ or ÔhomogenousÕ entities. However, the engineering fracture mechanics, which is based on continuum mechanics, cannot give a fracture criterion from the physical point of view. According to linear elastic continuum theory, the stress field around a crack tip is singular, evidently this is not the real case. However, the far field of a crack tip can indeed be described well by continuum mechanics. Moreover, based on the Griffith criterion for material brittle fracture [5], we can predict 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
Y.-F. Guo, Y.-C. Gao / Computational Materials Science 36 (2006) 432–439
all directions. But this theoretical prediction shows a clear deviation from the experimental observation that a crack may not cleavage along the low energy surface preferentially and sometimes a very pronounced anisotropy is observed for cleavage fracture [6,7]. Based on these, the continuum theory is no longer effective to describe the discrete properties at the crack tip region. In atomistic method, atomistic simulation is used to compute the individual motion of atoms from NewtonÕs equation of motion by using interatomic potential. It is a method of atomistic statistical mechanics, and good for studying the nonlinear or nonequilibrium properties at the crack tip region. By using atomistic simulation, both the atomic structure and dynamic processes of brittle cleavage and dislocation emission at the highly nonlinear region of crack tip are studied [8–13]. Moreover, lattice trapping effect [8,12,14,15] and anisotropic cleavage [14,15] are observed in atomistic simulations, which are caused by the discrete behavior in atomic-scale, and different from the results derived from the continuum analysis. The combined atomistic-continuum technique is particularly a method that spans materials properties from continuum to discrete. In the combined atomistic-continuum method, the boundary condition of the atomic region is usually given by the continuum elastic field of crack according to continuum mechanics. In order to remove the singularity of stress field at a crack tip, very big systems are often used to simulate the crack propagation behaviors, which contains about thousands and millions of atoms [16,17]. However, even the biggest atomistic simulation system that we can obtain today is still too small comparing with the size of real materials. Especially in ab initio study of the crack tip structure, because the system size is limited to tens of atoms, the size effect cannot be ignored [3,18]. Some discrete properties, such as lattice trapping effect, are also found to depend on the size of atomistic regions [19]. Large-scale atomistic simulations suggest that the local hyperelasticity near the crack tip region completely dominates crack dynamics if the size of the hyperelastic region approaches a characteristic length [20], and the sizedependent behavior of atomistic simulation for mode I crack are discussed [21,22]. In this paper, we firstly investigate the displacement distributions at a crack tip both with continuum and discrete method, and the truncated distance of continuum elastic field is discussed. Furthermore, the brittle cleavage fracture in bcc-iron crystal at low temperatures is studied with different sizes of system by using combined continuum-atomistic method. The minimum size of atomistic region for the study of brittle fracture processes is discussed, and we try to find the range of the system size that can bridge linear elastic continuum region with nonlinear discrete region.
433
2. Simulation model and procedures 2.1. Simulation model of crack As shown in Fig. 1, the model used in our study includes three regions: the interior atomistic region I containing the crack tip, the exterior continuum region II, and the boundary region III. The boundary region is a common zone between the atomistic region and the continuum region, and the continuum region can be taken into account in atomistic simulation through its effect on the common boundary region. In our simulation, {0 1 0}h1 0 0i and {0 1 0}h1 0 1i mode I (opening mode) crack in bcc-iron crystal are studied. The tensile stress is applied in h0 1 0i direction, which is perpendicular to the crack surface of {0 1 0} plane, and the crack front is oriented along h1 0 0i and h1 0 1i direction, respectively. We consider the plane strain problem, and apply the periodic boundary condition to crack front direction. 2.2. Continuum theory for crack tip According to anisotropic linear elastic continuum theory, the displacement field u and the stress field r of a straight crack for plane strain condition [23] are given as fug ¼ K I ðrÞ1=2 fgðhÞg
ð1Þ
and frg ¼
KI ðrÞ
1=2
ff ðhÞg
ð2Þ
Fig. 1. Model of crack: region I is the atomistic regions, region II is the continuum region, and region III is the common zone between atomistic and continuum.
434
Y.-F. Guo, Y.-C. Gao / Computational Materials Science 36 (2006) 432–439
where
3. Results and discussion T
frg ¼ frx ; ry ; rxy g ; T
fug ¼ fux ; uy g ;
T
ff g ¼ ffx ; fy ; fxy g ;
fgg ¼ fgx ; gy g
T
In Eqs. (1) and (2), KI is the stress intensity factor, and r labels the distance from the crack tip. f and g are functions of the angle h between ~ r and the crack plane, including the crystallographic orientation of crack system via the appropriate elastic constants [23]. In our study we use the stress intensity factor KI to express the corresponding external load r, it can be defined as ÔaKICÕ, where KIC is the critical stress intensity factor, a is the coefficient that is greater than zero. Considering the anisotropic of discrete lattice [23], KIC can be given as h p n 1=2 1=2 K IC ¼ 2c ðs11 s22 =2Þ ðs22 =s11 Þ i1=2 1=2 þð2s12 þ s66 Þ=2s11 ð3Þ where c is the free surface energy, and sij are the compliance coefficient corrected for plane strain conditions. 2.3. Simulation condition and procedure Before the atomistic simulation of crack propagation, we first get the starting configuration of atoms at the crack tip according to the elastic displacement field calculated from Eq. (1), where r and h are the polar coordinates of the atom for which the displacement is calculated. In order to study the size-dependent behavior of atomistic simulation, different sizes of atomistic system are taken, which are expressed as ÔLX*LYÕ, where LX and LY indicate the length of atomistic region separately in X- and Y-direction. In our simulation LX and LY always keep equal, thus we use LY to express the size of the system. In the periodic direction (Z-direction) six atomic layers are contained, and the length of initial ˚ . The Finnis– crack in atomistic region is about 40.5 A Sinclair N-body potential [24] is applied in atomistic simulation, which has been proved to be efficient for bcc crystal calculations. After the starting atomic configuration of a crack at a specified stress intensity KI is obtained, the atoms are fully relaxed under the fixed-displacement boundary condition. The outer most two layers of atoms in the atomistic region are identified as common region between atomistic and continuum region, and the thickness of the common region is bigger than the cut-off distance of the Finnis–Sinclair N-body potential. The temperature of the system has been kept constant at 5 K during the simulation, which is obtained by scaling all atoms instantaneous velocities with the appropriate Maxwell–Boltzmann distribution at the specified temperature. The magnitude of time steps is 5 · 1015 s.
3.1. Distributions of displacement at a crack tip In this section, we select {0 1 0}h1 0 0i crack to study the displacement distribution at a crack tip. In order to get the critical stress intensity factor KIC, the free surface energy for bcc-iron is calculated using the Finnis–Sinclair N-body potential. The calculated surface energy of (0 1 0) surface is 1.90 Jm2. The compliance coefficient s11 and s22 are 6.99 pPa1, s12 is 2.53 pPa1, and s66 is 8.20 pPa1. Thus the critical stress intensity factor KIC for {0 1 0}h1 0 0i crack we obtained 0.8334 MN m3/2. From the continuum Eq. (1) we get the starting atomic configuration of a {0 1 0}h1 0 0i crack at the critical stress intensity factor KIC, and it is shown in Fig. 2 ˚ simulation system is selected in as solid circles. 710 A ˚ . Subsequently this study where LX and LY are 710 A the atoms are fully relaxed for 4000 time steps under the fixed-displacement boundary condition, which means that the displacement of atoms in the common region III are fixed and the atoms in region I are fully relaxed. The relaxed atomic configuration is shown in Fig. 2 as open circles. We can see clearly in Fig. 2 that the atoms at the near crack tip region move away from their initial position after relaxation, such as atoms A, C, E and F. Whereas for atoms located far away from the crack tip, its movement cannot be observed obviously.
Fig. 2. Atomic configuration of crack at KI = KIC. Solid circles represent unrelaxed atoms, and open circles represent relaxed atoms.
Y.-F. Guo, Y.-C. Gao / Computational Materials Science 36 (2006) 432–439
435
Fig. 3. Influence of relaxation on atoms displacement for a {0 1 0}h1 0 0i crack. In (a) and (b) solid squares and open squares represent displacements for unrelaxed atoms and relaxed atoms, respectively; dY and dX represent the change of atoms displacement after relaxation in Y- and X-direction.
Fig. 3 gives atom displacement before and after relaxation. Because tensile loading is applied in Y-direction, we first discuss the displacement of an array of atoms that is perpendicular to crack surface; r represents the distance between the atom and the crack tip. UY and UX represent atoms displacement in Y- and X-direction, respectively, which are obtained by comparing atom position of a crack with its original position in a perfect crystal. In Fig. 3(a) and (b) we can find that UY and UX increase with r increasing, but the displacement of relaxed atoms is larger than that of unrelaxed atoms. The difference of displacement between relaxed and unrelaxed atoms is shown in Fig. 3(c) and (d), and we can see clearly that the absolute value of dy decreases ˚ rapidly with r increasing when r is smaller than 120 A ˚ or so; while r is greater than 120 A, dY almost changes no more, and approaches zero. These results indicate that the atoms relaxation process can only affect atoms displacement significantly within the crack tip region
˚ . Furthermore, because tenat the radius of about 120 A sile stress is applied in Y-direction, dX almost keeps zero except for a few atoms just at the near crack tip region. Because we get the displacement of unrelaxed atoms according to continuum elasticity theory, and the relaxed results calculated from discrete method, it can be concluded that the radius of the discrete nonlinear re˚ , and the truncated gion at the crack tip is about 120 A ˚ away distance of continuum elastic field is about 120 A from the crack tip. This result is also consistent with the previous study from the energy flux calculation that the hyperelastic zone near the crack tip is on the order of a few hundred atomic spacing [20]. Subsequently we test different sizes of system and investigate the influence of system size on the displacement field at a crack tip. In Fig. 4 we compare the results ˚ system with 423 A ˚ and 250 A ˚ system, and find of 710 A that the size of system cannot affect the atoms relaxation displacement significantly if the radius of atomistic
436
Y.-F. Guo, Y.-C. Gao / Computational Materials Science 36 (2006) 432–439
˚ system, solid squares in (a) indicate 423 A ˚ system, Fig. 4. Atoms displacement for different system sizes. Open squares in (a) and (b) indicate 710 A ˚ system. and solid squares in (b) indicate 250 A
region is greater than the truncated distance of continuum elastic field. Moreover, we calculate the displacement field of a {0 1 0}h1 0 1i crack, and results are shown in Fig. 5, which are exactly similar to the results for a {0 1 0}h1 0 0i crack. It indicates that the radius of the discrete region at the crack tip region is not dependent on the type of cracks, and the truncated distance ˚ away from of continuum elastic field is also about 120 A the crack tip. 3.2. Cleavage fracture at low temperature In this section, we focus on the study of the sizedependent behavior of crack brittle cleavage fracture by using different sizes of simulation system. Considering the lattice trapping effect [13], {0 1 0}h1 0 1i mode I cracks are selected in this section for studying brittle cleavage behavior in bcc-iron. The length of atomistic
region in X- and Y-direction (LX and LY) always keeps ˚ , just as that equal when LX and LY is greater than 140 A in Section 3.1. We must mention that, if the size of the system is too small, the crack tip may reach the boundary of atomistic region after the crack has moved along the cleavage plane. Thus in the present model, when the ˚ , we will keep LX size of system is smaller than 140 A ˚ invariant at the value of 140 A and only LY decreases with system size decreasing, and LY is used to express the size of system. The total number of atoms in atomistic region is therefore ranged from 3000 to 378,000, ˚ to 602 A ˚. corresponding LY changes from 28 A Crack propagations at a static relaxation condition are simulated firstly in this section. After we get the starting crack at different loads separately, the atoms of crack are fully relaxed for 4000 time steps while the outer most two layers of atoms in atomistic region are fixed. The simulation result indicates that the crack
Fig. 5. Influence of relaxation on atoms displacement for a {0 1 0}h1 0 1i crack. In (a) solid squares and open squares represent displacements for unrelaxed atoms and relaxed atoms, respectively; and in (b) dY represent the change of atoms displacement after relaxation in Y-direction.
Y.-F. Guo, Y.-C. Gao / Computational Materials Science 36 (2006) 432–439
437
˚ . As shown in Fig. 7, there also exis greater than 350 A ists an inflexion in the curve of the length of crack, just as the curve of atoms average energy in Fig. 6. When LY ˚ , the final length of the crack inis smaller than 250 A creases greatly with LY increasing. While LY is greater ˚ or so, the increasing of the length of crack than 300 A become more and more slowly. This indicates that the change of the size of atomistic region cannot influence the simulation result of cleavage fracture evidently if ˚. the size of atomistic region is greater than 300 A 3.3. Evolution of atomistic structure at the crack tip
Fig. 6. Atoms average energy for different system sizes at KI = 1.6KIC.
moves along cleavage plane and cleavage fracture occurs when the load is lower than 1.8KIC. Fig. 6 gives atoms average energy in atomistic region for different system sizes at 1.6KIC, and Fig. 7 shows the final length of the crack after 4000 steps relaxation. In Fig. 6 we find that the average energy of atoms decreases greatly with LY increasing when LY is smaller ˚ . As LY is greater than 200 A ˚ , the decrease than 200 A of atoms average energy become more and more slow, and the change has become very small when LY
Fig. 7. Length of crack for different system sizes at KI = 1.6KIC. Lc indicates the final length of crack after 4000 steps relaxation.
In our previous work, we studied that twin formed at the crack tip accompanying crack cleavage in brittle fracture at high loading rates or at low temperatures [13]. In this section we have tested different sizes of system, and try to find the influence of the size of atomistic region on twin formation in brittle cleavage fracture. The system size is just the same as which in Section 3.2. Before simulation, we at first get a starting crack at the critical stress intensity factor KIC, and atoms in the initial crack are fully relaxed. Subsequently the external loading is applied increasingly by re-scaling the displacement of atoms at the boundary region according to linear elasticity theory. For each size of system, the load increases incrementally from 1.0KIC to 2.0KIC, and the load step is 0.1KIC. At every loading process, the system is relaxed for 300 steps, thus the loading rate is 0.1KIC/300 steps. The atomic structure of the crack tip at different loads for different system sizes are analyzed below: when ˚ (about 30 layers of atoms in LY is smaller than 45 A [0 1 0] direction), crack propagation along {0 1 0} cleavage plane cannot be observed. Fig. 8(a)–(c) show the atomic configuration of the crack tip for the system that 20 layers of atoms are contained in [0 1 0] direction. In Fig. 8(a) we can find the breaking of bonds at the crack tip at 1.2KIC. But with the load increasing, the crack tip does not move any more. When the load increases to 1.6KIC, the crack tip is blunted as shown in Fig. 8(b), and the crack cleavage along [0 1 0] direction, which is perpendicular to crack surface. ˚ , crack propagation When LY is greater than 45 A along cleavage plane is observed. Fig. 8(d)–(i) show the atomic configurations of crack tips while 40 and 98 atomic layers are contained along [0 1 0] direction. As shown in Fig. 8(d)–(f), the breaking of bond at the crack tip is observed at 1.2KIC, and the crack tip moves continuously along {0 1 0} cleavage plane with the load increasing. Meanwhile, stacking fault is found to form at the crack tip while crack tip moves, but no twin strip observed. When the number of atomic layer is greater ˚ ), the crack begins to cleavage at than 98 (about 140 A 1.1KIC. As shown in Fig. 8(h) and (i), twin formation is observed at the crack tip accompanying crack cleavage
438
Y.-F. Guo, Y.-C. Gao / Computational Materials Science 36 (2006) 432–439
˚ , KI = 1.2KIC; (b) LY = 29 A ˚ , KI = 1.6KIC; (c) Fig. 8. Atomic configuration of crack at different loads for different system sizes: (a) LY = 29 A ˚ , KI = 1.9KIC; (d) LY = 57 A ˚ , KI = 1.2KIC; (e) LY = 57 A ˚ , KI = 1.6KIC; (f) LY = 57 A ˚ , KI = 1.9KIC; (g) LY = 140 A ˚ , KI = 1.1KIC; (h) LY = 29 A ˚ , KI = 1.6KIC; (i) LY = 140 A ˚ , KI = 1.9KIC. LY indicates the length of atomistic region in Y-direction. NZ = 140 A
propagation, which is agreed with the results that we obtained in our previous studies [13]. Moreover, we also compared the atomic configuration of crack tip at 1.9 KIC while LY is, respectively, at ˚ , 280 A ˚ and 360 A ˚ . We find, for the 140 A ˚ system, 140 A both the width and length of twin are a little smaller ˚ and 360 A ˚ systems. While LY is than those for 280 A ˚ , the width and length of twin at the greater than 280 A crack tip change no more. Additionally, it is also found that the length of the crack do not change anymore ˚ . It indicates that both when LY is greater than 280 A the atomic configuration and the length of crack do not change distinctly with the system size increasing if ˚ or so. Thus we can conclude LY is greater than 280 A that in order to reflect the real process of brittle fracture, the diameter of atomistic region cannot be smaller than ˚. 300 A
4. Conclusions In this paper, we first investigate the displacement distributions at a crack tip both with continuum and discrete method. The results indicate that the proper radius ˚ , that of the discrete region at a crack tip is about 120 A is also the truncated distance of continuum elastic field around a crack tip. Subsequently, the influence of the
size of atomistic region on brittle fracture for bcc-iron crystal is studied by using combined continuum-atomistic method. The atomistic model is used for the region near the crack tip while the continuum model is used for the outside domain, and the solution of the displacement field for opening crack according to anisotropic continuum fracture mechanics provides the boundary condition of the atomistic region. The results indicate that the necessary minimum size of atomistic model ˚ . When the size of for brittle fracture is about 300 A atomistic region is smaller than this value, the elastic field of continuum mechanics cannot reflect the atomic processes at the crack tip properly. Thus we assume the system size that can bridge continuum with discrete region is about 30 nm. The elastic continuum mechanics cannot describe materials brittle fracture behavior effectively when the system size is smaller than this value, and the discrete properties should be considered while the material behavior at nanometer scale is concerned.
Acknowledgments The research was supported by Chinese Nature Science Foundation (No. 10228204), NJTU Science Foundation (No. 2004RC006) and NJTU Paper
Y.-F. Guo, Y.-C. Gao / Computational Materials Science 36 (2006) 432–439
Foundation. The authors thank Prof. Yuesheng Wang and Dr. Suowen Gao for their helpful discussions.
References [1] P. Gumbsch, An atomistic study of brittle fracture: toward explicit failure criteria from atomistic modeling, J. Mater. Res. 10 (1995) 2897. [2] H. Rafii-Tabar, L. Hua, M. Cross, A multi-scale atomisticcontinuum modelling of crack propagation in a two-dimensional macroscopic plate, J. Phys. 10 (1998) 2375. [3] J.Q. Broughton, F.F. Abraham, N. Bernstein, E. Kaxiars, Concurrent coupling of length scales: methodology and application, Phys. Rev. B 60 (1999) 2391. [4] S.V. Bhenoy, R. Miller, E.B. Tadmor, D. Rodney, R. Phillips, M. Ortiz, An adaptive finite element approach to atomic-scale mechanics—the quasicontinuum method, J. Mech. Phys. Solids 47 (1999) 611. [5] A.A. Griffith, The phenomena of rupture and flow in solids, Philos. Trans. Roy. Soc. A 221 (1920) 163. [6] A. George, G. Michot, Dislocation loops at crack tips: nucleation and growth—an experimental study in silicon, Mater. Sci. Eng. A 164 (1993) 118. [7] J. Riedle, P. Gumbsch, H.F. Fischmeister, Cleavage anisotropy in tungsten single crystals, Phys. Rev. Lett. 76 (1996) 3594. [8] B. DeCelis, A.S. Argon, S. Yip, Molecular dynamics simulation of crack tip processes in alpha-iron and copper, J. Appl. Phys. 54 (1983) 4864. [9] K.S. Cheng, S. Yip, Brittle–ductile transition in intrinsic fracture behavior of crystals, Phys. Rev. Lett. 65 (1990) 2804. [10] F.F. Abraham, Dynamics of brittle fracture with variable elasticity, Phys. Rev. Lett. 77 (1996) 869. [11] F. Cleri, S. Yip, D. Wolf, S.R. Phillpot, Atomistic-scale mechanism of crack-tip plasticity: dislocation nucleation and crack tip shielding, Phys. Rev. Lett. 79 (1997) 1309.
439
[12] S. Kohlhoff, P. Gumbsch, H.F. Fischmeister, Crack propagation in bcc crystals studied with a combined finite-element and atomistic model, Philos. Mag. A 64 (1991) 851. [13] Y.F. Guo, C.Y. Wang, D.L. Zhao, Atomistic simulation of crack cleavage and blunting in bcc-Fe, Mater. Sci. Eng. A 349 (2003) 29. [14] P. Gumbsch, Modelling brittle and semi-brittle fracture processes, Mater. Sci. Eng. A 319–321 (2001) 1. [15] P. Gumbsch, R.M. Cannon, Atomistic aspects of brittle fracture, MRS Bull. 25 (2000) 15. [16] F.F. Abraham, D. Schneider, B. Land, D. Lifka, J. Skovira, J. Gerner, M. Rosenkrantz, Instability dynamics in three-dimensional fracture: an atomistic simulation, J. Mech. Phys. Solids 45 (1997) 1461. [17] V. Bulatov, F.F. Abraham, L. Kubin, B. Devincre, S. Yip, Connecting atomistic and mesoscale simulations of crystal plasticity, Nature 391 (1998) 669. [18] R. Perez, P. Gumbsch, An ab initio study of the cleavage anisotropy in silicon, Acta Mater. 48 (2000) 4517. [19] J.C.H. Spence, Y.M. Huang, O. Sankey, Lattice trapping and surface reconstruction for silicon cleavage on (1 1 1). Ab initio quantum molecular dynamics calculation, Acta Metall. Mater. 41 (1993) 2815. [20] M.J. Buehler, F.F. Abraham, H. Gao, Hyperelasticity governs dynamic fracture at a critical length scale, Nature 426 (2003) 141. [21] P. Espan˜ol, I. Zu´n˜iga, M.A. Rubio, Effect of boundary conditions in mode I fracture in brittle materials, Physica D 96 (1996) 375. [22] B.L. Holian, R. Revelo, Fracture simulations using large-scale molecular dynamics, Phys. Rev. B 51 (1995) 11275. [23] G.C. Sih, H. Liebowitz, Mathematical theories of brittle facture, Fracture: An Advanced Treatise, vol. 2, Academic, New York, 1968 (Chapter 2). [24] M.W. Finnis, J.E. Sinclair, A simple N-body potential for transition metals, Philos. Mag. A 50 (1984) 45; Erratum 53 (1986) 161.
Combined atomistic simulation and continuum mechanics: Size-dependent behavior of atomistic simulation for brittle fracture in bcc-iron Ya-Fang Guo *, Yu-Chen Gao Institute of Mechanics, Beijing Jiaotong University, Beijing 100044, China Received 21 December 2004; accepted 18 June 2005
Abstract The size-dependent behavior of atomistic simulation for brittle fracture in bcc-iron is studied using combined continuum-atomistic method. The result of displacement distribution at the crack tip indicates the radius of discrete region at the crack tip is about ˚ , and the truncated distance of continuum elastic field is about 120 A ˚ away from the crack tip. Further investigations of energy 120 A and atomic structure show that the elastic field of continuum mechanics cannot affect the crack tip processes effectively through the ˚ . We assume that the range of the system size that can common region if the diameter of the atomistic region is smaller than 300 A ˚. bridge continuum with discrete region is at the scale of about 300 A 2005 Elsevier B.V. All rights reserved. PACS: 71.15.Pd; 62.20.Mk Keywords: Atomistic simulation; Size-dependent behavior; Brittle fracture; Bcc-iron
1. Introduction The materials fracture is a macroscopic phenomenon, but it is originated from the atomic-scale processes, such as atoms bonds breaking at a crack tip, it is therefore a typical multiscale phenomenon in physics. Traditionally engineering fracture mechanics, just as a methodology, is used to study material fracture properties, and many useful criterions are obtained. In recent years, atomistic simulation has become an effective method to probe the atomic-scale processes of fracture with the development of computer technique, and the combination of contin-
*
Corresponding author. Fax: +86 10 51 68 2094. E-mail address: [email protected] (Y.-F. Guo).
0927-0256/$ - see front matter 2005 Elsevier B.V. All rights reserved. doi:10.1016/j.commatsci.2005.06.014
uum mechanic and atomistic simulation technique has become a powerful tool for the research of material fracture properties [1–4]. Continuum mechanics describes the macroscopically observable behavior of materials by viewing components as ÔcontinuousÕ or ÔhomogenousÕ entities. However, the engineering fracture mechanics, which is based on continuum mechanics, cannot give a fracture criterion from the physical point of view. According to linear elastic continuum theory, the stress field around a crack tip is singular, evidently this is not the real case. However, the far field of a crack tip can indeed be described well by continuum mechanics. Moreover, based on the Griffith criterion for material brittle fracture [5], we can predict 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
Y.-F. Guo, Y.-C. Gao / Computational Materials Science 36 (2006) 432–439
all directions. But this theoretical prediction shows a clear deviation from the experimental observation that a crack may not cleavage along the low energy surface preferentially and sometimes a very pronounced anisotropy is observed for cleavage fracture [6,7]. Based on these, the continuum theory is no longer effective to describe the discrete properties at the crack tip region. In atomistic method, atomistic simulation is used to compute the individual motion of atoms from NewtonÕs equation of motion by using interatomic potential. It is a method of atomistic statistical mechanics, and good for studying the nonlinear or nonequilibrium properties at the crack tip region. By using atomistic simulation, both the atomic structure and dynamic processes of brittle cleavage and dislocation emission at the highly nonlinear region of crack tip are studied [8–13]. Moreover, lattice trapping effect [8,12,14,15] and anisotropic cleavage [14,15] are observed in atomistic simulations, which are caused by the discrete behavior in atomic-scale, and different from the results derived from the continuum analysis. The combined atomistic-continuum technique is particularly a method that spans materials properties from continuum to discrete. In the combined atomistic-continuum method, the boundary condition of the atomic region is usually given by the continuum elastic field of crack according to continuum mechanics. In order to remove the singularity of stress field at a crack tip, very big systems are often used to simulate the crack propagation behaviors, which contains about thousands and millions of atoms [16,17]. However, even the biggest atomistic simulation system that we can obtain today is still too small comparing with the size of real materials. Especially in ab initio study of the crack tip structure, because the system size is limited to tens of atoms, the size effect cannot be ignored [3,18]. Some discrete properties, such as lattice trapping effect, are also found to depend on the size of atomistic regions [19]. Large-scale atomistic simulations suggest that the local hyperelasticity near the crack tip region completely dominates crack dynamics if the size of the hyperelastic region approaches a characteristic length [20], and the sizedependent behavior of atomistic simulation for mode I crack are discussed [21,22]. In this paper, we firstly investigate the displacement distributions at a crack tip both with continuum and discrete method, and the truncated distance of continuum elastic field is discussed. Furthermore, the brittle cleavage fracture in bcc-iron crystal at low temperatures is studied with different sizes of system by using combined continuum-atomistic method. The minimum size of atomistic region for the study of brittle fracture processes is discussed, and we try to find the range of the system size that can bridge linear elastic continuum region with nonlinear discrete region.
433
2. Simulation model and procedures 2.1. Simulation model of crack As shown in Fig. 1, the model used in our study includes three regions: the interior atomistic region I containing the crack tip, the exterior continuum region II, and the boundary region III. The boundary region is a common zone between the atomistic region and the continuum region, and the continuum region can be taken into account in atomistic simulation through its effect on the common boundary region. In our simulation, {0 1 0}h1 0 0i and {0 1 0}h1 0 1i mode I (opening mode) crack in bcc-iron crystal are studied. The tensile stress is applied in h0 1 0i direction, which is perpendicular to the crack surface of {0 1 0} plane, and the crack front is oriented along h1 0 0i and h1 0 1i direction, respectively. We consider the plane strain problem, and apply the periodic boundary condition to crack front direction. 2.2. Continuum theory for crack tip According to anisotropic linear elastic continuum theory, the displacement field u and the stress field r of a straight crack for plane strain condition [23] are given as fug ¼ K I ðrÞ1=2 fgðhÞg
ð1Þ
and frg ¼
KI ðrÞ
1=2
ff ðhÞg
ð2Þ
Fig. 1. Model of crack: region I is the atomistic regions, region II is the continuum region, and region III is the common zone between atomistic and continuum.
434
Y.-F. Guo, Y.-C. Gao / Computational Materials Science 36 (2006) 432–439
where
3. Results and discussion T
frg ¼ frx ; ry ; rxy g ; T
fug ¼ fux ; uy g ;
T
ff g ¼ ffx ; fy ; fxy g ;
fgg ¼ fgx ; gy g
T
In Eqs. (1) and (2), KI is the stress intensity factor, and r labels the distance from the crack tip. f and g are functions of the angle h between ~ r and the crack plane, including the crystallographic orientation of crack system via the appropriate elastic constants [23]. In our study we use the stress intensity factor KI to express the corresponding external load r, it can be defined as ÔaKICÕ, where KIC is the critical stress intensity factor, a is the coefficient that is greater than zero. Considering the anisotropic of discrete lattice [23], KIC can be given as h p n 1=2 1=2 K IC ¼ 2c ðs11 s22 =2Þ ðs22 =s11 Þ i1=2 1=2 þð2s12 þ s66 Þ=2s11 ð3Þ where c is the free surface energy, and sij are the compliance coefficient corrected for plane strain conditions. 2.3. Simulation condition and procedure Before the atomistic simulation of crack propagation, we first get the starting configuration of atoms at the crack tip according to the elastic displacement field calculated from Eq. (1), where r and h are the polar coordinates of the atom for which the displacement is calculated. In order to study the size-dependent behavior of atomistic simulation, different sizes of atomistic system are taken, which are expressed as ÔLX*LYÕ, where LX and LY indicate the length of atomistic region separately in X- and Y-direction. In our simulation LX and LY always keep equal, thus we use LY to express the size of the system. In the periodic direction (Z-direction) six atomic layers are contained, and the length of initial ˚ . The Finnis– crack in atomistic region is about 40.5 A Sinclair N-body potential [24] is applied in atomistic simulation, which has been proved to be efficient for bcc crystal calculations. After the starting atomic configuration of a crack at a specified stress intensity KI is obtained, the atoms are fully relaxed under the fixed-displacement boundary condition. The outer most two layers of atoms in the atomistic region are identified as common region between atomistic and continuum region, and the thickness of the common region is bigger than the cut-off distance of the Finnis–Sinclair N-body potential. The temperature of the system has been kept constant at 5 K during the simulation, which is obtained by scaling all atoms instantaneous velocities with the appropriate Maxwell–Boltzmann distribution at the specified temperature. The magnitude of time steps is 5 · 1015 s.
3.1. Distributions of displacement at a crack tip In this section, we select {0 1 0}h1 0 0i crack to study the displacement distribution at a crack tip. In order to get the critical stress intensity factor KIC, the free surface energy for bcc-iron is calculated using the Finnis–Sinclair N-body potential. The calculated surface energy of (0 1 0) surface is 1.90 Jm2. The compliance coefficient s11 and s22 are 6.99 pPa1, s12 is 2.53 pPa1, and s66 is 8.20 pPa1. Thus the critical stress intensity factor KIC for {0 1 0}h1 0 0i crack we obtained 0.8334 MN m3/2. From the continuum Eq. (1) we get the starting atomic configuration of a {0 1 0}h1 0 0i crack at the critical stress intensity factor KIC, and it is shown in Fig. 2 ˚ simulation system is selected in as solid circles. 710 A ˚ . Subsequently this study where LX and LY are 710 A the atoms are fully relaxed for 4000 time steps under the fixed-displacement boundary condition, which means that the displacement of atoms in the common region III are fixed and the atoms in region I are fully relaxed. The relaxed atomic configuration is shown in Fig. 2 as open circles. We can see clearly in Fig. 2 that the atoms at the near crack tip region move away from their initial position after relaxation, such as atoms A, C, E and F. Whereas for atoms located far away from the crack tip, its movement cannot be observed obviously.
Fig. 2. Atomic configuration of crack at KI = KIC. Solid circles represent unrelaxed atoms, and open circles represent relaxed atoms.
Y.-F. Guo, Y.-C. Gao / Computational Materials Science 36 (2006) 432–439
435
Fig. 3. Influence of relaxation on atoms displacement for a {0 1 0}h1 0 0i crack. In (a) and (b) solid squares and open squares represent displacements for unrelaxed atoms and relaxed atoms, respectively; dY and dX represent the change of atoms displacement after relaxation in Y- and X-direction.
Fig. 3 gives atom displacement before and after relaxation. Because tensile loading is applied in Y-direction, we first discuss the displacement of an array of atoms that is perpendicular to crack surface; r represents the distance between the atom and the crack tip. UY and UX represent atoms displacement in Y- and X-direction, respectively, which are obtained by comparing atom position of a crack with its original position in a perfect crystal. In Fig. 3(a) and (b) we can find that UY and UX increase with r increasing, but the displacement of relaxed atoms is larger than that of unrelaxed atoms. The difference of displacement between relaxed and unrelaxed atoms is shown in Fig. 3(c) and (d), and we can see clearly that the absolute value of dy decreases ˚ rapidly with r increasing when r is smaller than 120 A ˚ or so; while r is greater than 120 A, dY almost changes no more, and approaches zero. These results indicate that the atoms relaxation process can only affect atoms displacement significantly within the crack tip region
˚ . Furthermore, because tenat the radius of about 120 A sile stress is applied in Y-direction, dX almost keeps zero except for a few atoms just at the near crack tip region. Because we get the displacement of unrelaxed atoms according to continuum elasticity theory, and the relaxed results calculated from discrete method, it can be concluded that the radius of the discrete nonlinear re˚ , and the truncated gion at the crack tip is about 120 A ˚ away distance of continuum elastic field is about 120 A from the crack tip. This result is also consistent with the previous study from the energy flux calculation that the hyperelastic zone near the crack tip is on the order of a few hundred atomic spacing [20]. Subsequently we test different sizes of system and investigate the influence of system size on the displacement field at a crack tip. In Fig. 4 we compare the results ˚ system with 423 A ˚ and 250 A ˚ system, and find of 710 A that the size of system cannot affect the atoms relaxation displacement significantly if the radius of atomistic
436
Y.-F. Guo, Y.-C. Gao / Computational Materials Science 36 (2006) 432–439
˚ system, solid squares in (a) indicate 423 A ˚ system, Fig. 4. Atoms displacement for different system sizes. Open squares in (a) and (b) indicate 710 A ˚ system. and solid squares in (b) indicate 250 A
region is greater than the truncated distance of continuum elastic field. Moreover, we calculate the displacement field of a {0 1 0}h1 0 1i crack, and results are shown in Fig. 5, which are exactly similar to the results for a {0 1 0}h1 0 0i crack. It indicates that the radius of the discrete region at the crack tip region is not dependent on the type of cracks, and the truncated distance ˚ away from of continuum elastic field is also about 120 A the crack tip. 3.2. Cleavage fracture at low temperature In this section, we focus on the study of the sizedependent behavior of crack brittle cleavage fracture by using different sizes of simulation system. Considering the lattice trapping effect [13], {0 1 0}h1 0 1i mode I cracks are selected in this section for studying brittle cleavage behavior in bcc-iron. The length of atomistic
region in X- and Y-direction (LX and LY) always keeps ˚ , just as that equal when LX and LY is greater than 140 A in Section 3.1. We must mention that, if the size of the system is too small, the crack tip may reach the boundary of atomistic region after the crack has moved along the cleavage plane. Thus in the present model, when the ˚ , we will keep LX size of system is smaller than 140 A ˚ invariant at the value of 140 A and only LY decreases with system size decreasing, and LY is used to express the size of system. The total number of atoms in atomistic region is therefore ranged from 3000 to 378,000, ˚ to 602 A ˚. corresponding LY changes from 28 A Crack propagations at a static relaxation condition are simulated firstly in this section. After we get the starting crack at different loads separately, the atoms of crack are fully relaxed for 4000 time steps while the outer most two layers of atoms in atomistic region are fixed. The simulation result indicates that the crack
Fig. 5. Influence of relaxation on atoms displacement for a {0 1 0}h1 0 1i crack. In (a) solid squares and open squares represent displacements for unrelaxed atoms and relaxed atoms, respectively; and in (b) dY represent the change of atoms displacement after relaxation in Y-direction.
Y.-F. Guo, Y.-C. Gao / Computational Materials Science 36 (2006) 432–439
437
˚ . As shown in Fig. 7, there also exis greater than 350 A ists an inflexion in the curve of the length of crack, just as the curve of atoms average energy in Fig. 6. When LY ˚ , the final length of the crack inis smaller than 250 A creases greatly with LY increasing. While LY is greater ˚ or so, the increasing of the length of crack than 300 A become more and more slowly. This indicates that the change of the size of atomistic region cannot influence the simulation result of cleavage fracture evidently if ˚. the size of atomistic region is greater than 300 A 3.3. Evolution of atomistic structure at the crack tip
Fig. 6. Atoms average energy for different system sizes at KI = 1.6KIC.
moves along cleavage plane and cleavage fracture occurs when the load is lower than 1.8KIC. Fig. 6 gives atoms average energy in atomistic region for different system sizes at 1.6KIC, and Fig. 7 shows the final length of the crack after 4000 steps relaxation. In Fig. 6 we find that the average energy of atoms decreases greatly with LY increasing when LY is smaller ˚ . As LY is greater than 200 A ˚ , the decrease than 200 A of atoms average energy become more and more slow, and the change has become very small when LY
Fig. 7. Length of crack for different system sizes at KI = 1.6KIC. Lc indicates the final length of crack after 4000 steps relaxation.
In our previous work, we studied that twin formed at the crack tip accompanying crack cleavage in brittle fracture at high loading rates or at low temperatures [13]. In this section we have tested different sizes of system, and try to find the influence of the size of atomistic region on twin formation in brittle cleavage fracture. The system size is just the same as which in Section 3.2. Before simulation, we at first get a starting crack at the critical stress intensity factor KIC, and atoms in the initial crack are fully relaxed. Subsequently the external loading is applied increasingly by re-scaling the displacement of atoms at the boundary region according to linear elasticity theory. For each size of system, the load increases incrementally from 1.0KIC to 2.0KIC, and the load step is 0.1KIC. At every loading process, the system is relaxed for 300 steps, thus the loading rate is 0.1KIC/300 steps. The atomic structure of the crack tip at different loads for different system sizes are analyzed below: when ˚ (about 30 layers of atoms in LY is smaller than 45 A [0 1 0] direction), crack propagation along {0 1 0} cleavage plane cannot be observed. Fig. 8(a)–(c) show the atomic configuration of the crack tip for the system that 20 layers of atoms are contained in [0 1 0] direction. In Fig. 8(a) we can find the breaking of bonds at the crack tip at 1.2KIC. But with the load increasing, the crack tip does not move any more. When the load increases to 1.6KIC, the crack tip is blunted as shown in Fig. 8(b), and the crack cleavage along [0 1 0] direction, which is perpendicular to crack surface. ˚ , crack propagation When LY is greater than 45 A along cleavage plane is observed. Fig. 8(d)–(i) show the atomic configurations of crack tips while 40 and 98 atomic layers are contained along [0 1 0] direction. As shown in Fig. 8(d)–(f), the breaking of bond at the crack tip is observed at 1.2KIC, and the crack tip moves continuously along {0 1 0} cleavage plane with the load increasing. Meanwhile, stacking fault is found to form at the crack tip while crack tip moves, but no twin strip observed. When the number of atomic layer is greater ˚ ), the crack begins to cleavage at than 98 (about 140 A 1.1KIC. As shown in Fig. 8(h) and (i), twin formation is observed at the crack tip accompanying crack cleavage
438
Y.-F. Guo, Y.-C. Gao / Computational Materials Science 36 (2006) 432–439
˚ , KI = 1.2KIC; (b) LY = 29 A ˚ , KI = 1.6KIC; (c) Fig. 8. Atomic configuration of crack at different loads for different system sizes: (a) LY = 29 A ˚ , KI = 1.9KIC; (d) LY = 57 A ˚ , KI = 1.2KIC; (e) LY = 57 A ˚ , KI = 1.6KIC; (f) LY = 57 A ˚ , KI = 1.9KIC; (g) LY = 140 A ˚ , KI = 1.1KIC; (h) LY = 29 A ˚ , KI = 1.6KIC; (i) LY = 140 A ˚ , KI = 1.9KIC. LY indicates the length of atomistic region in Y-direction. NZ = 140 A
propagation, which is agreed with the results that we obtained in our previous studies [13]. Moreover, we also compared the atomic configuration of crack tip at 1.9 KIC while LY is, respectively, at ˚ , 280 A ˚ and 360 A ˚ . We find, for the 140 A ˚ system, 140 A both the width and length of twin are a little smaller ˚ and 360 A ˚ systems. While LY is than those for 280 A ˚ , the width and length of twin at the greater than 280 A crack tip change no more. Additionally, it is also found that the length of the crack do not change anymore ˚ . It indicates that both when LY is greater than 280 A the atomic configuration and the length of crack do not change distinctly with the system size increasing if ˚ or so. Thus we can conclude LY is greater than 280 A that in order to reflect the real process of brittle fracture, the diameter of atomistic region cannot be smaller than ˚. 300 A
4. Conclusions In this paper, we first investigate the displacement distributions at a crack tip both with continuum and discrete method. The results indicate that the proper radius ˚ , that of the discrete region at a crack tip is about 120 A is also the truncated distance of continuum elastic field around a crack tip. Subsequently, the influence of the
size of atomistic region on brittle fracture for bcc-iron crystal is studied by using combined continuum-atomistic method. The atomistic model is used for the region near the crack tip while the continuum model is used for the outside domain, and the solution of the displacement field for opening crack according to anisotropic continuum fracture mechanics provides the boundary condition of the atomistic region. The results indicate that the necessary minimum size of atomistic model ˚ . When the size of for brittle fracture is about 300 A atomistic region is smaller than this value, the elastic field of continuum mechanics cannot reflect the atomic processes at the crack tip properly. Thus we assume the system size that can bridge continuum with discrete region is about 30 nm. The elastic continuum mechanics cannot describe materials brittle fracture behavior effectively when the system size is smaller than this value, and the discrete properties should be considered while the material behavior at nanometer scale is concerned.
Acknowledgments The research was supported by Chinese Nature Science Foundation (No. 10228204), NJTU Science Foundation (No. 2004RC006) and NJTU Paper
Y.-F. Guo, Y.-C. Gao / Computational Materials Science 36 (2006) 432–439
Foundation. The authors thank Prof. Yuesheng Wang and Dr. Suowen Gao for their helpful discussions.
References [1] P. Gumbsch, An atomistic study of brittle fracture: toward explicit failure criteria from atomistic modeling, J. Mater. Res. 10 (1995) 2897. [2] H. Rafii-Tabar, L. Hua, M. Cross, A multi-scale atomisticcontinuum modelling of crack propagation in a two-dimensional macroscopic plate, J. Phys. 10 (1998) 2375. [3] J.Q. Broughton, F.F. Abraham, N. Bernstein, E. Kaxiars, Concurrent coupling of length scales: methodology and application, Phys. Rev. B 60 (1999) 2391. [4] S.V. Bhenoy, R. Miller, E.B. Tadmor, D. Rodney, R. Phillips, M. Ortiz, An adaptive finite element approach to atomic-scale mechanics—the quasicontinuum method, J. Mech. Phys. Solids 47 (1999) 611. [5] A.A. Griffith, The phenomena of rupture and flow in solids, Philos. Trans. Roy. Soc. A 221 (1920) 163. [6] A. George, G. Michot, Dislocation loops at crack tips: nucleation and growth—an experimental study in silicon, Mater. Sci. Eng. A 164 (1993) 118. [7] J. Riedle, P. Gumbsch, H.F. Fischmeister, Cleavage anisotropy in tungsten single crystals, Phys. Rev. Lett. 76 (1996) 3594. [8] B. DeCelis, A.S. Argon, S. Yip, Molecular dynamics simulation of crack tip processes in alpha-iron and copper, J. Appl. Phys. 54 (1983) 4864. [9] K.S. Cheng, S. Yip, Brittle–ductile transition in intrinsic fracture behavior of crystals, Phys. Rev. Lett. 65 (1990) 2804. [10] F.F. Abraham, Dynamics of brittle fracture with variable elasticity, Phys. Rev. Lett. 77 (1996) 869. [11] F. Cleri, S. Yip, D. Wolf, S.R. Phillpot, Atomistic-scale mechanism of crack-tip plasticity: dislocation nucleation and crack tip shielding, Phys. Rev. Lett. 79 (1997) 1309.
439
[12] S. Kohlhoff, P. Gumbsch, H.F. Fischmeister, Crack propagation in bcc crystals studied with a combined finite-element and atomistic model, Philos. Mag. A 64 (1991) 851. [13] Y.F. Guo, C.Y. Wang, D.L. Zhao, Atomistic simulation of crack cleavage and blunting in bcc-Fe, Mater. Sci. Eng. A 349 (2003) 29. [14] P. Gumbsch, Modelling brittle and semi-brittle fracture processes, Mater. Sci. Eng. A 319–321 (2001) 1. [15] P. Gumbsch, R.M. Cannon, Atomistic aspects of brittle fracture, MRS Bull. 25 (2000) 15. [16] F.F. Abraham, D. Schneider, B. Land, D. Lifka, J. Skovira, J. Gerner, M. Rosenkrantz, Instability dynamics in three-dimensional fracture: an atomistic simulation, J. Mech. Phys. Solids 45 (1997) 1461. [17] V. Bulatov, F.F. Abraham, L. Kubin, B. Devincre, S. Yip, Connecting atomistic and mesoscale simulations of crystal plasticity, Nature 391 (1998) 669. [18] R. Perez, P. Gumbsch, An ab initio study of the cleavage anisotropy in silicon, Acta Mater. 48 (2000) 4517. [19] J.C.H. Spence, Y.M. Huang, O. Sankey, Lattice trapping and surface reconstruction for silicon cleavage on (1 1 1). Ab initio quantum molecular dynamics calculation, Acta Metall. Mater. 41 (1993) 2815. [20] M.J. Buehler, F.F. Abraham, H. Gao, Hyperelasticity governs dynamic fracture at a critical length scale, Nature 426 (2003) 141. [21] P. Espan˜ol, I. Zu´n˜iga, M.A. Rubio, Effect of boundary conditions in mode I fracture in brittle materials, Physica D 96 (1996) 375. [22] B.L. Holian, R. Revelo, Fracture simulations using large-scale molecular dynamics, Phys. Rev. B 51 (1995) 11275. [23] G.C. Sih, H. Liebowitz, Mathematical theories of brittle facture, Fracture: An Advanced Treatise, vol. 2, Academic, New York, 1968 (Chapter 2). [24] M.W. Finnis, J.E. Sinclair, A simple N-body potential for transition metals, Philos. Mag. A 50 (1984) 45; Erratum 53 (1986) 161.