Nanoscale oscillatory fracture propagation in metallic glasses

Physica A 388 (2009) 1978–1984

Contents lists available at ScienceDirect

Physica A journal homepage: www.elsevier.com/locate/physa

Nanoscale oscillatory fracture propagation in metallic glasses Y. Braiman a,b,∗ , T. Egami a,c,d a

Oak Ridge National Laboratory, Oak Ridge, TN 37831, United States

b

Department of Mechanical, Aerospace, and Biomedical Engineering, University of Tennessee, Knoxville, TN 37996, United States

c

Department of Materials Science and Engineering, University of Tennessee, Knoxville, TN 37996, United States

d

Department of Physics and Astronomy, University of Tennessee, Knoxville, TN 37996, United States

article

info

Article history: Received 3 July 2008 Received in revised form 24 October 2008 Available online 20 January 2009

a b s t r a c t We describe the oscillatory crack propagation for small propagation velocities at the atomistic scale that was recently observed for brittle metallic glasses [G. Wang, Y.T. Wang, Y.H. Liu, M.X. Pan, D.Q. Zhao, W.H. Wang, Appl. Lett. 89 (2006) 121909; G. Wang, D.Q. Zhao, H.Y. Bai, M.X. Pan, A.L. Xia, B.S. Han, X.K. Xi, Y. Wu, W.H. Wang, Phys. Rev. Lett. 98 (2007) 235501]. Based on a simple model of crack propagation [Y. Braiman, T. Egami, Phys. Rev. E, 77 (2008) 065101(R)], we derived and analyzed expressions for the feature size, oscillation period, and maximum strain accumulated in the material. © 2009 Elsevier B.V. All rights reserved.

Glasses represent a good testing ground for a theory of fracture because of the absence of microstructures, such as grain boundaries, elastic anisotropy and lattice dislocations, which complicate the analysis of crystalline materials. Metallic glasses are particularly interesting, since they can be ductile as well as brittle, while most oxide glasses are brittle. Even when metallic glasses are macroscopically brittle, at the atomic scale their behavior is still controlled by viscoelasticity, while oxide glasses are atomically brittle. Surprisingly, nano-scale periodic morphology was recently discovered on the fracture surface of a brittle metallic glass [1]. To initiate fracture, a small seed crack was introduced at the edge of the glassy alloy, and the sample was fractured by a tensile stress resulting in crack propagation. The characteristic feature size observed on the fracture surface varied from 30 to approximately 200 nm. Understanding crack propagation in materials constitutes an interesting and very challenging problem. A minimal model to describe oscillatory fracture propagation at the mesoscopic scale was proposed by Blumenfeld [2]. In Blumenfeld’s model [2] the velocity of crack oscillates from some minimum (nonzero) to a maximum value for an oscillatory behavior and is constant for a steady state case. Velocity oscillations in rapid fracture in thin brittle gels were recently observed by Livne et al. [3]. Models of crack propagation were recently proposed by several groups [4–9]. However, none of these models addresses the oscillatory crack propagation at the atomistic to nanometer scale. Recently we proposed a model that characterizes time dependence of crack propagation and accounts for the tip dynamics in response to external uniform driving with a constant velocity [10]. We adopted a simple one-dimensional model for simplicity, while the real structure can be more complex. Our main purpose was to demonstrate that the periodicity in crack propagation is strongly correlated to the ratio of the elastic-to-viscous properties and provide a path for experimental confirmation of our results. The details of the model are presented in the Ref. [10]. Here we only briefly summarize the main features of the model. We consider crack propagation under a tensile stress, the so-called Mode One fracture. Until the applied stress reaches a critical value the system responds elastically. Once the critical value is exceeded the crack starts to propagate. We describe the motion of an infinitesimally small volume of material in the vicinity of the tip. As the tip moves, it deforms the vicinity by creating a deformation zone φ (see Fig. 1 and Refs. [11,12]) that is time-dependent [12]. We derived [10] the following

∗

Corresponding address: Oak Ridge National Laboratory, 1 Bethel Valley Rd, P.O. Box 2008, MS 6015, Oak Ridge, TN 37831, United States. E-mail address: [email protected] (Y. Braiman).

0378-4371/$ – see front matter © 2009 Elsevier B.V. All rights reserved. doi:10.1016/j.physa.2008.12.072

Y. Braiman, T. Egami / Physica A 388 (2009) 1978–1984

1979

Fig. 1. A schematic drawing of the tip created deformation zone.

equation to describe the motion of the infinitesimal mass in the direction of the crack propagation:

Γ x˙ = kel (ν0 t − x) − Fb

(1)

where Γ x˙ describes energy dissipation due to crack propagation, ν0 is the average crack propagation velocity that can be estimated from νt = ν0 tan θ , νt is the crack opening velocity and θ is the opening angle. The resistance force Fb can be written as: Fb = kel δ(1 − e−φ/δ (1 + φ/δ)), here δ is the characteristic length of the viscous flow zone, δ ≈ O(5–10 nm) involving approximately 20 atoms, corresponding to the width of the shear band [13] and kel is the elastic constant of the material. The proposed form of the resistance force Fb is consistent with the expression for separation potential suggested by Xu and Needelman [14] and with assumptions of the elastic chain models [15]. The expression for Fb has some similarities with the expression for the normal and shear work separation between the interfacial surfaces [16,17]. We do not take into consideration the effects of nonlinearity (hyper-plasticity) recently suggested by Buehler et al. [18] that may not be significant for brittle materials on the nanometer scale. Substituting the expression for Fb into Eq. (1) leads to:

Γ x˙ = kel (ν0 t − x) − F0 (1 − e−φ/δ (1 + φ/δ)),

(2)

where F0 ≡ kel δ . The effect of viscosity is taken in consideration through the equation for the deformation zone (DZ) φ and is related to a Kelvin–Voigt model for viscoelasticity. The deformation zone is related to the volume (per unit length) where the viscoelastic flow process to resist the tip motion is taking place. In our model, the size of the DZ φ is the diameter of the cross section of the deformation volume (see Fig. 1) and has the units of the length. The size of DZ is time-dependent and its evolution depends on the elastic and viscous properties of the material and on the average crack propagation velocity. Experimentally, time evolution of the DZ for the nanoscale plastic deformation and fracture of the polymers was studied by the in situ nanoindentation in a transmission electron microscope [12]. The deformation zone size as a function of time for plastic deformation was experimentally measured (see Fig. 2 in Ref. [12]). We derived the equation for φ based on the φ following considerations [10]. We presumed a generic form φ˙ = − τ + ν0 and τ1 = τ1 + τ 1 . Here τm is the characteristic φ

φ

relaxation time of the material expressed by the Maxwell relaxation time, τm =

η

G

m

tip

, where η is the viscosity, G is the shear

modulus. τtip is the relaxation time due to the tip motion and can be approximated as τ −1 tip ≈ 1

|˙x| a

, where |˙x| is the velocity

of the tip and a is the characteristic interatomic length of the material. Defining B ≡ τ , we obtained the following equation m for φ :

  |˙x| ˙φ = − B + φ + ν0 . a

(3)

It is interesting to note that the form of the Eq. (3) can be related to the Kelvin–Vogt model for viscoelastic fluids with B corresponding to the ratio of the shear modulus G to the viscosity coefficient η. A similar form of equations has been used to describe dry friction [[19,20] and references therein], although there are also significant differences (see, for example Ref. [19] and Eqs. (2) and (3) in our model). Moreover, for fracture propagation the mass of the tip is negligible while the effect of inertia on frictional dynamics is often considerable. We estimated the range of the parameters in the Ref. [10]. Here we assume vt ≈ O(10−7 ) m/s, θ ≈ 1, thus v0 ≈ O(10−5 ) m/s. The value of the lattice constant a is approximated as a ≈ 2.5 × 10−10 m and we assume δ ≈ O(5 × 10−9 ) m.

1980

Y. Braiman, T. Egami / Physica A 388 (2009) 1978–1984

Fig. 2. A typical behavior of the tip ξ for variety of β values. The parameter values are: ξ˙0 = 1.657 × 10−4 ; δ0 = 100.6; D = 1.04 × 10−3 .

The value of kel can be estimated as kel ≈ E δ , E is the Young modulus. Assuming E ≈ O(5 × 1010 ) Pa, we obtain kel ≈ O(2.5×102 ) N/m. We estimate kcr (due to plastic deformation round the tip) as kcr = κ kel where κ is the stress concentration factor. We assume κ ≈ 60 thus kcr ≈ O(1.5 × 104 ) N/m. We then estimate ξ0 = kel /F0 = 1/κδ ≈ O(4. × 106 ) m−1 , the characteristic length 1/ξ0 ≈ O(0.25 × 10−6 ) m, and the value of Γ can be estimated as: Γ = E ε 2 δ/2νcr , here ε is the strain far from the crack and νcr is the maximum crack propagation velocity [21] that we estimate to be of the same order of magnitude as the sound velocity (νc ≈ O(102 m/s)). We also assumed that the elastic energy released is much higher than the energy needed to create a new fracture surface [22]. We assume ε ≈ 2 × 10−2 and vcr ≈ 2 × 102 m/s thus Γ ≈ 2.5 × 10−4 J s/m2 . We now rewrite the Eqs. (2) and (3) in the dimensionless units:

ξ˙ = (ξ˙0 τ − ξ ) − (1 − (1 + ϕ/δ0 )e−ϕ/δ0 )   |ξ˙ | ϕ˙ = 1 − β + ϕ

(4) (5)

D

where ξ = ξ0 x τ = ω0 t; ξ0 = ξ˙

ξ ν

kel F0

ω0 =

kel

Γ

ξ

; ξ˙0 = ω0 ν0 = 0

Γ F0

ν0 ; ϕ =

ω0 φ ; δ0 ν0

=

ω0 δ; ν0

β =

B

ω0

=

BΓ kel

D = aξ0 =

akel ; F0

and β0 = 0B 0 = F elB ν0 . 0 A typical behavior of the crack tip is demonstrated in the Fig. 2. A typical behavior of φ is shown in the Fig. 3. When the deformation rate is low, the crack tip motion corresponds to a stick-slip motion (see Fig. 2). Thus, we first look at the ‘‘static’’ part of the solution (dξ /dt = 0). Eq. (4) then takes the following form: k

ξ = ξ˙0 τ − (1 − e−ϕ/δ0 (1 + ϕ/δ0 )).

(6)

Y. Braiman, T. Egami / Physica A 388 (2009) 1978–1984

1981

Fig. 3. A typical behavior of the response function φ for variety of β values. The other parameters values are as in the Fig. 1 and ω0 /ν0 = 2.51 × 1010 m−1 .

Substituting for ϕ in the Eq. (5) the expression for ξ derivative obtained by the direct differentiation of the Eq. (6) we obtain: 1+

ϕ˙ =

˙

 −β ϕ

ξ0 D

1+

(7)

ϕ 2 −ϕ/δ0 e δ02 D

and the time required for the tip to complete one period of motion is given by: ϕ 2 −ϕ/δ0 e δ02 D

1+

ϕmax

Z T =

1+

0

˙

ξ0 D

 dϕ. −β ϕ

(8)

The integral in the Eq. (8) can be further rewritten as: ϕmax

Z

dϕ

T = 1+

0

Define g ≡

ξ˙0 D

ξ0 D

+



−β ϕ

ϕmax

Dδ02

0

ϕ 2 e−ϕ/δ0 ˙  dϕ. ξ 1 + D0 − β ϕ

(9)

− β . Then the Eq. (9) will read as:

ϕmax

Z

˙

Z

1

dϕ

T =

1 + gϕ

0

Dδ

ϕmax

Z

1

+

2 0

0

ϕ 2 e−ϕ/δ0 dϕ. 1 + gϕ

(10)

This integral can be split into four more simple integrals as: T = T1 + T2 + T3 + T4 , where ϕmax

Z T1 =

1 + gϕ

0

T4 =

dϕ ϕmax

Z

1 Dg δ

2 0

e

;

T2 =

− δϕ

0

dϕ

1 + gϕ

0

ϕmax

Z

1 Dg δ

2 0

ϕe

− δϕ

0

dϕ;

T3 = −

0

ϕmax

Z

1 Dg δ

2 0

e

− δϕ

0

dϕ;

0

. 

− ϕmax

ϕ



1 One can show that in a good approximation T ≈ ξ˙ −β ) . 1 − e δ0 (1 + max δ0 D 0 The condition for ϕmax follows from dξ /dt = 0 thus ϕ˙ max = 1 − βϕmax . Combining this equation with the Eq. (6) we obtain:

ϕ˙ max =

1+

˙

1+

ξ0 D

 − β ϕmax

2 ϕmax e−ϕmax /δ0 δ02 D

ϕ˙ max = 1 − βϕmax .

(11)

1982

Y. Braiman, T. Egami / Physica A 388 (2009) 1978–1984

The solution of the Eq. (9) is:

  ϕmax ϕmax = ln(1 − βϕmax ) + ln − ln(ξ˙0 δ0 ). δ0 δ0

(12)

Eqs. (8) and (12) provide an accurate description of the tip periodicity T and the strain ϕ/δ0 developed in the material provided the stick-slip motion is stable. We now provide an approximationsolutionfor the case of low strain when the ratio of ϕ to the characteristics length of ϕ δ0

viscous flow δ0 is much less then one

 1 . This condition may occur close to the transition to smooth propagation. In

that case the Eq. (4) can be written as:

ξ˙ = (ξ˙0 τ − ξ ) − (1 − e−ϕ/δ0 ).

(13)

The Eq. (5) will not be changed. Eq. (6) will read as:

ξ = ξ˙0 τ − (1 − e−ϕ/δ0 ).

(14)

The Eq. (7) will read as: 1+

ϕ˙ =

˙

ξ0 D

 −β ϕ

ϕ

1 + δ D e−ϕ/δ0 0

.

(15)

And the equation for the propagation period, Eq. (8) will read as: ϕ

ϕmax

Z

1 + δ D e−ϕ/δ0 0

T =

˙

ξ0

1+

0

D

 dϕ. −β ϕ

(16)

We first analyze the case of β = 0. This expression for the period T will read as: ϕmax

Z T =

ϕ

1 + δ D e−ϕ/δ0 0 1+

0

ξ˙0 D

ϕ

dϕ.

(17)

Since the major contribution to the integral (17) comes from the values of ϕ ≈ ϕmax , we can approximate ϕ −ϕ/δ0 e  1. Therefore, δ D

ξ˙0 ϕmax D

 1 and

0

ϕmax

Z T = 0

ϕ −ϕ/δ0 e δ0 D ξ˙0 D

ϕ

dϕ =

1

ξ˙0

− δ0 .

(18)

And ξ˙0

ϕ˙ =

ϕ

D

δ0 D

ϕ

= ξ˙0 δ0 e

e−ϕ/δ0

− δϕ

0

thus ϕ ≈ −δ0 ln(1 − ξ˙0 t )

(19)

where t < T and (1 − ξ˙0 t ) > 0. ϕ We now examine a more general case of δ  1 and β > 0. 0 The Eq. (11) will read as:

ϕ˙ max = ϕ˙ max

1+

˙

ξ0 D

 − β ϕmax

ϕ

1 + δmax e−ϕmax /δ0 D = 1 − βϕmax .

(20)

Leading to:

ϕmax = ln(1 − βϕmax ) − ln(ξ˙0 δ0 ). δ0

(21)

We now investigate few asymptotic cases. First we assume βϕmax  1. Then ln(1 − βϕmax ) ≈ −βϕmax and thus (ξ0 δ0 ) ϕmax ≈ − δ01ln+βδ . The oscillation period: 0 Z ϕmax 1 + ϕ e−ϕ/δ0 Z  1 ϕmax −ϕ/δ0 δ0  − ϕmax δ0 D δ0 d ϕ ≈ T = e d ϕ = 1 − e . ξ˙ A 0 A 0 1 + D0 ϕ ˙

Y. Braiman, T. Egami / Physica A 388 (2009) 1978–1984

1983

Fig. 4. Numerical (red dots) and analytical (Eq. (12), blue line) curves for ϕmax as a functions of β in the logarithmic and linear (the Inset) scales. The parameters are the same as in the Figs. 1 and 2. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

˙

 − β . Inserting the expression for ϕmax we obtain: D     ln(ξ˙0 δ0 ) ln(ξ˙0 δ0 ) 1 δ0 1 − e 1+βδ0 ≈ T = 1 − e 1+βδ0 . A ξ˙0

Here A ≡ δ0 D

ξ0

Next we investigate the case of βϕmax ≈ O(1). ln(1 − βϕmax )  ϕmax

therefore ϕmax ≈

1 − ξ˙0 δ0

β

.

The length of the feature size can be estimated as the length of the jump during the slip event ∆x = ν0 T where the value of the period T is given by the Eq. (8). Higher values of β correspond to smaller viscosity (higher temperature) while lower values of β correspond to higher viscosity (lower temperature). In the Fig. 4 we show numerical (red dots) and analytical (Eq. (12), blue line) curves for ϕmax as a function of β in the logarithmic and linear (the Inset) scales. Stability of the solutions has not been taken into consideration in the Fig. 4 and will be considered later on. For a given set of parameters, our numerical results show an excellent agreement with the analytical results. An increase in β results in decrease of ϕmax and the amount of strain that can be accumulated in glass goes down. In Fig. 5 we show numerical (red dots) and analytical (Eq. (8), blue line) curve for T as a function of β in the logarithmic and linear (the inset) scales. Stability of the solutions has not been taken into consideration in Fig. 5 and is considered in the Ref. [10]. For a given set of parameters, our numerical results show an excellent agreement with the analytical results. An increase in β leads to decrease in the period T and for large enough β transition from the stick-slip to a continuous character of the motion of the tip occurs. A typical feature size for a given values of the parameters depends on the value of β and is the maximum for β = 0. We can estimate the feature size as ∆x ≈ ∆ξ /ξ0 . As for β = 0, ∆ξ ≈ 1 we estimate ∆xmax ≈ 1/ξ0 ≈ 2.4 × 10−7 m. The feature size ∆x = T ν0 is the maximum for β = 0 that corresponds to the limit of a very high viscosity (low temperature limit). It is then decreases with increasing β and discountinuously drops to zero at some value of β = βc that corresponds to the high temperature limit. For given values of the parameters, the smallest feature size corresponds approximately to ∆xmin ≈ δ/5 = 8 × 10−10 m. We have performed linear stability analysis of the Eqs. (4) and (5) (see Ref. [10]). As β decreases,the smooth fracture propagation where the instantaneous crack propagation velocity x˙ is equal to the external velocity ν0 (and x(t ) is a linear function of time) becomes unstable and being replaced by a jerky fracture 1

propagation. Stability analysis leads to the following equation for the transition point value βc : A2 δ 2 D(1 + A) = e− Aδ , A ≡ βc + ν/D resulting in βc ≈ 0.11. This value for βc is an excellent agreement with the numerically calculated transition point. In conclusion, we have proposed a simple model to account for the oscillatory crack propagation at the nanoscale observed for brittle metallic glasses based upon the nonlinear response of a viscoelastic medium. Our model is based upon the non-linear viscoelasticity of metallic glasses, and may have wider implications on the fracture behavior in viscoelastic media in general. The results may have wide implications on the fracture behavior of viscoelastic materials in general.

1984

Y. Braiman, T. Egami / Physica A 388 (2009) 1978–1984

Fig. 5. Numerical (red dots) and analytical (Eq. (8), blue line) curves for T as a function of β in th logarithmic and linear (the Inset) scales. The parameters are the same as in the Figs. 1 and 2. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Acknowledgement This work was supported by the Division of Materials Science and Engineering, Office of Basic Energy Sciences(LLH), US Department of Energy under contract DE-AC05-00OR-22725 with UT-Battele, LLC. References [1] G. Wang, Y.T. Wang, Y.H. Liu, M.X. Pan, D.Q. Zhao, W.H. Wang, Appl. Phys. Lett. 89 (2006) 121909; G. Wang, D.Q. Zhao, H.Y. Bai, M.X. Pan, A.L. Xia, B.S. Han, X.K. Xi, Y. Wu, W.H. Wang, Phys. Rev. Lett. 98 (2007) 235501. [2] R. Blumenfeld, Phys. Rev. Lett. 76 (1996) 3703; and in Theoret. Appl. Fract. Mech. 30 (1998) 209. [3] A. Livne, O. Ben-David, J. Fineberg, Phys. Rev. Lett. 98 (2007) 124301. [4] J.S. Langer, Phys. Rev. A 46 (1992) 3123. [5] A.E. Lobkovsky, J.S. Langer, Phys. Rev. E 58 (1998) 1568. [6] I.S. Aranson, V.A. Kalatsky, V.M. Vinokur, Phys. Rev. Lett. 85 (2000) 118. [7] A. Karma, A.E. Lobkovsky, Phys. Rev. Lett. 92 (2004) 245510. [8] J. Kierfeld, V.M. Vinokur, Phys. Rev. Lett. 96 (2006) 175502. [9] R. Spatschek, M. Hartmann, E. Brener, H. Muller-Krumbhaar, K. Kassner, Phys. Rev. Lett. 96 (2006) 0155502. [10] Y. Braiman, T. Egami, Phys. Rev. E 77 (2008) 065101(R). [11] X.F. Zhu, Y.P. Li, G.P. Zhang, J. Tan, Y. Liu, Appl. Phys. Lett. 92 (2008) 161905. [12] J. Zhou, K. Komvopoulos, A.M. Minor, Appl. Phys. Lett. 88 (2006) 181908. [13] Y. Zhang, A.L. Greer, Appl. Phys. Lett. 89 (2006) 071907. [14] X.-P. Xu, A. Needleman, Modelling Simul. Mater. Sci. Eng. 1 (1993) 111. [15] A.M. Balk, A.V. Cherkaev, L. Slepyan, J. Mech. Phys. Solids 49 (2001) 131. [16] E. Van der Diessen, A. Needleman, Interface Sci. 11 (2003) 291. [17] V.S. Deshpande, A. Needleman, E. Van der Giessen, Acta Mater. 50 (2002) 831. [18] M.J. Buehler, F. Abraham, H. Gao, Nature 426 (2003) 141. [19] B.N.J. Persson, Phys. Rev. B 55 (1997) 8004. [20] Y.F. Lim, K. Chen, Phys. Rev. E 58 (1998) 5637. [21] E. Sharon, J. Fineberg, Nature 397 (1999) 333. [22] O. Pla, F. Guinea, E. Louis, S.V. Ghaisas, L.M. Sander, Phys. Rev. B 57 (1998) R13981.