1.021 , 3.021, 10.333, 22.00 I ntroduc tion to Modeling and Simulation

Spring 2011

Part I – C ontinuum and partic le me thods

Application to modeling brittle materials

Lecture 7

Markus J. Buehler

Laboratory for Atomistic and Molecular Mechanics Department of Civil and Environmental Engineering Massachusetts Institute of Technology

Content overview

I. Particle and continuum me thods

1. Atoms, molecul e s, chemistry

2. Continuum modeling approac hes and solution approaches

3. Statistical mechanics

4. Molecular dynamics, Monte Carlo

5. Visualization and data analysis

6. Mechanical proper ties – applic ation: how things fail (and how to prevent it)

7. Multi-scale modeling par adigm

8. Biological systems (simulation in biophysics) – h ow proteins work and how to model them

II. Quantum mechanical methods

1. It’s A Q uantum World: T he Theory of Quantum Mechanics

2. Quantum Mechanics: Practice Makes Perfect

3. The Many-Body Problem: Fr om Many-Body to Single- Particle

4. Quantum modeling of materials

5. From Atoms to Solids

6. Basic pr operties of mater i als

7. Advanced proper ties of materials

8. What else can we do?

Lectures 2-13

Lectures 14-26

Overview: Material covered so far…

 Lecture 1: Broad introduction to IM / S

 Lecture 2 : Introduction to atomis tic and continuum modeling (multi- scale modeling paradigm, differenc e between continuum and atomistic approach, case study: diffusion)

 Lecture 3 : Basic statistical mechan ics – p roperty calculation I (property calculation: microscopic states vs. macroscopic properties, ensembles, probability density and partition function)

 Lecture 4 : Property calculation II (Monte Carlo, advanced property calculation, introduction to chemical interactions)

 Lecture 5: How to model chemical interactions I (example: movie of copper deformation/disloc ations, etc.)

 Lecture 6: How to model chemical interactions II (pair potentials, fracture – introduction)

 Lecture 7: Application – M D simu lation of materials failure

Lecture 7: Application to modeling brittle materials

Outline:

1. Basic deformation mechanism in brittle materials - c rack extension

2. Atomistic modeling of fracture

2.1 Physical properties of atomic lattices

2.2 Applic ation

3. Bond order force fields - h ow to model chemical reactions

Goal of today’s lecture:

 Apply our tools to model a particu lar material phenomena – b rittle fracture ( useful for pset #2 )

 Learn bas ics in fracture of brittle materials

 Learn how to build a model to describ e brittle fracture (from scratch)

1. Basic deformation mechanism in brittle materials - crack extension

Brittle fracture

 Materials: glass, silicon, many ceramics, rocks

 At large loads, rather than accommodating a shape change, materials break

Imag e court e sy of quin n.anya . Li ce nse : CC - B Y .

Basic fracture process: dissipation of elastic energy

Undeformed Stretching=store elastic energy Release elastic energy

dissipated into breaking chemical bonds

Continuum description of fracture

 Fracture is a diss ipative process in which elastic energy is dissipated to break bonds (and to heat at large crack speeds)

 Energy to break bonds = surface energy  s (energy necessary to create new surface, dimens ions: energy/area, Nm/m 2 )

"Relaxed" element



"Strained" element

a ~

 a

a ~

a

(2)



(1)

Image by MIT OpenCou r seWare.





2 E

1  2 ~ ! ~

  a B

2 E

 2  s B a

change of elastic (potential) energy = G

energy to create

surfaces 8

Griffith condition for fracture initiation

 Energy release rate G , that is, the elastic energy released per unit crack advance must be equal or lar ger than the energy necessary to create new surfaces

G : 

1  2 

2 E

 2 y s

 Provides criterion to predict fail ure initiation

 Calculation of G can be complex, but straight forward for thin st rips as shown above

 Approach to calculate G based on “ stress intensity factor ” ( s e e further literature, e.g. Br oberg, Anderson, Freund, Tada)

Brittle fracture mechanisms

 Once nuc leated, cracks in brittle materials spread rapidly, on the order of sound speeds

 Sound speeds in materials (=wave speeds):

 Ray leigh-wave speed c R (speed of surface waves)

 shear wave speed c s (speed of shear waves)

 longitudinal wave speed c l (speed of longitudinal waves )

 Maximum speeds of cracks is given by sound speeds, depending on mode of loading (mode I, II, II I )

Linear elastic continuum theory

Brittle fracture loading conditions

 Commonly cons ider a single crack in a material geometry, under three types of loading: mode I, mode II and mode III

Mode I

Mode II

Mode III

T ensile load, focus of this lectur e

Image by MIT OpenCou r seWare.

Limiting speeds of cracks: linear elastic

continuum theory

Mother-daughter mechanism

Super- Sub-Rayleigh Rayleigh

Mode II

Intersonic

Supersonic

Mode I

C r C s

C l

Limiting speed v

Subsonic

Supersonic

Mode III

C s

C l

Limiting speed v

Linear Nonlinear

Image by MIT OpenCou r seWare.

• C racks can not exceed the limiting speed given by the corresponding wave speeds unless materia l behavior is nonlinear

• C racks that exceed limiting s peed would produce energy (physically impos sible - linear elastic continuum theory )

Sound speeds in materials: overview

Wave speeds are calc ulated based on elastic properties of material

 = shear modulus

E  8 / 3 

  3 / 8 E

Brittle fracture mechanisms: fracture is a multi- scale phenomenon, from nano to macro

Imag e removed du e t o copyrigh t rest r i c t ion s. Fig. 6.2 in Buehler, Ma r kus J. At omisti c Modeling of Ma t e rials Failure . Spring er , 2008.

Physical reason for crack limiting speed

 Phys ical (mathematical) reason for th e limiting speed is that it becomes increasingly difficult to increase the speed of the crack by adding a larger load

 When the crack approaches the limiting speed, the resistance to fracture diverges to infinity (= dynamic fracture toughness )

divergence

Driving force (applied load)

1  2

Imag e removed du e t o co pyrigh t restric t ion s.

Pleas e s e e: Fig. 6.15 in Bu eh ler, Markus J. Atomistic Modeling of Materials Fa ilure . Sprin ger, 2008.

G   

2 E

f ( y s , v )

 2

G ~

E 15

2. Atomistic modeling of fracture

16

What is a model?

Mike Ashby (Cambridge University):

 A model is an idealization. Its relationship to the real problem i s like that of the map of the London tube trains to the real tube systems: a gross simplification, but one that captures certain essentials.

“Physical situation” “Model”

© G oo gle, I n c. All righ t s res e rved. This c o n t en t is excluded fr om our Creat i ve Co mmo ns li cense . Fo r m o re inf o rma t io n , see h ttp:/ /ocw.mi t . e du/fai ruse .

© M a s s a c h u s e tts B a y T r a n s p o r t a ti o n Au t h o r i ty . A l l r i g h ts res e rved. This c o nt ent is exclud ed from our Creat i ve Commons lic e n s e. For more infor m at ion, s ee http:/ /ocw. m it.edu/fai ruse . 17

A “simple” atomistic model: geometry

  1 k

2 0

( r  r 0 )

stable configuration of pair potential

(store energy)

2D triangular (hexagonal)

Y

X

I x

a

r

3 1/2 r

o

W eak fracture layer

I y

elasticity

lattice

Pair potential to describe atomic interactions Confine crack to a 1D path (weak fracture layer):

Define a pair potential

whose bonds never break (bulk) and a potential whose bonds break

Image by MIT OpenCou r seWare.

Harmonic and harmonic bond snapping potential



  1 k ( r  r ) 2

2 0 0



 1

 2

   1

0 ( r  r 0 )

r  r break



 2 0

( r break

 r 0 )

r  r break

2.1 Physical properties of atomic lattices

How to calculate elastic properties and fracture surface energy – parameters to link with continuum theory of fracture

  (  )

free energy dens ity (energy per unit volume)

 2  (  )

Stress    

Young’s modulus

E    2

  E 

1D example – “Cauchy-Born rule”

 Impose homogeneous strain field on 1D string of atoms

 Then obtain

  E 

from that

r  ( 1   )  r 0

 (  ) 

1

r 0  D

 ( r ) 

1

r 0  D

 ( ( 1   )  r 0 )

Strain energy density function

interatomic potential

    (  )

 

 2  (  )

E    2

out-of-plane area

r 0  D

atomic volume 21

2D hexagonal lattice

[010]

r 2 r 2

r 2 r 2

[100]

r 1

r 1

r 1

r 1

y

r 2 r 2

r 2 r 2

x

 ' '  k

Image by MIT OpenCou r seWare.

 (  ij ) 

 ' '  3  2

8

 2 

11  22 

2  ( 

  12

) 2 

    (  ij )

 2  (  )

c  ij

ij

ij

ijkl

  ij   kl

22

2D triangular lattice, LJ potential

[010]

r 2 r 2

r 2 r 2

[100]

r 1

r 1

r 1

r 1

y

r 2 r 2

r 2 r 2

x

Image by MIT OpenCou r seWare.

 = 1

 = 1

12:6 LJ potential

Image by MIT OpenCourseWare. 23

2D triangular lattice, harmonic potential

Str ess vs. Strain and Poisson Ratio

4

50

3 45

T angent Modulus

40

2

35

1

30

0

0 0.02 0.04

Strain

0.06 0.08

25

0 0.02 0.04

Strain

0.06 0.08

Elastic properties of the triangular lattice with harmonic interactions, stress versus strain (left) and

tangent moduli E x and E y (right). The stress state is uniaxial tension, that is the stress in the direction orthogonal to the loading is relaxed and zero.

Image by MIT OpenCou r seWare. 24

Elastic properties

Enables to calc ulate wav e speeds:

 ' '  k

Surface energy calculation

y

fracture plane

x

fracture plane

r o

3 1/2 r

o

Image by MIT OpenCou r seWare.

Harmonic potential with bond snapping distance r break

Note: out-of-plane unity thickness

2.2 Application

Focus: effects of material nonlinearities (reflected in choice of model)

Coordinate system and atomistic model

Y

X

I x

a

r

3 1/2 r

o

W eak fracture layer

I y

Image by MIT OpenCou r seWare.

Pair potential to describe atomic interactions Confine crack to a 1D path (weak fracture layer)

Linear versus nonlinear elasticity=hyperelasticity

Hyperela sticity

Strain

Image by MIT OpenCou r seWare.

Linear elasticity: Young’s modulus (stiffness) does not change with deformation

Nonlinear elasticity = hyperelasticity : Young’s m odulus (stiffn e ss) changes with deformation

Subsonic and supersonic fracture

 Under certain conditions, materi al nonlinearities (that is, the behavior of materials under large deformation = hyperelasticity) becomes important

 This can lead to different limiting speeds

than described by the model introduced above

 ( r ) ~ 1

Strain

Deformation field near a crack

small deformation

Image by MIT OpenCou r seWare.

Limiting speeds of cracks

Image by MIT OpenCou r seWare.

• U nder presence of hyperelastic effects, cracks can exceed the conventional barrier given by the wave speeds

• T his is a “local” e ffect due to enhancement of energy flux

• Subsonic fracture due to loc a l soft ening, that is, reduction of energy flux

31

Stiffening vs. softening behavior

real materials

Hyperela sticity

Strain

“linear elasticity”



2 E

Increased/decreased wav e speed

Image by MIT OpenCou r seWare.

MD model development: biharmonic potential

Metals

Polymers, rubber

• Stiffness change under deformation , with different strength

• A tomic bonds break at critical atomic separation

• Want: simple set of parameters that c ontrol these properties (as few as possible, to gain generic insight)

Biharmonic potential – c ontrol parameters

 ' k 1

k 0

k ratio

 k 1

/ k 0

 break

  r  r 0

r

r on r break r 0

Biharmonic potential definition

The biharmonic potential is defined as:

where is the critical atomic separ ation for the onset of th e hyperelastic effec t, and

and

are found by continuity conditions of the potential at .

The values and refer to the small- and large-strain spring constants.

Energy filtering: visualization approach

 Only plot atoms associated with higher energy

 Enables determination of crack tip (atoms at surface have higher energy , as they have fewer nei ghbors than atoms in the bulk)

N

En erg y of atom i

U i 

  ( r ij )

200

180

160

140

120

100

80

60

40

20

0

a

0 50 100 150 200 250 300 350 400

x- coordinate

j  1

Image by MIT OpenCou r seWare. 36

MD simulation results: confirms linear continuum theory

a ( t )

crack speed

=time derivative of a

v ( t ) =da/dt

Image by MIT OpenCou r seWare.

Virial stress field around a crack

Imag es removed du e t o copyrigh t restric t ion s.

Pleas e s e e: Fig. 6.30 in Bu eh ler, Markus J. Atomistic Modeling of Materials Fa ilure . Sprin ger, 2008.

large stresses at crack tip

- i nduce bond fail ure

Biharmonic potential – bilinear elasticity

Softening Stiffening

2 2

1.5

E t =66

1.5

1

0.5

E t =33

1

0.5

E t =66

E t =33

0

0 0.01 0.02 0.03

Str ain (  )

0

0 0.01 0.02 0.03

Str ain (  )

Image by MIT OpenCou r seWare.

Harmonic system

Subsonic fracture

 Softening material behav ior leads to subsonic fracture, that is, t he crack can never attain its theoretical limiting speed

 Materials: metals, ceramics

< 3.4 = c R

Supersonic fracture

 Stiffening material behav ior leads to subsonic fracture, that is, the crack can exceed its theoretical limiting speed

 Materials: polymers

> 3.4 = c R

Different ratios of spring constants

Reduc ed crack speed c / c R

1. 5

1. 3

Stiffening

Intersonic fracture

k 2 / k 1= 1/ 4 k 2 / k 1= 1/ 2 k2/ k 1 = 2 k2/ k 1 = 3 k2/ k 1 = 4

1. 1

Shear wave speed

0. 9

0. 7

Sub-Rayleigh fracture

Softening

1. 12 1. 13 1. 14 1. 15 1. 16 1. 17 1. 18

r on

1. 19 1. 2

Stiffenin g an d softenin g effect : Increas e o r reductio n o f crac k s pee d

Physical basis for subsonic/supersonic fracture

 Changes in energy flow at the crac k tip due to changes in local wave speed (energy flux higher in materials with higher wave speed)

 Controlled by a characteristic length scale 

R e print e d by permis sio n from Ma cmill an Publish e rs L t d: Nature.

S o urc e : Bueh ler, M., F. Abraham, and H. Gao. "H ypere l a s t i c i ty Governs Dynamic Fractur e at a Critica l L e ngth S ca l e." Nature 42 6 (2003): 1 41-6. © 20 03.

Supersonic fracture: mode II (shear)

Pleas e s e e: Bu eh ler, Markus J., Farid F. Abraham, and Huaj ian Ga o. "Hypere l as t i c i t y Governs Dynamic

Fractur e at a Critica l L e ngth S ca l e." Nature 42 6 (November 13, 2003): 141-6. 45

Theoretical concept: energy flux reduction/enhancement

K dominance zone

Chara cteristic energy length

Hyperelastic zone

Fra cture process zone

Classical

New

Energy Flux Related to Wave Speed: High local w a v e speed, high energy flux, cr ack can mo v e faster (and rev erse f or low local w a v e speed).

Image by MIT OpenCou r seWare.

Energy flux related to wave speed: high local wave speed, high energy flux, crack can move faster (and re verse for low local wave speed)

3. Bond order force fields - h ow to model chemical reactions

Challenge: chemical reactions

sp3 sp2

Energy

 stretch

C-C distance

r

CHA R MM-type potential can not describe chemical reactions

Why can not model chemical reactions with

|spring-like potentials?

 stretch

 1 k

2

stretch

( r  r 0 )

Set of parameters only valid for particular molec ule type / type of chemical bond

 bend

 1 k

2

bend ( 

  0 )

k  k

stretch, sp 2 stretch, sp 3

Reactive potentials or reactive force fields overcome these limitations

Key features of reactive potentials

 How can one accurately describe the transition energies during chemic al reactions?

 Us e computationally more efficient descriptions than relying on purely quantum mechanical (QM) methods (see part II, methods limited to 100 atoms )

H  C

 C 

 H 2

H  H 2  C

 C  H 2

inv o lves processes with electrons

A

A- -B

B

Key features of reactive potentials

 Molecular model that is capable o f describing chemical reactions

 Continuous energy landscape during reactions ( k ey to enable integration of equations)

 No typing necessary, that is, atoms can be sp, sp2, sp3… w/o further “tags” – only element types

 C o m p utation ally efficient (that is, should involve finite range interactions), so that large sys tems can be treated (> 10,000 atoms)

 Parameters with physical meaning ( such as for the LJ potential )

Theoretical basis: bond order potential

Concept: Use pair potential that depen ds on atomic environment (similar to EAM, here app lied to covalent bonds)