Multiscale modeling of fracture of MgO: Sensitivity of interatomic potentials

Theoretical and Applied Fracture Mechanics 53 (2010) 74–79

Contents lists available at ScienceDirect

Theoretical and Applied Fracture Mechanics journal homepage: www.elsevier.com/locate/tafmec

Multiscale modeling of fracture of MgO: Sensitivity of interatomic potentials J. Chen, J.D. Lee * Department of Mechanical and Aerospace Engineering, The George Washington University, Washington, DC 20052, USA

a r t i c l e

i n f o

Article history: Available online 16 December 2009 Keywords: Multiscale modeling Crack propagation Crack closure Sensitivity of interatomic potential

a b s t r a c t Three different interatomic potentials, namely, B-G I Model, B-G II Model and L-C Model, are used in multiscale modeling and simulation of a center-cracked specimen made of magnesia subjected to monotonically increasing loading. The specimen is decomposed into a far field, a near field and a crack-tip region. The analytical solution in the far field from linear elastic fracture mechanics (LEFM) is utilized. The solution of the near field is based on a multiscale field theory. In the crack-tip region, molecular dynamics (MD) simulation is employed. These methodologies are integrated to simulate mixed mode fracture of magnesia (MgO). Three different interatomic potentials are examined and the interatomic potential and interatomic force between Mg–Mg, Mg–O and O–O are shown. The numerical results of crack propagation demonstrate that (1) crack closure is witnessed in B-G I Model but not in B-G II Model and L-C Model, (2) B-G II Model and L-C Model diverge in the early stage. The cause of instability and the remedy are also discussed. Ó 2009 Elsevier Ltd. All rights reserved.

1. Introduction Multiscale modeling provides the means to study local physical phenomena in a large structure with microstructural features. It is worthwhile to mention several different multiscale methods: variational multiscale methods [1], homogenized Dirichlet projection method [2], domain decomposition methods [3] and Quasicontinuum (QC) method [4]. Quasicontinuum (QC) method is a mixed continuum and atomistic method to significantly reduce the numbers of atoms of which the DOF must be explicitly considered. It was originally developed by Tadmor et al. to study single crystal mechanics in 1996 [5] and later extended by Shenoy et al. to treat polycrystals and poly phase material in 1999 [6]. The complete computation domain is divided by several representative atoms. Those representative atoms will be sorted as two kinds of atoms. One is nonlocal atoms, whose energies are computed by an explicit consideration of all its neighbors and another one is local atoms, which are computed from the local deformation gradient using Cauchy–Born rule [7]. After that, the total energy will be the sum of the energy of local atoms and nonlocal atoms. However, due to the atoms in the transient zone, the total energy should be corrected by a ghost force function. To state the system in equilibrium, it demands the gradient of the potential energy is zero, which implies only static cases are concerned. However, among those multiscale approaches that use interatomic potential model, none of them studies the sensitivity of

* Corresponding author. E-mail address: [email protected] (J.D. Lee). 0167-8442/$ - see front matter Ó 2009 Elsevier Ltd. All rights reserved. doi:10.1016/j.tafmec.2009.12.007

interatomic potentials. It is necessary to discuss whether slightly different interatomic potentials lead to similar results or not. In this work, the center-cracked specimen is decomposed into three parts: (1) far field, modeled by linear elastic fracture mechanics (LEFM), (2) near field (cf. Fig. 1), modeled by a multiscale field theory [8] and analyzed by a generalized finite element method (GFEM) [9] and (3) crack-tip region, modeled by molecular dynamics. The exact and analytical solution of the far field from LEFM is utilized to specify boundary conditions at the interface between the near field and the far field. Full-blown interatomic forces are employed between the near field and the crack-tip region. 2. Theories for regions of different length scales Let the center-cracked specimen be decomposed into three parts (cf. Fig. 1). The far field is governed by the linear elastic fracture mechanics (LEFM). At the crack tip, it is modeled by MD simulation. In between, named as the near field, it is simulated by GFEM [9]. In LEFM, the exact and analytical Sneddon’s solution of stress field and displacement field for an infinite plate with a centercrack, in 2D plane stress or plane strain, can be expressed as (cf. Fig. 2) [10–12]:     K1 h þ h2 r 1 a2 3  txx ¼ pffiffiffiffiffiffiffiffiffiffiffiffi r1 cos h1  sin h1 sin ðh þ h2 Þ 2 parr2 2 rr2     K2 h þ h2 r1 a2 3 þ pffiffiffiffiffiffiffiffiffiffiffiffi 2r1 sin h1   sin h1 cos ðh þ h2 Þ þ v  r; 2 parr2 2 rr2 ð1Þ

75

J. Chen, J.D. Lee / Theoretical and Applied Fracture Mechanics 53 (2010) 74–79

   K1 h þ h2 h þ h2  2r 21 sin h1 sin h1  4lux ¼ pffiffiffiffiffiffiffiffiffiffiffiffi ðj  1Þrr2 cos parr 2 2 2    K2 h þ h2 h þ h2 þ pffiffiffiffiffiffiffiffiffiffiffiffi ðj þ 1Þrr2 sin þ 2r21 sin h1 cos h1  parr2 2 2  0:5r 1 cos h1 ðr  vÞðj þ 1Þ; ð4Þ



  K1 h þ h2 h þ h2 4luy ¼ pffiffiffiffiffiffiffiffiffiffiffiffi ðj þ 1Þrr 2 sin  2r 21 sin h1 cos h1  parr2 2 2    K2 h þ h2 h þ h2 þ pffiffiffiffiffiffiffiffiffiffiffiffi ð1  jÞrr 2 cos  2r 21 sin h1 sin h1  parr2 2 2  0:5r1 sin h1 ðr  vÞðj  3Þ; ð5Þ pffiffiffiffiffiffi pffiffiffiffiffiffi where l is the shear modulus; K 1 ¼ r pa; K 2 ¼ s pa;

(

j¼

Fig. 1. Finite element mesh for near field and atomic region. Only the tip region of the right half plane is modeled. In the center, it is the crack-tip region, which is modeled by MD simulation. Outside the center, it is the near field, which is modeled by the generalized finite element method. The displacement imposed on the boundary of near field is obtained from the Sneddon’s solution.

    K1 h þ h2 r 1 a2 3 t yy ¼ pffiffiffiffiffiffiffiffiffiffiffiffi r 1 cos h1  sin h1 sin ðh þ h2 Þ  2 parr2 2 rr 2 K 2 r 1 a2 3 þ pffiffiffiffiffiffiffiffiffiffiffiffi sin h1 cos ðh þ h2 Þ; 2 parr2 rr2

ð2Þ

dðv þ Dv a Þ ¼f dt

a

þ /a þ rx  t a½kinetic þ rya  sa½kinetic ;

ð7Þ

where the superscript a refers to the ath atom in the lattice cell; the terms on the RHD are respectively interatomic force density, body force density, and kinetic part of homogeneous and inhomogeneous stresses due to temperature. The kinetic part of stresses indicates that temperature is an independent variable in this atom-embedded continuum theory, which is different from MD simulation where temperature can be computed by velocity field. In this work, only the interatomic force is considered. The weak form of GFEM is [14] Ng Z X

K1 3 t xy ¼ pffiffiffiffiffiffiffiffiffiffiffiffi sin h1 cos ðh þ h2 Þ 2 parr2     K2 h þ h2 r 1 a2 3  þ pffiffiffiffiffiffiffiffiffiffiffiffi r 1 cos h1  sin h1 sin ðh þ h2 Þ ; 2 parr2 2 rr 2

ð6Þ

for plane stress:

It is noticed that near the crack tip there is only plane strain because the distance r, as it approaches zero, is always smaller than any plate thickness. Therefore, plane strain always prevails. This is why ASTM only recognizes plane strain fracture toughness value. That is KIC. In this work, Sneddon’s solution is only used for the far field. In other words, the displacements along the interface between the near field and the far field are used as the boundary conditions imposed on the near field (cf. Fig. 1). The near field solution can be obtained through GFEM based on the multiscale field theory [8,9], which can be reduced to MD simulation if the size of the element is reduced to the lattice constant. On one hand, it is an atom-embedded continuum theory and therefore, detailed atomic motion can be observed in the continuum region. On the other hand, it provides the accuracy close to MD simulation but in a much more computationally effective way. The balance law of linear momentum can be written as [8,13]:

qa

Fig. 2. Geometry and loading condition of infinite center-cracked specimen with a being half of the crack size [15].

3  4t for plane strain; 3t 1þt

K¼1

mðKÞ Na X X € a ðKÞ  f a gdua ðKÞdXK þ € i ðKÞ  Fi gdui ¼ 0: fqcða;KÞ u fmcðiÞ u

XK a¼1

ð8Þ

i¼1

The crack-tip region is modeled by MD simulation: 2

ð3Þ

mi

d ri dt

2

¼ Fi ¼

X

F ij ;

ð9Þ

j

Table 1 Three different interatomic potential models of MgO [15]. Mg–Mg

Mg–O

O–O

L-C

B-G I

B-G II

L-C

B-G I

B-G II

L-C

B-G I

B-G II

Aij

0

0

0

30.193

47.2

34.165

836.5

350.88

178.97

Bij

0

0

0

0.61265

0.56635

0.5632

C ij

0

0

0

0

0

0

0.28157 46.659

0.41415 54.09

0.50456 128.86

76

J. Chen, J.D. Lee / Theoretical and Applied Fracture Mechanics 53 (2010) 74–79

ij Zi Zj r þ Aij e Bij  C ij ðrij Þ6 ; rij ( ) @U ij Z i Z j Aij  rijij ij ij 8 ij B F ¼ i ¼ ½ri  rj : þ e þ 6C ðr Þ @r ðr ij Þ3 Bij r ij

U ij ¼

ð10Þ

ð11Þ

Eq. (9) is Newton’s second law, which governs all the detailed atomic motion. It is noticed that the interatomic force is derived from the interatomic potential, e.g., the Coulomb–Buckingham potential in Eq. (10). Eq. (11) gives the interatomic force from the Coulomb–Buckingham potential. Full-blown interatomic forces are employed between the near field and the crack-tip region.

Fig. 3. Interatomic potential of three models for Mg–Mg interaction.

Fig. 6. Interatomic force of three models for Mg–Mg interaction.

Fig. 4. Interatomic potential of three models for Mg–O interaction.

Fig. 7. Interatomic force of three models for Mg–O interaction.

Fig. 5. Interatomic potential of three models for O–O interaction.

Fig. 8. Interatomic force of three models for O–O interaction.

77

J. Chen, J.D. Lee / Theoretical and Applied Fracture Mechanics 53 (2010) 74–79 Table 2 Results from three different interatomic potentials of MgO. L-C Model

B-G I Model

B-G II Model

After 4000 step relaxation

10,000th step

15,000th step

20,000th step

25,000th step

32,500th step

47,500th step

(continued on next page)

78

J. Chen, J.D. Lee / Theoretical and Applied Fracture Mechanics 53 (2010) 74–79

Table 2 (continued) L-C Model

B-G I Model

B-G II Model

62,500th step

Three different models of Coulomb–Buckingham Potential are utilized in this work and listed in Table 1 in atomic units [15,16]. In L-C and B-G I Models, formal charges (2e) are assumed while in B-G II Model partial charges (1:7e) are assumed. The interatomic potential (1 Hartree ¼ 4:3597482  1018 Joule) vs. atomic distance (1 Bohr ¼ 5:291772108  1011 m) of all the three models for Mg–Mg, Mg–O, and O–O interactions are shown in Figs. 3– 5, respectively. The interatomic force (1 Hartree=Bohr ¼ 8:2387225  108 Newton) vs. interatomic distance (Bohr) of all the three models for Mg–Mg, Mg–O, and O–O interactions are shown in Figs. 6–8, respectively. In Fig. 5(Fig. 8), there is a potential(force) barrier to prevent the collision of two oxygen atoms at a critical distance. This means if the interatomic distance is less than the critical distance, the interatomic force becomes attractive from being repulsive. The attractive force will make the interatomic distance even smaller and thus induces instability. It is seen that from Eq. (10), when the interatomic distance is getting smaller, the C ij ðrij Þ6 term dominates and approaches to negative infinity. It should be emphasized that the strength (height) of the barrier in B-G I Model is almost fourfold of that in L-C Model and B-G II Model. 3. Simulation results and remedy The material used in this work is magnesia. It has eight atoms in a conventional unit cell. The entire solution domain, including MD simulation region, is 141  101  3 lattice points. The lattice constant is equal to 7.937 Bohr (1 Bohr ¼ 5:291772108  1011 m). The relaxation time is equal to 4000 time steps. Each time step is 40 atomic units, which is 40  ð2:418884326505 1017 sÞ  9:68  1016 s. The half of crack size, a, is 10,000 Bohr, i.e., 10; 000  ð5:291772108  1011 Þ  5:3  107 m. In mixed mode fracture, both shear stress and normal stress, s and r, are applied; the shear stress will dislocate the atoms and the normal stress will pull the atoms apart. This combination results in cracks. Table 2 shows the numerical results of atomic motions in the neighborhood of crack tip in three models of magnesia. In B-G I Model, crack closes up in the wake of advancing crack tip; in L-C Model and B-G II Model, the crack does not close up. After 4000 relaxation steps, the crystal will have its own optimized lattice constant for each different model. At 17349-th step, B-G II Model diverges. At 32734-th step, L-C Model diverges. At 62500-th step (last step in the simulation), B-G I Model does not diverge but the crack closes up in the wake of the advancing crack tip. It proves that the MD simulation is very sensitive to interatomic potentials; even for the same material, slightly different potential leads to completely different results. Similar to radiation damage, in fracture mechanics, atoms can get very close to one another. To prevent the disturbance of the energy level of the atomic shells, i.e., to activate the repulsion at very small distances, two approaches can be used. In a low energy problem, a r12 term can simply be added with the same coefficient as

the r 6 term [16]. In a high energy problem, Ziegler–Biersak–Litmark (ZBL) potential form is suggested and can be represented universally [17]. The ZBL potential is a screening electrostatic potential of nucleon–nucleon interaction as

U ij ¼

Zi Zj /ðr ij Þ; r ij

/ðr ij Þ ¼

4 X

Am e

ð12Þ

bm rij

;

ð13Þ

aBohr ;

ð14Þ

au

m¼1

au ¼

0:8854 þ Z 0:23 Z 0:23 i j

where aBohr is the Bohr radius, equals to 0.539177 angstroms (1 Angstrom ¼ 1010 m). 4. Conclusion In this work, the specimen is idealized as an infinite plate with a center-crack solved by different approaches for different regions of length scales. Between the classical fracture mechanics and the molecular dynamics (MD) simulation, a multiscale field theory bridges the gap. From the beginning to the end, the sight of atoms is always in sight. From Figs. 4–6, the difference of three interatomic potential models is witnessed. That is the fundamental reason why the results of atomic motions are different for the same material in the same physical problem. Furthermore, in O–O interaction, when the interatomic distance is getting smaller and over the barrier, the dominating term will approach to negative infinity meaning highly attractive force, which induces the instability of the MD simulation. From Table 2, the detailed comparison of atomic motions is shown. B-G I Model gives crack closure but other two models, BG Model II and L-C Model, do not agree with this result but instead diverge in the early but different stages. The problem of crack propagation is similar to that of radiation damage, where the atoms are forced to be close to one another. Such a behavior can be simulated by adding another term or to use ZBL potential. To reiterate, adjustment of the interatomic potential becomes necessary for simulating the real physical behavior, particularly in regions in the order of atomic distance such as the neighborhood of crack tip. References [1] T.J.R. Hughes, Multiscale phenomena: Green’s function, the Dirichlet-toNeumann formulation, subgrid scale models, bubbles and the origins of stabilized methods, Computer Methods in Applied Mechanics and Engineering 127 (1995) 387–401.

J. Chen, J.D. Lee / Theoretical and Applied Fracture Mechanics 53 (2010) 74–79 [2] T.I. Zohdi, J.T. Oden, G.J. Rodin, Hierarchical modeling of heterogeneous bodies, Computer Methods in Applied Mechanics and Engineering 138 (1996) 273– 298. [3] S.P. Xiao, T. Belytschko, A bridging domain method for coupling continua with molecular dynamics, Computer Methods in Applied Mechanics and Engineering 193 (2004) 1645–1669. [4] R. Miller, E.B. Tadmor, R. Phillips, M. Ortiz, Quasicontinuum simulation of fracture at the atomic scale, Modeling and Simulation in Material Science and Engineering 6 (1998) 607–638. [5] E.B. Tadmor, R. Phillips, M. Ortiz, Quasicontinuum analysis of defects in solids, Philosophical Magazine A 73 (1996) 1529–1563. [6] V.B. Shenoy, R. Miller, E.B. Tadmor, D. Rodney, R. Phillips, M. Ortiz, An adaptive finite element approach to atomic-scale mechanics – the quasicontinuum method, Journal of the Mechanics and Physics of Solids 47 (1999) 611–642. [7] M. Born, K. Huang, Dynamical Theory of Crystal Lattices, Oxford University Press, 1954. [8] Y. Chen, J.D. Lee, Atomistic formulation of a multiscale theory for nano/micro physics, Philosophical Magazine 85 (2005) 4095–4126.

79

[9] J.D. Lee, X. Wang, Y. Chen, Multiscale computation for nano/micro materials, Journal of Engineering Mechanics 135 (3) (2009) 192–202. [10] I.N. Sneddon, The distribution of stress in the neighborhood of a crack in an elastic solid, Proceedings, Royal Society of London A-187 (1946) 229–260. [11] I.N. Sneddon, M. Lowengrub, Crack Problems in the Classical Theory of Elasticity, John Wiley and Sons, 1969. [12] Y. Chen, J.D. Lee, A. Eskandarin, Meshless Methods in Solid Mechanics, Springer, Washington, DC, 2006. [13] Y. Chen, J.D. Lee, Conservation laws at nano/micro scales, Journal of Mechanics of Materials and Structures 1 (2006) 681–704. [14] J.D. Lee, X.Q. Wang, Y. Chen, Multiscale material modeling and its application to a dynamic crack propagation, Theoretical and Applied Mechanics 51 (2009) 33–40. [15] G. Henkelman, B.P. Uberuaga, D.J. Harris, J.H. Harding, N.L. Allan, MgO addimer diffusion on MgO(1 0 0): a comparison of ab initio and empirical models, Physical Review B 72 (2005) 115437. [16] R. Grimes, Imperial College, UK, Private Communication. [17] J.F. Ziegler, J.P. Biersack, U. Littmark, The Stopping and Range of Ions in Matter, Pergamon, New York, 1985.