**UFRJ**

# MODELING DYNAMIC FRACTURE USING LARGE-SCALE atomistic simulations

(Parte **1** de 6)

CHAPTER 1

Markus J. Buehler

Massachusetts Institute of Technology,

Department of Civil and Environmental Engineering 7 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-misthackle” 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 [2]. Abraham and coworkers reported molecular-dynamics simulations of fracture in systems up to 500,0 atoms, which was a significant number at the time. In these atomistic calculations, a Lennard- Jones (LJ) potential [23] was used. The results in [2, 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 3

to numerous further studies in the following years (e.g. [26]) |

perhaps severe plastic deformation near the macroscopic crack tip. Such phenomena were observed at similar velocities in both experiments and modeling [25]. Since the molecular-dynamics simulations are performed in atomically perfect lattices, it was concluded that the dynamic instabilities are a universal property of cracks, which have been subject

The last few years have witnessed ultra large-scale atomistic modeling of dynamical fracture with system sizes exceeding one billion atoms [4, 5, 12, 13, 27, 28]. Many aspects of fracture have been investigated, including crack limiting speed [10, 29-31], dynamic fracture toughness [32] and dislocation emission processes at crack tips and during nanoindentation [3, 34]. Recent progresses also include systematic atomistic-continuum studies of fracture [29-31, 35-38], investigations of the role of hyperelasticity in dynamic fracture [30, 39] and the instability dynamics of cracks [2, 24, 39]. A variety of numerical models have been proposed, including concurrent multi-scale schemes that combine atomistic and continuum domains within a single model [40-48].

1.2 Outline of this chapter

physical aspects of dynamic fracture |

In this chapter, we mainly focus on the work involved by the authors in using simplistic interatomic potentials to probe crack dynamics in model materials, with an aim to gain broad insights into fundamental,

A particular focus of our studies is on understanding the effect of hyperelasticity on dynamic fracture. Most existing theories of fracture assume a linear elastic stress–strain law. However, the relation between stress and strain in real solids is strongly nonlinear due to large deformations near a moving crack tip, a phenomenon referred to here as hyperelasticity. Our studies strongly suggest that hyperelasticity, in contrast to most of the classical linear theories of fracture, indeed has a major impact on crack dynamics.

M.J. Buehler and H. Gao 4

The plan of this chapter is as follows. First, we review atomistic modeling techniques, in particular our approach of using simple model potentials to study dynamic fracture. We will then cover three topics: (i) confined crack dynamics along weak layers in homogeneous materials, focusing on the crack limiting speed (Section 3), (i) instability dynamics of fracture, focusing on the critical speed for the onset of crack instability (Section 4), and (i) dynamics of cracks at interfaces of dissimilar materials (Section 5). Whereas cracks are confined to propagate along a

Fig. 1. Subplot (a): A schematic illustration of the simulation geometry used in our large-scale atomistic studies of fracture. The geometry is characterized by the slab width xl and slab length yl, and the initial crack length a. We consider different loading cases, including mode I, mode I (indicated in the figure) and mode I (not shown). Subplots (b) and (c) illustrate two different possible crack orientations. Configuration (b) is used for the studies discussed in Sections 3 and 5, whereas configuration (c) is used for the studies described in Section 4. The orientation shown in subplot (c) has lower fracture surface energy than the orientation shown in subplot (b).

Modeling dynamic fracture using large-scale atomistic simulations 5

outlook to future research in this area |

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

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].

on the order of a few femtoseconds |

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Δ

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- 5]). 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

where 0rdenotes 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

4)( r 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 62/2/a=σ.

The reduced temperature is ε/kTB 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

2, 24, 31, 40, 56-58] |

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,

2.3 Model geometry and simulation procedure

We typically consider a crack in a two-dimensional geometry with slab size yxll× and initial crack length a, as schematically shown in Fig. 1. The slab size is chosen large enough such that waves reflected

between 2 and 4 times larger than the size orthogonal to the initial crack |

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

of our studies are carried out in a two-dimensional, hexagonal lattice |

periodic boundary conditions in the out-of-plane z direction |

The slab is initialized at zero temperature prior to simulation. Most The three-dimensional studies are performed in a FCC lattice, with 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 I 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 I), 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]

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 yposition 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 8

performed within a region of constant stress intensity factor |

intervals to eliminate very high frequency fluctuations. To obtain the steady state velocity of the crack, crack speed measurements are

3. Constrained cracks in homogeneous materials: How fast can cracks propagate?

In this section we focus on the limiting speed of cracks, discussing recent results on cracks faster than both shear and longitudinal wave speeds in an elastically stiffening material [10, 12, 30, 31].

The elasticity of a solid clearly depends on its state of deformation.

Metals will weaken, or soften, and polymers may stiffen as the strain approaches the state of materials failure. It is only for infinitesimal deformation that the elastic moduli can be considered constant and the elasticity of the solid linear. However, many existing theories model fracture using linear elasticity. Certainly, this can be considered questionable since material is failing at the tip of a dynamic crack

(Parte **1** de 6)