Insights on Dynamic Fracture in Brittle Materials: Molecular Dynamics, Notas de estudo de Engenharia de Produção

Baixe Insights on Dynamic Fracture in Brittle Materials: Molecular Dynamics e outras Notas de estudo em PDF para Engenharia de Produção, somente na Docsity! 1 CHAPTER 1 MODELING DYNAMIC FRACTURE USING LARGE-SCALE ATOMISTIC SIMULATIONS Markus J. Buehler Massachusetts Institute of Technology, Department of Civil and Environmental Engineering 77 Massachusetts Avenue Room 1-272, Cambridge, MA., 02139, USA E-mail: mbuehler@MIT.EDU Huajian Gao Max Planck Institute for Metals Research Heisenbergstrasse 3, D-70569 Stuttgart, Germany We review a series of large-scale molecular dynamics studies of dynamic fracture in brittle materials, aiming to clarify questions such as the limiting speed of cracks, crack tip instabilities and crack dynamics at interfaces. This chapter includes a brief introduction of atomistic modeling techniques and a short review of important continuum mechanics concepts of fracture. We find that hyperelasticity, the elasticity of large strains, can play a governing role in dynamic fracture. In particular, hyperelastic deformation near a crack tip provides explanations for a number of phenomena including the “mirror-mist- hackle” instability widely observed in experiments as well as supersonic crack propagation in elastically stiffening materials. We also find that crack propagation along interfaces between dissimilar materials can be dramatically different from that in homogeneous materials, exhibiting various discontinuous transition mechanisms (mother-daughter and mother-daughter-granddaughter) to different admissible velocity regimes. M.J. Buehler and H. Gao 2 1. Introduction Why and how cracks spread in brittle materials is of essential interest to numerous scientific disciplines and technological applications [1-3]. Large-scale molecular dynamics (MD) simulation [4-13] is becoming an increasingly useful tool to investigate some of the most fundamental aspects of dynamic fracture [14-20]. Studying rapidly propagating cracks using atomistic methods is particularly attractive, because cracks propagate at speeds of kilometers per second, corresponding to time-and length scales of nanometers per picoseconds readily accessible within classical MD methods. This similarity in time and length scales partly explains the success of MD in describing the physics and mechanics of dynamic fracture. 1.1 Brief review: MD modeling of fracture Atomistic simulations of fracture were carried out as early as 1976, in first studies by Ashurst and Hoover [21]. Some important features of dynamic fracture were described in that paper, although the simulation sizes were extremely small, comprising of only 64x16 atoms with crack lengths around ten atoms. A later classical paper by Abraham and coworkers published in 1994 stimulated much further research in this field [22]. Abraham and coworkers reported molecular-dynamics simulations of fracture in systems up to 500,000 atoms, which was a significant number at the time. In these atomistic calculations, a Lennard- Jones (LJ) potential [23] was used. The results in [22, 24] were quite striking because the molecular-dynamics simulations were shown to reproduce phenomena that were discovered in experiments only a few years earlier [25]. An important classical phenomenon in dynamic fracture was the so-called “mirror-mist-hackle” transition. It was known since 1930s that the crack face morphology changes as the crack speed increases. This phenomenon is also referred to as dynamic instability of cracks. Up to a speed of about one third of the Rayleigh-wave speed, the crack surface is atomically flat (mirror regime). For higher crack speeds the crack starts to roughen (mist regime) and eventually becomes very rough (hackle regime), accompanied by extensive crack branching and Modeling dynamic fracture using large-scale atomistic simulations 5 prescribed path in Section 3, they are completely unconstrained in Section 4. Section 5 contains studies on both constrained and unconstrained crack propagation. We conclude with a discussion and outlook to future research in this area. 2. Large-scale atomistic modeling of dynamic fracture: A fundamental viewpoint 2.1 Molecular dynamics simulations Our simulation tool is classical molecular-dynamics (MD) [23, 49]. For a more thorough review of MD and implementation on supercomputers, we refer the reader to other articles and books [19, 23, 50-52]. Here we only present a very brief review. MD predicts the motion of a large number of atoms governed by their mutual interatomic interaction, and it requires numerical integration of Newton’s equations of motion usually via a Velocity Verlet algorithm [23] with time step tΔ on the order of a few femtoseconds. 2.2 Model potentials for brittle materials: A simplistic but powerful approach The most critical input parameter in MD is the choice of interatomic potentials [23, 49]. In the studies reported in this article, the objective is to develop an interatomic potential that yields generic elastic behaviors at small and large strains, which can then be linked empirically to the behavior of real materials and allows independent variation of parameters governing small-strain and large-strain properties. Interatomic potentials for a variety of different brittle materials exist, many of which are derived form first principles (see, for example [53- 55]). However, it is difficult to identify generic relationships between potential parameters and macroscopic observables such as the crack limiting or instability speeds when using such complicated potentials. We deliberately avoid these complexities by adopting a simple pair potential based on a harmonic interatomic potential with spring constant M.J. Buehler and H. Gao 6 0k . In this case, the interatomic potential between pairs of atoms is expressed as ( )2000 2 1)( rrkar −+=φ , (1) where 0r denotes the equilibrium distance between atoms, for a 2D triangular lattice, as schematically shown in the inlay of Figure 1. This harmonic potential is a first-order approximation of the Lennard-Jones 12:6 potential [23], one of the simplest and most widely used pair potentials defined as ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ ⎥⎦ ⎤ ⎢⎣ ⎡−⎥⎦ ⎤ ⎢⎣ ⎡= 612 4)( rr r σσεφ . (2) We express all quantities in reduced units, so that lengths are scaled by the LJ-parameter σ which is assumed to be unity, and energies are scaled by the parameter ε , the depth of the minimum of the LJ potential. We note that corresponding to choosing 1=ε in eq. (2), we choose 10 −=a in the harmonic approximation shown in eq. (1). Further, we note that 12246.1260 ≈=r (see Figure 1). Note that the parameter σ is coupled to the lattice constant a as 6 2/2/a=σ . The reduced temperature is ε/k TB with Bk being the Boltzmann constant. The mass of each atoms in the models is assumed to be unity, relative to the reference mass *m . The reference time unit is then given by εσ /2** mt = . For example, when choosing electron volt as reference energy ( 1=ε eV), Angstrom as reference length ( 1=σ Å), and the atomic mass unit as reference mass ( 1* =m amu), the reference time unit corresponds to 14E018051* -.t ≈ seconds. Although the choice of a simple harmonic interatomic potential can not lead to quantitative calculations of fracture properties in a particular material, it allows us to draw certain generic conclusions about fundamental, material-independent mechanisms that can help elucidating the physical foundations of brittle fracture. The harmonic potential leads to linear elastic material properties, and thus serves as the starting point when comparing simulation results to predictions by the classical linear theories of fracture. Modeling dynamic fracture using large-scale atomistic simulations 7 Various flavors and modifications of the simplistic, harmonic model potentials as described in eq. (1) are used for the studies reviewed in this article. These modifications are discussed in detail in each section. The concept of using simplistic model potentials to understand the generic features of fracture was pioneered by Abraham and coworkers [10, 12, 22, 24, 31, 40, 56-58]. 2.3 Model geometry and simulation procedure We typically consider a crack in a two-dimensional geometry with slab size yx ll × and initial crack length a , as schematically shown in Fig. 1. The slab size is chosen large enough such that waves reflected from the boundary do not interfere with the propagating crack at early stages of the simulation. The slab size in the crack direction is chosen between 2 and 4 times larger than the size orthogonal to the initial crack. The slab is initialized at zero temperature prior to simulation. Most of our studies are carried out in a two-dimensional, hexagonal lattice. The three-dimensional studies are performed in a FCC lattice, with periodic boundary conditions in the out-of-plane z direction. The slab is slowly loaded with a constant strain rate of xε& (corresponding to tensile, mode I loading), and/or xyε& (corresponding to shear, mode II loading). We establish a linear velocity gradient prior to simulation to avoid shock wave generation from the boundaries. While the loading is applied, the stress ijσ (corresponding to the specified loading case) steadily increases leading to a slowly increasing crack velocity upon fracture initiation. In the case of anti-plane shear loading (mode III), the load is applied in a similar way. It can be shown that the stress intensity factor remains constant in a strip geometry inside a region of crack lengths a [59] )( 4141 xyx llal −<< . (3) Accurate determination of crack tip velocity is crucial as we need to be able to measure small changes in the propagation speed. The crack tip position a is determined by finding the surface atom with maximum y position in the interior of a search region inside the slab using a potential energy criterion. This quantity is averaged over small time M.J. Buehler and H. Gao 10 speed observed in some experiments and many computer simulations. However, it is not generally accepted that hyperelasticity should play a significant role in dynamic fracture. One reason for this belief stems from the fact that the zone of large deformation in a loaded body with a crack is highly confined to the crack tip, so that the region where linear elastic theory does not hold is extremely small compared to the extensions of the specimen [14, 15]. This is demonstrated in Figure 2(a). We have performed large-scale molecular-dynamics studies in conjunction with continuum mechanics concepts to demonstrate that hyperelasticity can be crucial for understanding dynamic fracture. Our study shows that local hyperelasticity around the crack tip can significantly influence the limiting speed of cracks by enhancing or reducing local energy flow. This is true even if the zone of hyperelasticity is small compared to the specimen dimensions. The hyperelastic theory completely changes the concept of the maximum Fig. 3. Concept of increasingly strong hyperelastic effect (subplot (a)). The increasingly strong hyperelastic effect is modeled by using biharmonic potentials (subplot (b)). The bilinear or biharmonic model allows to tune the size of the hyperelastic region near a moving crack, as indicated in subplots (c) and (d) [30]. The local increase of elastic modulus and thus wave speeds can be tuned by changing the slope of the large-strain stress-strain curve (“local modulus” [30, 62, 63]). Modeling dynamic fracture using large-scale atomistic simulations 11 crack velocity in the classical theories. For example, the classical theories [14, 15] clearly predict that mode I cracks are limited by the Rayleigh wave speed and mode II cracks are limited by the longitudinal wave speed. In mode III, theory predicts that the limiting speed of cracks is the shear wave speed [14, 15]. In contrast, super-Rayleigh mode I and supersonic mode II cracks are allowed by hyperelasticity and have been seen in computer simulations [13, 31]. We have also observed mode III cracks faster than the shear wave speed [66]. 3.2 Strategy of investigation Our strategy of study is to first establish harmonic reference systems with behaviors perfectly matching those predicted by the existing linear elastic theories of fracture. Using a biharmonic model potentials [30], we then introduce increasingly stronger nonlinearities and show that under certain conditions, the linear elastic theory breaks down as the material behaviour near crack tips deviates increasingly from the linear elastic approximation. This is further visualized in Figure 3. In comparison with experiments, a major advantage of computer simulations is that they allow material behaviors to be fine tuned by introducing interatomic potentials that focus on specific aspects, one at a time, such as the large-strain elastic modulus, smoothing at bond breaking and cohesive stress. From this point of view, MD simulations can be regarded as computer experiments that are capable of testing theoretical concepts and identifying controlling factors in complex systems. 3.3 Elastic properties: The link between the atomistic scale and continuum theory 3.3.1 The virial stress and strain Stresses are calculated according to the virial theorem [67, 68]. The atomistic stress is given by (note that Ω is the atomic volume) M.J. Buehler and H. Gao 12 ⎟⎟ ⎟ ⎠ ⎞ ⎜⎜ ⎜ ⎝ ⎛ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎣ ⎡ ∂ ∂ − Ω = = ∑ αβ βα φσ rr jiij rrrr, 1 2 11 . (4) We also use the concept of an atomic strain to link the atomistic simulation results with continuum theory. For a 2D triangular lattice, we define the left Cauchy-Green strain tensor of an atom l [69] ( )( )[ ] ⎭ ⎬ ⎫ ⎩ ⎨ ⎧ −−= ∑ = N k k j l j k i l i l ij xxxxr b 1 2 03 1 , (5) with 6=N nearest neighbors, and the variable lix denoting the position vector of atom l . Unlike the virial stress, the virial strain is valid instantaneously in space and time. The maximum principal strain 1b is obtained by diagonalization of the strain tensor ijb . The maximum principal engineering strain is then given by 111 −= bε . The definition of atomic stress and strain allows an immediate comparison of MD results with the prediction by continuum theory, such as changes in deformation field as a function of increasing crack speed [14]. This has been discussed extensively in the literature [29, 37, 51, 60, 61, 70, 71]. Some of these results will be reviewed later. Fig. 4. Elastic properties of the two-dimensional model system with harmonic potential, including (a) stress versus strain and Poisson’s ratio, and (b) tangent Young’s modulus as a function of strain. The figure shows the results for two different crystal orientations (see left part of the figure) [70]. Modeling dynamic fracture using large-scale atomistic simulations 15 from continuation conditions. The elastic properties associated with this potential are shown in Figure 5 for uniaxial stress in the x and y- directions. This potential allows to smoothly interpolate between harmonic potentials and strongly nonlinear potentials by changing the parameter onr and/or the ratio of the spring constants 01 / kk . 3.3.4 Fracture surface energy The fracture surface energy γ is an important quantity for nucleation and propagation of cracks. It is defined as the energy required to generate a unit distance of a pair of new surfaces (cracks can be regarded as sinks for energy, where elastic energy is converted into surface fracture energy). The Griffith criterion predicts that the crack tip begins to propagate when the crack tip energy release rate G reaches the fracture surface energy γ2=G [72]. Knowledge of the atomic lattice and the interatomic potential can also be used to define the fracture surface energy. This quantity depends on the crystallographic directions and is calculated to be d φγ Δ−= , (13) where φΔ the energy necessary to break atomic bonds as the crack advances a distance d . For the harmonic bond snapping potential as described above, the fracture surface energy is given by 0332.0≈γ for the direction of high surface energy (Figure 1(b)), and 0288.0≈γ for the other direction (Figure 1(c); about 15 % smaller than in the other direction). 3.4 Harmonic reference systems After briefly defining the atomistic models in the previous sections, we now discuss simulation results using the harmonic potentials as the reference system. M.J. Buehler and H. Gao 16 Our results show that the harmonic systems behave as predicted by classical, linear elastic continuum theories of fracture [14]. For a mode I tensile crack, linear theory predicts that the energy release rate vanishes for all velocities in excess of the Rayleigh wave speed, implying that a mode I crack cannot move faster than the Rayleigh wave speed. This prediction is indeed confirmed in systems with the harmonic potential. The crack velocity approaches the Rayleigh wave speed independent of the slab size, provided that the applied strain is larger than 1.08 percent and the slab width is sufficiently large ( 000,1>xl ). Results for a mode I crack and two different spring constants are shown in Figure 6. In other studies, we find that the crack limiting speed of mode II cracks is the longitudinal wave speed. Mode III cracks are limited by the shear wave speed (details and results of those simulations are not Fig. 6. Crack tip history (upper part) as well as the crack speed history (lower part) for a soft as well as a stiff harmonic material [29, 35]. Modeling dynamic fracture using large-scale atomistic simulations 17 described here). These results agree well with linear elastic fracture mechanics theories [14]. 3.5 Biharmonic simulations: Elastic stiffening Now we describe a series of computer experiments carried out with the biharmonic interatomic potential introduced in eq. (4) and (5). The results indicate that a local hyperelastic zone around the crack tip (as schematically shown in Figure 2(a)) can have significant effect on the velocity of the crack. We consider hyperelastic effect of different strengths by using a biharmonic potential with different onset strains governed by the parameter onr . The parameter onr governs the onset strain of the hyperelastic effect 00 /)( rrronon −=ε . The simulations reveal crack propagation at super-Rayleigh velocities in steady-state with a local stiffening zone around the crack tip. Fig. 7. Change of the crack speed as a function of εon. The smaller the εon, the larger the hyperelastic region, and the larger the crack speed [30]. These results indicate that the crack speed of mode I cracks is not limited by the Rayleigh-wave speed. M.J. Buehler and H. Gao 20 The problem of a super-Rayleigh mode I crack in an elastically stiffening material is somewhat analogous to Broberg’s [73] problem of a mode I crack propagating in a stiff elastic strip embedded in a soft matrix. The geometry of this problem is shown in Fig. 11. Broberg [73] has shown that, when such a crack propagates supersonically with respect to the wave speeds of the surrounding matrix, the energy release rate can be expressed in the form ),,( 21 2 ccvf E hG σ= (14) where σ is the applied stress, h is the half width of the stiff layer and f is a non-dimensional function of crack velocity and wave speeds in the strip and the surrounding matrix ( 21,cc ). The dynamic Griffith energy balance requires γ2=G , indicating that crack propagation velocity is a function of the ratio χh where 2/~ σγχ E can be defined as a characteristic length scale for local energy flux. By dimensional analysis, the energy release rate of our hyperelastic stiffening material is expected Fig. 10. Sequence of snapshots showing supersonic propagation of a crack under shear loading. A small localized hyperelastic region (not shown here; see reference [30]) at the crack tip leads to crack speeds faster than both shear and longitudinal wave speeds in the material [30]. Modeling dynamic fracture using large-scale atomistic simulations 21 to have similar features except that Broberg’s strip width h should be replaced by a characteristic size of the hyperelastic region Hr . Therefore, we introduce the concept of a characteristic length 2σ γβχ E= (15) for local energy flux near a crack tip. The coefficient β may depend on the ratio between hyperelastic and linear elastic properties. We have simulated the Broberg problem and found that the mode I crack speed reaches the local Rayleigh wave speed as soon as χh reaches unity. Numerous simulations verify that the scaling law in eq. (7) holds when γ , E and σ is changed independently. The results are shown in Figure 12. From the simulations, we estimate numerically 87≈β and the characteristic energy length scale 750≈χ . The existence of a characteristic length for local energy flux near the crack tip has not been discussed in the literature and plays the central role in understanding the effect of hyperelasticity. Under a particular experimental or simulation condition, the relative importance of hyperelasticity is determined by the ratio χ/Hr . For small χ/Hr , the crack dynamics is dominated by the global linear elastic properties since much of the energy transport necessary to sustain crack motion occurs in the linear elastic region. However, when χ/Hr approaches unity, as is Fig. 11. Geometry of the Broberg problem [30, 73], consisting of a crack embedded in a stiff layer (high wave velocities, Young’s modulus E1) surrounded by soft medium (low wave velocities, Young’s modulus E0). M.J. Buehler and H. Gao 22 the case in some of our molecular dynamics simulations, the dynamics of the crack is dominated by local elastic properties because the energy transport required for crack motion occurs within the hyperelastic region. The concept of energy characteristic length χ immediately provides an explanation how the classical barrier for transport of energy over large distances can be undone by rapid transport near the tip. 3.7 Discussion and Conclusions We have shown that local hyperelasticity has a significant effect on the dynamics of brittle crack speeds and have discovered a characteristic length associated with energy transport near a crack tip. The assumption of linear elasticity fails if there is a hyperelastic zone in the vicinity of the crack tip comparable to the energy characteristic length. Therefore, we conclude that hyperelasticity is crucial for understanding and Fig. 12. Calculation results of the geometry defined in the Broberg problem [30]. The plot shows results of different calculations where the applied stress, elastic properties and fracture surface energy are independently varied. In accordance with the concept of the characteristic energy length scale, all points fall onto the same curve and the velocity depends only on the ratio χ/h . Modeling dynamic fracture using large-scale atomistic simulations 25 Most existing theories of fracture assume a linear elastic stress-strain law. However, from a hyperelastic point of view, the relation of stress and strain in real solids is strongly nonlinear near a moving crack tip. Using massively parallel large-scale atomistic simulations, we show that hyperelasticity also plays an important role in dynamical crack tip instabilities. We find that the dynamical instability of cracks can be regarded as a competition between different instability mechanisms Fig. 14. Crack propagation in a Lennard-Jones (see eq. (2)) system as reported earlier [22, 24]. Subplot (a) shows the σxx-field and indicates the dynamical mirror-mist-hackle transition as the crack speed increases. The crack velocity history (normalized by the Rayleigh-wave speed) is shown in subplot (b). M.J. Buehler and H. Gao 26 controlled by local energy flow and local stress field near the crack tip. We hope the simulation results can help explain controversial experimental and computational results. 4.1 Introduction Here we focus on the instability dynamics of rapidly moving cracks. In several experimental [25, 64, 74, 75] and computational [20, 22, 26, 76, 77] studies, it was established that the crack face morphology changes as the crack speed increases, a phenomenon which has been referred to as the dynamical instability of cracks. Up to a critical speed, newly created crack surfaces are mirror flat (atomistically flat in MD simulations), whereas at higher speeds, the crack surfaces start to roughen (mist regime) and eventually becomes very rough (hackle regime). This is found to be a universal behavior that appears in various brittle materials, including ceramics, glasses and polymers. This dynamical crack instability has also been observed in computer simulation [20, 22, 26, 76, 77]. The result of a large-scale MD simulation illustrating the mirror-mist-hackle transition is shown in Fig. 14. Despite extensive studies in the past, previous results have led to numerous discrepancies that remain largely unresolved up to date. For example, linear elasticity analyses first carried out by Yoffe [78] predict that the instability speed of cracks is about 73% of the Rayleigh-wave speed [14, 15], as the circumferential hoop stress exhibits a maximum at an inclined cleavage plane for crack speeds beyond this critical crack speed. This is in sharp contrast to observations in several experiments and computer simulations. Experiments have shown that the critical instability speed can be much lower in many materials. In 1992, Fineberg et al. [25, 64] observed an instability speed at about one third of the Rayleigh wave speed, which significantly deviates from Yoffe’s theory [78]. Similar observations were made in a number of other experimental studies in different classes of materials (e.g. crystalline silicon, polymers such as PMMA). The mirror-mist-hackle transition at about one third of the Rayleigh-wave speed was also observed in the large-scale MD simulations carried out by Abraham et al. in 1994 [22] Modeling dynamic fracture using large-scale atomistic simulations 27 (see Fig. 14). The instability issue is further complicated by recent experiments which show stable crack propagation at speeds even beyond the shear wave speed in rubber-like materials [79]. The existing theoretical concepts seem insufficient to explain the various, sometimes conflicting, results and observations made in experimental and numerical studies. The dynamical crack instability has been a subject of numerous theoretical investigations in the past decades, and several theoretical explanations have been proposed. First, there is Yoffe’s model [78] which shows the occurrence of two symmetric peaks of normal stress on inclined cleavage planes at around 73% of the Rayleigh-wave speed. Later, Gao [80] proved that Yoffe’s model is consistent with a criterion of crack kinking into the direction of maximum energy release rate. Eshelby [81] and Freund [82] made an argument that the dynamic energy release rate of a rapidly moving crack allows the possibility for the crack to split into multiple branches at a critical speed of about 50% of the Raleigh speed. Marder and Gross [26] presented an analysis which Fig. 15. The concept of hyperelastic softening close to bond breaking, in comparison to the linear elastic, bond-snapping approximation. M.J. Buehler and H. Gao 30 that the crack tip instability speeds for a wide range of materials behaviors fall between the predictions of these two models. In the spirit of “model materials” as introduced in Section 3, we develop a new, simple material model which allows a systematic transition from linear elastic to strongly nonlinear material behaviors, with the objective to bridge different existing theories and determine the conditions of their validity. By systematically changing the large-strain elastic properties while keeping the small-strain elastic properties constant, the model allows us to tune the size of hyperelastic zone and to probe the conditions under which the elasticity of large strains governs the instability dynamics of cracks. In the case of linear elastic model with bond snapping, we find that the instability speed agrees well with the predicted value from Yoffe's model. We then gradually tune up the hyperelastic effects and find that the instability speed increasingly agrees with Gao’s model. In this way, we achieve, for the first time, a unified treatment of the instability problem leading to a generalized model that bridges Yoffe’s linear elastic branching model to Gao’s hyperelastic model. 4.2 Atomistic modelling Although simple pair potentials do not allow drawing conclusions for unique phenomena pertaining to specific materials, they enable us to understand universal, generic relationships between potential shape and fracture dynamics in brittle materials; in the present study we use a simple pair potential that allows the hyperelastic zone size and cohesive stress to be tuned. The potential is composed of a harmonic function in combination with a smooth cut-off of the force based on the Fermi-Dirac (F-D) distribution function to describe smooth bond breaking. We do not include any dissipative terms. The force versus atomic separation is expressed as 1 00 1exp)()( − ⎥ ⎦ ⎤ ⎢ ⎣ ⎡ +⎟⎟ ⎠ ⎞ ⎜⎜ ⎝ ⎛ Ξ− Ξ −= breakr rrrkr dr dφ (17) Modeling dynamic fracture using large-scale atomistic simulations 31 Assuming that the spring constant 0k is fixed, the potential has two additional parameters, breakr and Ξ . The parameter breakr (corresponding to the Fermi energy in the F-D-function) denotes the critical separation for breaking of the atomic bonds and allows tuning the breaking strain as well as the cohesive stress at bond breaking. In particular, we note that breakcoh r~σ . (18) The parameter Ξ (corresponding to the temperature in the F-D- function) describes the amount of smoothing at the breaking point. In addition to defining the small-strain elastic properties (by changing the parameter 0k ), the present model allows us to control the two most critical physical parameters describing hyperelasticity, (i) cohesive stress (by changing the parameter breakr ), and (ii) the strength of softening close to the crack tip (by changing the parameter Ξ ). Figure 16 depicts force versus atomic separation of the interatomic potential used in our study. The upper part shows the force versus separation curve with respect to changes of breakr , and the lower part shows the variation in shape when Ξ is varied. For small values of Ξ (around 50), the softening effect is quite large. For large values of Ξ (beyond 1,000), the amount of softening close to bond breaking becomes very small, and the solid behaves like one with snapping bonds. The parameter breakr allows the cohesive stress cohσ to be varied Fig. 17. Crack propagation in a homogeneous harmonic solid. When the crack reaches a velocity of about 73 percent of Rayleigh wave speed, the crack becomes unstable in the forward direction and starts to branch (the dotted line indicates the 60° plane of maximum hoop stress) [60]. M.J. Buehler and H. Gao 32 independently. This model potential also describes the limiting cases of material behavior corresponding to Yoffe’s model (linear elasticity with snapping bonds) and Gao’s model (strongly nonlinear behavior near the crack tip). Yoffe’s model predicts that the instability speed only depends on the small-strain elasticity. Therefore, the instability speed should remain constant at 73% of the Rayleigh-wave speed, irregardless of the choices of the parameters breakr and Ξ . On the other hand, Gao’s model predicts that the instability speed is only dependent on the cohesive stress cohσ (and thus breakr ): ρρ σ breakcohGao inst r v ∝= (19) (note that ρ denotes the density, as defined above). According to Gao’s model, variations in the softening parameter Ξ should not influence the crack instability speed. In most of earlier computational studies, the analysis was performed only for a single interatomic potential, such as done for a LJ potential by Abraham et al. [22]. It remains unclear how the nature and shape of the atomic interaction affects the instability dynamics. Here we propose systematic numerical studies based on continuously varying potential Fig. 18 Comparison between hoop stresses calculated from molecular-dynamics simulation with harmonic potential and those predicted by linear elastic theory [14] for different reduced crack speeds v/cR [29, 60]. The plot clearly reveals development of a maximum hoop stress at an inclined angle at crack speeds beyond 73% of the Rayleigh-wave speed. Modeling dynamic fracture using large-scale atomistic simulations 35 The predictions by Yoffe’s model are included in Fig. 19 as the red line, and the predictions by Gao’s model are plotted as the blue points. As seen in Fig. 19, we observe that for any choice of breakr and Ξ , the instability speed lies in between the prediction by Gao’s model and that by Yoffe’s model. Whether it is closer to Gao’s model or to Yoffe’s model depends on the choice of breakr and Ξ . For small values of breakr and Ξ , we find that the instability speed depends on the cohesive stress, which is a feature predicted by Gao’s model. Further, the instability speed seems to be limited by the Yoffe speed (as can be confirmed in Fig. 19 for 150=Ξ and 300=Ξ for large values of breakr ). Whereas the observed limiting speeds increase with breakr for 22.1<breakr , they saturate at the Yoffe speed of 73% of Rayleigh-wave speed for larger values of 22.1≥breakr (see Fig. 19, bottom curve for 300=Ξ ). In this case, the instability speed is independent of breakr and independent of Ξ . This behavior is reminiscent of Yoffe’s deformation controlled instability mechanism and suggests a change in governing mechanism for the instability speed: The instability speed becomes governed by a Yoffe-like deformation field mechanism for 22.1≥breakr , whereas the instability seems to be influenced by the cohesive-stress for 22.1<breakr . We find a similar transition for different choices of Ξ ranging from 50=Ξ to 500,1=Ξ , with different values of breakr at which the transition is observed. These results indicate that the instability speed depends on the strength of softening (parameter Ξ ) and on the cohesive stress close to bond breaking (parameter breakr ). 4.4 The modified instability model We observe that the first derivative of the instability speed with respect to the cohesive stress in our MD simulations agrees reasonably well with Gao’s model. However, the observed instability speed differs from Gao’s model by a constant value which seems to depend on the softening parameter Ξ . We measure the deviation from Gao’s model by a shift parameter Gaoinst MD instshift vvv −= which is a function of Ξ . The curve is shown in Figure 20. M.J. Buehler and H. Gao 36 The physical interpretation of shiftv is that it accounts for the relative importance of hyperelastic softening close to the crack tip: Gao’s model corresponds to the limiting case when the softening region is large and completely dominates the energy flow, and it therefore constitutes the lower limit for the instability speed. Indeed, we find that for very strong softening, that is, when 0→Ξ , shiftv vanishes (Figure 20). In contrast, shiftv assumes larger values in the case of vanishing softening as ∞→Ξ . The physical significance of this parameter can be understood from the perspective of the characteristic energy length scale of dynamic fracture 2~ σ γχ E (20) proposed earlier (see Section 3) [30]. The characteristic energy length scale χ describes the region from which energy flows to the crack tip to drive its motion. If the size of the softening region is comparable to χ , i.e., Fig. 20 Change of shift parameter shiftv as a function of the smoothing parameter Ξ . Small values of Ξ lead to cases where hyperelasticity dominates the instability onset ( 1/ ≈χHr ), and larger values of Ξ correspond to the case when the importance of hyperelasticity vanishes and the instability onset is governed by changes in the deformation field ( 0/ →χHr ) Modeling dynamic fracture using large-scale atomistic simulations 37 1>> χ Hr , (21) hyperelasticity dominates energy flow, thus 0→shiftv , and the predictions of Gao’s model should be valid. In contrast, if the size of the softening region is smaller than χ , i.e. 0≈ χ Hr , (22) hyperelasticity plays a reduced role in the instability dynamics and the purely hyperelastic model becomes increasingly approximate, so that shiftv takes larger values and eventually Yoffe’s model of deformation field induced crack kinking begins to play a governing role. The shift parameter )(Ξshiftv should therefore depend on the relative size of the hyperelastic region compared to the characteristic energy flow ⎟⎟ ⎠ ⎞ ⎜⎜ ⎝ ⎛ = χ H shift rfv . (23) With the new parameter shiftv , the instability speed can be expressed as ρ σ χ cohH shiftinst rvv +⎟⎟ ⎠ ⎞ ⎜⎜ ⎝ ⎛ = (24) Based on the results of our numerical experiments, we propose that eq. (19) and the Yoffe model should be combined to predict the onset of instability at a critical speed. The critical crack tip instability speed is then given by ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ +⎟⎟ ⎠ ⎞ ⎜⎜ ⎝ ⎛ = ρ σ χ cohH shiftYOFFEinst r vvv ,minmod (25) where RYOFFE cv 73.0≈ is a constant, independent on the hyperelastic properties. We refer to this model as the modified instability model. Fig. 21 compares the predictions by eq. (25) with the MD simulation results (dashed curve). M.J. Buehler and H. Gao 40 hyperelastic effect). As can be clearly seen in Fig. 23 for 4/ 01 =kk , the instability speed can even be super-Rayleigh and approach intersonic speeds. This is inconsistent with the classical theories but can be understood from a hyperelasticity point of view. This observation suggests that the stiffening materials behavior tends to have a stabilizing effect on straight crack motion. 4.6 Discussion and conclusions Our simulation studies strongly suggest that the large-deformation elastic properties play a critical role in the instability dynamics of cracks. Our results indicate that the onset of instability can be understood as a competition between energy flow governed instability (Gao’s model Fig. 23. MD simulation results of instability speed for stiffening materials behaviour, showing stable super-Rayleigh crack motion as observed in recent experiment (see the legends in the lower left). Such observation is in contrast to any existing theories, but can be explained based on the hyperelastic viewpoint. Modeling dynamic fracture using large-scale atomistic simulations 41 [62, 63]) and stress field governed instability (Yoffe’s model [78]), as summarized in eq. (12). To the best of our knowledge, Eq. (25) describes for the first time this transition between Yoffe’s and Gao’s model of instability dynamics (see Fig. 21). We hypothesize that the transition between the two mechanisms depends on the relative importance of hyperelasticity around the crack tip, as described by the ratio of the size of the hyperelastic region and the characteristic energy length scale χ/Hr . In most experiments and computers simulations, materials show a strong softening effect. Therefore, hyperelastic softening near the crack tip may explain the reduced instability speed observed in experiments. In the case of a locally stiffening hyperelastic region, we find that the critical speeds for crack instability could be higher since the local hyperelastic zone with local higher wave speed allows for faster energy transport and higher local wave speeds; in this case, the small-strain Yoffe speed is no longer a barrier due to local stiffening. This concept provides a feasible explanation of recent experimental observations of stable mode I cracking at crack velocities beyond the shear wave speed Fig. 24. Importance of accurate large-strain properties of bond breaking during fracture in silicon [85-87]. The results suggest that the large-strain properties close to bond breaking is critical in determining how and under which conditions cracks propagate. Subplot (a) shows results without correct large-strain elasticity, whereas subplot (b) depicts the results with correct large-strain elasticity at bond breaking. The difference of bond breaking is shown in subplot (c), indicating the incorrect behavior of the Tersoff potential at bond breaking. M.J. Buehler and H. Gao 42 in homogeneous materials that show a strong stiffening hyperelastic effect (see Fig. 23). The notion that large-strain elastic properties close to bond breaking are important has recently also been demonstrated in studies of fracture in silicon [86]. Figure 24 shows a comparison of fracture dynamics by changing the description of bond breaking from an empirical Tersoff- type [88] to a first principles description using reactive potentials [89]. Only if a correct quantum mechanical derived description of bond breaking is used, the MD model is capable to reproduce important experimental features of fracture of silicon. The work reported here, together with the results on the maximum crack speed (see Section 3), strongly suggest that hyperelasticity can play governing roles in dynamic fracture. 5. Dynamic crack propagation along interfaces between dissimilar materials: Mother-, daughter- and granddaughter cracks Cracking along interfaces has received significant attention in the past decades, both due to technological relevance for example in composite materials, as well as because of scientific interest. In this section, we review recent progress in large-scale atomistic studies of crack propagation along interfaces of dissimilar materials [37, 66, 90]. We consider two linear-elastic material blocks bound together with a weak potential whose bonds snap early upon a critical atomic separation, as shown schematically in Fig. 25. This approach confines crack motion along the interface. In the two blocks, atoms interact with harmonic potentials with different spring constants adjacent to the interface. An initial crack is introduced along the interface and subjected to tensile (section 5.2) and shear (section 5.3) dominated displacement loading along the upper and lower boundaries of the sample. Under tensile loading, we observe that upon initiation at a critical load, the crack quickly approaches a velocity a few percent larger than the Rayleigh-wave speed of the soft material. After a critical time, a secondary crack is nucleated a few atomic spacings ahead of the crack. This secondary crack propagates at the Rayleigh-wave speed of the stiff material. If the elastic mismatch is sufficiently large (e.g. ten as in our Modeling dynamic fracture using large-scale atomistic simulations 45 To model a crack along an interface between two elastically dissimilar materials, we study two crystal blocks described by harmonic interatomic potentials with different spring constants 01 kk > . For simplicity, we consider two-dimensional triangular lattices with isotropic elastic properties. The harmonic potential is defined as in eq. (1). We choose 57.283/36 30 ≈=k with the corresponding Rayleigh- wave speed 4.3≈Rc , shear wave speed 68.3≈sc and longitudinal wave speed 36.6≈lc . The difference in Young’s modulus of the two materials can be related to the ratio 0 1 k k =Ξ . (26) The wave speeds Ic (with ),,( lsRI = ) of the two materials thus differ by a factor Ξ : 0,1, II cc ⋅Ξ= (27) Atomic bonds across the weak interface have spring constant 0k and break upon a critical separation 17.1=breakr . Fig. 28. Stress fields σxx , σyy and σxy for a crack at an interface with elastic mismatch Ξ =10 after the secondary crack is nucleated. In all stress fields, the two Mach cones in the soft material are seen. The red color corresponds to large stresses, and the blue color to small stresses. M.J. Buehler and H. Gao 46 The simulation slab is subjected to an applied tensile and/or shear strain rate along the upper and lower boundaries, as shown in Fig. 25. The system size is given by 298,2≈xl and 398,4≈xl . Further details on modeling fracture based on harmonic potentials, as well as a detailed analysis on the elastic and fracture properties can be found elsewhere [13, 30, 31, 38] (see also discussion in previous sections). In homogeneous materials, it has been shown that the simulated results of limiting crack speeds, the energy flow fields and the stress fields near rapidly propagating cracks are in good agreement with continuum mechanics predictions (see results discussed in Section 3). The geometry of our atomistic model, the loading as well as the crystal orientation of crack propagation is described in Fig. 25. The upper (left) block of the simulation slab is relatively stiff with higher Young’s modulus and higher wave speeds than the lower (right) block. Fig. 29. The potential energy field for a crack at a biomaterial interface with elastic mismatch Ξ =10. Two Mach cones in the soft solid can clearly be observed. The mother (A) and daughter crack (B) are marked in the figure. The mother crack is represented by a surface wave in the soft material after nucleation of the secondary crack. Modeling dynamic fracture using large-scale atomistic simulations 47 5.2 Tensile dominated cracks at bimaterial interfaces: Mother- daughter mechanisms and supersonic fracture In all cases studied in this paper, the cracked bi-crystal is loaded under tensile loading with a strain rate 1,000.0=xxε& (the orientation of the loading is shown schematically in Figure 25) [37, 90]. In some cases, the loading is stopped as indicated in each case. Although there is usually no pure mode I interfacial cracks due to complex stress intensity factors, the stress field near the crack tip is expected to be mode I dominated under pure tensile loading. 5.2.1 Mode I dominated interfacial cracks: Constrained simulations Fig. 30. Atomic details of nucleation of the secondary crack under tensile dominated loading. The plot shows the shear stress field σxy near the crack tip. Atoms with the energy of a free surface are drawn as larger, bluish atoms. The plot suggests that a maximum peak of the shear stress ahead of the crack tip leads to breaking of atomic bonds and creation of new crack surfaces. After the secondary crack is nucleated, it coalesces with the mother crack and moves supersonically through the material. M.J. Buehler and H. Gao 50 snapshots are taken, the crack propagates at a (slightly) super-Rayleigh speed with respect to the soft material. Since the crack motion is still subsonic, no shock front is observed. Figure 28 shows the stress field and the particle velocity field after nucleation of the secondary crack. The secondary crack propagates supersonically with respect to the soft material and the Mach cones in the right half space (soft material) are clearly visible. Figure 29 shows the potential energy field for a crack after the secondary, supersonic crack is nucleated. The two mach cones can clearly be seen. Similar mother-daughter mechanism for mode I dominated interfacial fracture has also been observed for smaller elastic mismatch ratios 2=Ξ , 5=Ξ , and 7=Ξ (however, not in all cases crack motion is supersonic with respect to the soft material since the Rayleigh- wave speed of the stiff material may be smaller than the longitudinal wave speed of the soft material). What is governing nucleation of the secondary crack? MD simulations can be particularly helpful in investigating the atomistic details of the mechanism of nucleation of secondary cracks. In Figure 30, we plot the shear stress field xyσ for several instants in time during nucleation of the secondary crack (note we chose quite small time intervals between the snapshots). Atoms with the energy of a free surface are colored blue and highlighted by larger spheres. The plot suggests that a peak of shear stress ahead of the crack tip may have caused breaking of atomic bonds. After the secondary crack is nucleated, it coalesces with the mother crack and moves supersonically with respect to the soft material. This mechanism via a peak in shear stress is quite similar to the understanding of nucleation of secondary cracks in shear loaded cracks in homogeneous systems as discussed in [38]. We note that the location of the maximum tensile stress xxσ does not coincide with the location of nucleation of the secondary crack; rather, the latter was found to correlate with a peak in shear stress ahead of the mother crack. We there conclude that there exists a peak in shear stress ahead of a mode I dominated interfacial crack moving at the Rayleigh wave speed of the soft material and this peak shear stress leads to subsequent nucleation of a secondary crack which breaks the sound barrier at the soft Raleigh speed. Modeling dynamic fracture using large-scale atomistic simulations 51 5.2.2 Unconstrained simulations In this section, we relieve the constraint that crack propagation be restricted to the interface between the stiff and the soft material (glue). We report the results of three cases: (i) Crack motion is constrained along the interface until the secondary daughter crack is nucleated, and from then on crack motion is left unconstrained (we relieve the constraint at 126=t ). (ii) Crack motion is unconstrained from the beginning of the simulation; the condition for breaking is given by 17.1break =r in both stiff and soft solid (therefore, the fracture surface energy is lower in the right half space than in the left half space). (iii) Crack motion is unconstrained from the beginning of the simulation; in contrast to case (ii), here we adjust breakr in the left material such that the fracture surface energy is equal to that in the right material so that the slab has a homogeneous fracture energy: )(1 00 right break left break rrrr − Ξ −= (28) Fig. 33. Crack dynamics at a bimaterial interface with homogeneous fracture energy everywhere in the slab (case (iii)). In this study, the fracture surface energy in both the stiff and soft materials is equal. We observe that the crack continuously branches off into the soft solid and creates numerous small cracks. This indicates a general tendency of the crack to branch off into the soft material. M.J. Buehler and H. Gao 52 and thus for 17.1=rightbreakr we assume that 13749.1= left breakr . In case (i), we observe that the crack continues to move along its original direction and does not branch off into the soft or the stiff material until the bi-crystal is completely broken. In case (ii), we observe that before the daughter crack is nucleated, the crack starts to wiggle and is thus limited in speed. Figure 29 shows the simulation snapshots of case (ii). We find that the crack branches off into the soft material and continues to propagate there. The plot on the right hand side in Figure 31 shows a zoom into the crack tip region, showing the initial crack surface roughness (A) and the large dominant branch (B) into the soft material. Fig. 34. Crack tip history for a shear loaded crack propagating along an interface with stiffness ratio 3=Ξ . This crack tip history plot reveals a mother-daughter- granddaughter mechanism: After a critical time, a secondary crack (the daughter crack) is nucleated ahead of the mother crack. Shortly afterwards, a tertiary crack (the granddaughter crack) is nucleated. The transition in crack speed is from the daughter to the granddaughter is more continuous compared to that from the mother to the daughter (see also the velocity history plot depicted in Fig. 35). Modeling dynamic fracture using large-scale atomistic simulations 55 reminiscent of the mother-daughter mechanism of mode II cracks in homogeneous materials [31, 38]. Shortly after the daughter crack is nucleated, another change in propagation speed is observed at time t ≈ 108 (see Figures 34 and 35). This is due to initiation of a tertiary granddaughter crack ahead of the daughter crack. The nucleation of the granddaughter crack occurs at a short time (Δt ≈ 8) after nucleation of the daughter crack (Figs. 34 and 35). The transitions from mother to daughter and daughter to granddaughter cracks are quite sharp, and the crack motion stabilizes after the granddaughter crack is nucleated. Figure 36 shows a few snapshots of the potential energy field from the initial configuration until the birth of the granddaughter crack. A similar sequence is shown in Figure 37, depicting the σyy stress field. Fig. 37. The plot shows the σyy stress field near a shear loaded interface crack with stiffness ratio 3=Ξ (color code: red corresponds to high stresses), for a similar sequence as shown in Figure 36. The Mach cones can clearly be seen in snapshot (4). M.J. Buehler and H. Gao 56 In Figure 38, we show a blowup plot of the potential energy field close to the crack tip, with crack surfaces highlighted with “larger” atoms in red color. The mother (A), daughter (B) and granddaughter (C) cracks can be clearly identified. The result suggests that shear dominated cracks along a bimaterial interface can reach the longitudinal wave speed of the stiffer material and can become supersonic with respect to the wave speeds of the soft material. Similar observations have been made in experiments. Rosakis et al. [96] reported crack speeds exceeding the longitudinal wave speed of the softer material and, at least on one occasion, reaching the longitudinal wave speed of the stiffer material. In another experimental study by Wu and Gupta [97], a crack at Nb-saphire interface was found to approach the longitudinal wave speed of the stiff material. Our results suggest three critical wave speeds for shear dominated interface cracks: the Rayleigh-wave speed of the soft material, the longitudinal wave Fig. 38. The plot shows the potential energy field in the vicinity of the crack tip, shortly after nucleation of the granddaughter crack (color code: red corresponds to high potential energy and blue corresponds to the potential energy of atoms in a perfect lattice). The mother, daughter and granddaughter cracks can be identified. In the blow-up on the right, the mother (A), daughter (B), and granddaughter (C) are marked. The mother and daughter cracks are represented by waves of stress concentration in the soft material. For the daughter crack (B), the Mach cone associated with the shear wave speed in the soft material can clearly be observed. The granddaughter crack (C) carries two Mach cones in the soft material. Modeling dynamic fracture using large-scale atomistic simulations 57 speed of the soft material and the longitudinal wave speed of the stiff material. Atomistic studies are particularly suitable for studying the details of the processes close to the crack tip. Figure 39 shows a time sequence of the shear stress field very close to the crack tip during the nucleation of the secondary and tertiary cracks. The results suggest that a peak (concentration) in shear stress ahead of the crack causes nucleation of the daughter and granddaughter cracks. The shear stress peak ahead of the mother crack moves with the shear wave speed of the softer material (see Fig. 39, snapshot (2)). In Fig. 39, snapshot (3), the daughter crack has just nucleated, but there appears another peak in shear stress a few atomic spacings ahead of the daughter crack moving close to the shear wave speed of the stiffer material. In Fig. 39, snapshot (4), the granddaughter crack has appeared, as indicated by two shock fronts in the softer material. These observations are, at least qualitatively, in agreement with the mother-daughter mechanism observed for mode II cracks in homogeneous materials [31, 38]. The corresponding times of the snapshots are given in the caption of Fig. 39, and can be compared with the crack tip history and crack speed history shown in Figs. 34 and 35. 5.4 Discussion and conclusions: dynamics of interfacial cracks The studies reported in this paper shows that cracking along a bimaterial interface is a rich phenomenon. We show that, barring hyperelastic effects, the limiting speed of an interfacial crack is the longitudinal wave speed of the stiffer material. A shear-dominated interface crack can propagate supersonically with respect to the softer material. We find that the crack speed changes discontinuously as the loading is increased, with more detailed analysis revealing a mother- daughter (tensile dominated loading) and mother-daughter- granddaughter mechanism (shear dominated loading) to achieve the ultimate limiting speed (see Figs. 25, 34 and 35). A clear jump in crack velocity can be observed when the daughter crack is nucleated. However, the nucleation of the granddaughter crack occurs with a more continuous velocity change. M.J. Buehler and H. Gao 60 bimaterial interfaces. Combining large-scale atomistic simulations with continuum mechanics is a powerful approach in modern studies of dynamic fracture mechanics. 6. Summary and conclusion We have presented a review of some of the recent results on large scale MD modelling of dynamic fracture. We have reported progresses in three areas of focus: (i) crack limiting speed, (ii) crack tip instabilities, and (iii) crack dynamics at interfaces. The main conclusion is that hyperelasticity, the elasticity of large deformation, is crucial in order to form a clear picture of the failure process near a rapidly moving crack tip: Both the maximum crack speed and the dynamic instability are strongly influenced by the large-strain elastic properties. Phenomena such as intersonic mode I fracture (Figure 9) and supersonic mode II fracture (Figure 10) can not at all be addressed by classical, linear elastic theories. The observation of stable cracks in homogeneous materials at speeds beyond the Rayleigh-wave speed, as shown in Figure 23, can only be understood from a hyperelastic point of view. Further, we have shown that hyperelastic softening can significantly reduce the critical instability speed in cracks, as shown for instance in Figure 19. Classical linear elastic theories such as the Yoffe- model [78] fail when significant softening is present at the crack tip. Our new concept of a characteristic energy length scale χ helps to understand the relative significance of hyperelasticity in dynamic fracture (see Figure 13), both for the limiting speed of cracks (Figure 12) and for the instability dynamics (Figure 20). Further, we find that interfaces and geometric confinement can play an important role in the dynamics of cracks. Crack propagation constrained along interfaces can significantly change the associated maximum speeds of crack motion. This is illustrated for instance by the studies using the concept of a weak fracture layer where the Rayleigh- wave speed of cracks can be attained by cracks (see Figure 6), versus the studies of cracks in homogeneous materials where the crack starts to wiggle at 73% of the theoretical limiting speed (see Figure 17). If cracks propagate along interfaces of elastically dissimilar materials, the Modeling dynamic fracture using large-scale atomistic simulations 61 maximum crack speed can significantly change and new mechanisms of crack propagation such as daughter and granddaughter cracks appear, as illustrated in Figures 24, 32 and 33, leading to supersonic fracture relative to the soft material. The results discussed in this article suggest that the definition of wave speeds according to the small-strain elastic properties is questionable in many cases, and a more complete picture of dynamic fracture should include local wave speeds near the cohesive failure of materials near a crack tip. Acknowledgments The simulations were carried out at the Garching Supercomputer Center of the Max Planck Society and at the MARS Linux cluster at the Max Planck Institute for Metals Research. We acknowledge continuing discussions with Dr. Farid F. Abraham on modeling dynamic fracture with molecular-dynamics simulations. References 1. Field, J.E., Brittle fracture-its study and application. Contemp. Phys., 1971. 12: p. 1-31. 2. Rice, J.R. and R.M. Thomson, Ductile versus brittle behavior of crystals. Phil. Mag., 1974. 29: p. 73-97. 3. Argon, A., Brittle to ductile transition in cleavage failure. Acta Metall., 1987. 35: p. 185-196. 4. Vashishta, P., R.K. Kalia, and A. Nakano, Large-scale atomistic simulations of dynamic fracture. Comp. in Science and Engrg., 1999: p. 56-65. 5. Kadau, K., T.C. Germann, and P.S. Lomdahl, Large-Scale Molecular- Dynamics Simulation of 19 Billion particles. Int. J. Mod. Phys. C, 2004. 15: p. 193. 6. Buehler, M.J., et al., The dynamical complexity of work-hardening: a large- scale molecular dynamics simulation. Acta Mechanica Sinica, 2005. 21(2): p. 103-111. 7. Zhou, S.J., et al., Large-scale molecular-dynamics simulations of three- dimensional ductile failure. Phys. Rev. Lett., 1997. 78: p. 479-482. 8. Zhou, S.J., et al., Large-scale molecular-dynamics simulations of dislocation interactions in copper. Science, 1998. 279: p. 1525-1527. M.J. Buehler and H. Gao 62 9. Buehler, M.J., et al., Atomic Plasticity: Description and Analysis of a One- Billion Atom Simulation of Ductile Materials Failure. Comp. Meth. in Appl. Mech. and Engrg., 2004. 10. Abraham, F.F., How fast can cracks move? A research adventure in materials failure using millions of atoms and big computers. Advances in Physics, 2003. 52(8): p. 727-790. 11. Bulatov, V., et al., Connecting atomistic and mesoscale simulations of crystal plasticity. Nature, 1998. 391: p. 669-672. 12. Abraham, F.F., et al., Simulating materials failure by using up to one billion atoms and the world's fastest computer: Work-hardening. P. Natl. Acad. Sci. USA, 2002. 99(9): p. 5783-5787. 13. Abraham, F.F., et al., Simulating materials failure by using up to one billion atoms and the world's fastest computer: Brittle Fracture. P. Natl. Acad. Sci. USA, 2002. 99(9): p. 5788-5792. 14. Freund, L.B., Dynamic Fracture Mechanics. 1990: Cambridge University Press, ISBN 0-521-30330-3. 15. Broberg, K.B., Cracks and Fracture. 1990: Academic Press. 16. Slepyan, L.I., Models and Phenomena in Fracture Mechanics. 2002: Springer, Berlin. 17. Cheung, K.S. and S. Yip, A molecular-dynamics simulation of crack tip extension: the brittle-to-ductile transition. Modelling Simul. Mater. Eng., 1993. 2: p. 865-892. 18. Marder, M. and J. Fineberg, How things break. Phys. Today, 1996. 49(9): p. 24-29. 19. Marder, M., Molecular Dynamics of Cracks. Computing in Science and Engineering, 1999. 1(5): p. 48-55. 20. Holland, D. and M. Marder, Ideal brittle fracture of silicon studied with molecular dynamics. Phys. Rev. Lett., 1998. 80(4): p. 746. 21. Ashurst, W.T. and W.G. Hoover, Microscopic fracture studies in 2- dimensional triangular lattice. Phys. Rev. B, 1976. 14(4): p. 1465-1473. 22. Abraham, F.F., et al., Instability dynamics of fracture: A computer simulation investigation. Phys. Rev. Lett., 1994. 73(2): p. 272-275. 23. Allen, M.P. and D.J. Tildesley, Computer Simulation of Liquids. 1989: Oxford University Press. 24. Abraham, F.F., et al., A Molecular Dynamics Investigation of Rapid Fracture Mechanics. J. Mech. Phys. Solids, 1997. 45(9): p. 1595-1619. 25. Fineberg, J., et al., Instability and dynamic fracture. Phys. Rev. Lett., 1991. 67(4): p. 457-460. 26. Marder, M. and S. Gross, Origin of crack tip instabilities. J. Mech. Phys. Solids, 1995. 43(1): p. 1-48. Modeling dynamic fracture using large-scale atomistic simulations 65 61. Buehler, M.J. and H. Gao, Biegen und Brechen im Supercomputer. Physik in unserer Zeit, 2004. 35(1). 62. Gao, H., A theory of local limiting speed in dynamic fracture. J. Mech. Phys. Solids, 1996. 44(9): p. 1453-1474. 63. Gao, H., Elastic waves in a hyperelastic solid near its plane-strain equibiaxial cohesive limit. Philosphical Magazine Letters, 1997. 76(5): p. 307-314. 64. Fineberg, J., et al., Instability in the propagation of fast cracks. Phys. Rev. B, 1992. 45(10): p. 5146-5154. 65. Sharon, E. and J. Fineberg, Confirming the continuum theory of dynamic brittle fracture for fast cracks. Nature, 1999. 397: p. 333. 66. Buehler, M.J., Atomistic simulations of brittle fracture and deformation of ultra thin copper films, in Chemistry. 2004, University of Stuttgart: Stuttgart. 67. Tsai, D.H., Virial theorem and stress calculation in molecular-dynamics. J. of Chemical Physics, 1979. 70(3): p. 1375-1382. 68. Zhou, M., Equivalent continuum for dynamically deforming atomistic particle systems. Phil. Mag. A, 2002. 82(13). 69. Zimmermann, J.A., Continuum and atomistic modeling of dislocation nucleation at crystal surface ledges. 1999, Stanford University. 70. Buehler, M.J., H. Gao, and Y. Huang, Continuum and Atomistic Studies of the Near-Crack Field of a rapidly propagating crack in a Harmonic Lattice. Theoretical and Applied Fracture Mechanics, 2004. 41: p. 21-42. 71. Buehler, M.J., A. Hartmaier, and H.J. Gao, Atomistic and continuum studies of crack-like diffusion wedges and associated dislocation mechanisms in thin films on substrates. Journal Of The Mechanics And Physics Of Solids, 2003. 51(11-12): p. 2105-2125. 72. Griffith, A.A., The phenomenon of rupture and flows in solids. Phil. Trans. Roy. Soc. A, 1920. 221: p. 163-198. 73. Broberg, K.B., Dynamic crack propagation in a layer. Int. J. Solids Struct., 1995. 32(6-7): p. 883-896. 74. Cramer, T., A. Wanner, and P. Gumbsch, Energy dissipation and path instabilities in dynamic fracture of silicon single crystals. Phys. Rev. Lett., 2000. 85: p. 788-791. 75. Cramer, T., A. Wanner, and P. Gumbsch, Crack Velocities during Dynamic Fracture of Glass and Single Crystalline Silicon. Phys. Status Solidi A, 1997. 164: p. R5. 76. Gumbsch, P., S.J. Zhou, and B.L. Holian, Molecular dynamics investigation of dynamic crack instability. Phys. Rev. B, 1997. 55: p. 3445. 77. Deegan, R.D., et al., Wavy and rough cracks in silicon. Phys. Rev. E, 2003. 67(6): p. 066209. 78. Yoffe, E.H., The moving Griffith crack. Phil. Mag., 1951. 42: p. 739-750. 79. Petersan, P.J., et al., Cracks in rubber under tension break the shear wave speed limit. Phys. Rev. Letters, 2004. 93: p. 015504. M.J. Buehler and H. Gao 66 80. Gao, H., Surface roughening and branching instabilities in dynamic fracture. J. Mech. Phys. Solids, 1993. 41(3): p. 457-486. 81. Eshelby, J.D., Elastic field of a crack extending non-uniformly under general anti-planbe loading. J. Mech. Phys. Solids, 1969. 17(3): p. 177-199. 82. Freund, L.B., Crack propagation in an elastic solid subjected to general loading, II. Nonuniform rate of extension. J. Mech. Phys. Solids, 1972. 20: p. 141-152. 83. Abraham, F.F., Unstable crack motion is predictable. Advances in Physics, 2005. 53: p. 1071-1078. 84. Heizler, S.I., D.A. Kessler, and H. Levine, Mode I fracture in a nonlinear lattice with viscoelastic forces. Phys. Rev. E, 2002. 6: p. 016126. 85. Buehler, M.J., et al., The Computational Materials Design Facility (CMDF): A powerful framework for multiparadigm multi-scale simulations. Mat. Res. Soc. Proceedings, 2006. 894: p. LL3.8. 86. Buehler, M.J., A.C.T.v. Duin, and W.A. Goddard, Multi-paradigm multi-scale modeling of dynamical crack propagation in silicon using the ReaxFF reactive force field. Submitted for publication to Phys. Rev. Lett., 2005. 87. Buehler, M.J., A.C.T.v. Duin, and W.A. Goddard, Multi-paradigm multi-scale modeling of dynamical crack propagation in silicon using the ReaxFF reactive force field. Mat. Res. Soc. Proceedings 904E, 2006: p. BB4.28. 88. Tersoff, J., Empirical interatomic potentials for carbon, with applications to amorphous carbon. Phys. Rev. Lett., 1988. 61(25): p. 2879-2883. 89. Duin, A.C.T.v., et al., ReaxFF SiO: Reactive Force Field for Silicon and Silicon Oxide Systems. J. Phys. Chem. A, 2003. 107: p. 3803-3811. 90. Buehler, M.J. and H. Gao, A mother-daughter-granddaughter mechanism of supersonic crack growth of shear dominated intersonic crack motion along interfaces of dissimilar materials. Journal of the Chinese Institute of Engineers, 2004. 27(6): p. 763-769. 91. Xu, X.P. and A. Needleman, Numerical simulations of dynamic interfacial crack growth allowing for crack growth away from the bond line. International Journal of Fracture, 1996. 74(3): p. 253-275. 92. Xu, X.P. and A. Needleman, Numerical simulations of dynamic crack growth along an interface. International Journal of Fracture, 1996. 74(4): p. 289-324. 93. Rosakis, A.J., O. Samudrala, and D. Coker, Cracks faster than the shear wave speed. Science, 1999. 284(5418): p. 1337-1340. 94. Burridge, R., Admissible speeds for plain-strain self-similar cracks with friction but lacking cohesion. Geophys. J. Roy. Astron. Soc., 1973. 35: p. 439-455. 95. Andrews, D.J., Rupture velocity of plane strain shear cracks. J. Geophys. Res., 1976. 81: p. 5679-5687. 96. Rosakis, A.J., Intersonic crack propagation in bimaterial systems. J. Mech. Phys. Solids, 1998. 6(10): p. 1789-1813. Modeling dynamic fracture using large-scale atomistic simulations 67 97. Wu, J.X. and V. Gupta, Observations of Transonic Crack Velocity at a Metal/Ceramic Interface. Journal of the Mechanics and Physics of Solids, 2000. 48(3): p. 609-619. 98. Zak, A.R. and M.L. Williams, Crack point singularities at a bi-material interface. J. Appl. Mech., 1963. 30: p. 142-143. 99. Williams, M.L., The stresses areound a fault or crack in dissimilar media. Bull. Seismol. Soc. America, 1959. 49: p. 199-204. 100. Liu, C., J. Lambros, and A.J. Rosakis, Highly transient elastodynamic crack growth in a bimaterial interface: higher order asymptotic analysis and experiments. J. Mech. Phys. Solids, 1993. 41: p. 1887-1954. 101. Liu, C., Y. Huang, and A.J. Rosakis, Shear dominated transonic interfacial crack growth in a bimaterial-II. Asymptotic fields and favorable velocity regimes. J. Mech. Phys. Solids, 1993. 43(2): p. 189-206.