XFEM–dislocation dynamics multi-scale modeling of plasticity and fracture

Computational Materials Science 104 (2015) 98–107

Contents lists available at ScienceDirect

Computational Materials Science journal homepage: www.elsevier.com/locate/commatsci

XFEM–dislocation dynamics multi-scale modeling of plasticity and fracture Amirreza Keyhani a, Mohsen Goudarzi a, Soheil Mohammadi a,⇑, Reza Roumina b a b

High Performance Computing Lab (HPC), School of Civil Engineering, College of Engineering, University of Tehran, Tehran, Iran School of Metallurgical and Materials Engineering, College of Engineering, University of Tehran, Tehran, Iran

a r t i c l e

i n f o

Article history: Received 2 November 2014 Received in revised form 21 February 2015 Accepted 24 March 2015

Keywords: XFEM–DD Extended multi-scale Plasticity Cohesive crack Precipitate

a b s t r a c t An extended three-dimensional multi-scale framework is presented to model plasticity in precipitate containing materials in the presence of cracks. The present framework is constructed in two scales. At the micro scale, plasticity is computed via the line dislocation dynamics (DD) methodology in which the penetrable and impenetrable precipitates are modeled. At the macro scale, application of the extended finite element method (XFEM) allows for accurate analysis of mixed-mode cohesive crack propagation without the need for any remeshing procedure. This is a significant improvement especially when two different scales are involved. While the former simulations have so far been limited to the study of only mode-I crack propagation to avoid the expensive remeshing procedure, in the present work, the effects of loading rates and precipitate density on general crack propagations are studied by utilizing the XFEM–DD multi-scale framework to illustrate the efficiency and capability of the proposed approach. Ó 2015 Elsevier B.V. All rights reserved.

1. Introduction The fracture behavior of materials is mainly controlled by two mechanisms, separation of atomic bonds in brittle materials and plastic deformation in ductile specimens. Dislocations motion is directly responsible for plastic dissipation near the crack tip region in a ductile material. The dissipative mechanisms of dislocations are only rate controlling if there is a possibility of spreading plasticity by dislocations motion around the crack tip before the crack propagates. Interfacial properties of fracture surfaces can also have a controlling effect on cracking behavior [1]. In addition, the microstructural anisotropy near the crack tip region affects the crack propagation. For instance, random distribution of precipitates or arrangement of dislocations with respect to the crack direction may create a mixed-mode crack propagation state even if the crack is loaded purely in mode I. Therefore, a multi-scale analysis, which is governed by the micro scale dislocations motion within the finite element modeling of cracking, is of great importance as it can provide more insight into the role of plasticity and microstructural features on fracture properties of crystalline materials.

⇑ Corresponding author. Tel.: +98 21 6111 2258; fax: +98 21 6640 3808. E-mail addresses: [email protected] (A. Keyhani), [email protected] (M. Goudarzi), [email protected] (S. Mohammadi), [email protected] (R. Roumina). http://dx.doi.org/10.1016/j.commatsci.2015.03.032 0927-0256/Ó 2015 Elsevier B.V. All rights reserved.

Dislocation dynamics (DD), basically designed for modeling evolutions of dislocations [2–4], has been widely used to model plasticity-related phenomena at the micro scale [5–7]. Based on the fact that in dislocation dynamics computations, the elastic fields of dislocations are usually considered in an infinite homogenous domain, DD has been coupled with the finite element method (FEM) or the boundary element method (BEM) in order to treat boundary conditions properly [8–12]. The combined framework has been used to analyze plasticity in multi-scale frameworks [11,13,14] besides its development to study the crack tip plasticity and its effects on crack propagation [15–24]. For a review of advances in discrete dislocation modeling, see [25,26]. Several authors have adopted a dislocation-independent cohesive crack model to analyze crack plasticity problems. Cleveringa et al. [16,17] and Deshpande et al. [18–20] utilized a cohesive crack model based on the universal binding law of Rose et al. [27] for analyzing the stationary crack tip plasticity, mode-I crack propagation and fatigue crack growth by dislocation dynamics in a single crystal. More recently, Olarnrithinun et al. [21] and Shishvan and van der Giessen [23,24] utilized a similar cohesive law to predict the plastic response near the stationary crack tip in anisotropic crystals and to simulate mode-I crack growth. Former studies on multi-scale modeling of crack tip plasticity and crack propagation were mainly two dimensional and limited to only mode-I fracture due to some computational restrictions. For instance, a number of studies used symmetry to model half

99

A. Keyhani et al. / Computational Materials Science 104 (2015) 98–107

of a cohesive crack and adopted axial springs to represent cohesive stresses [15–24]. Such modeling procedures did not allow for simulation of mixed-mode crack propagation because of limitations of the conventional finite element method. Therefore, mixed-mode crack propagations were practically avoided by symmetry assumptions. Modeling mixed-mode crack propagation by the finite element method requires complicated remeshing algorithms for arbitrary orientation of crack segments [28,29], which becomes extremely difficult and expensive for coupled two-scale problems. Among the available remedies to avoid these restrictions, meshfree methods, the boundary element method and even the isogeometric analysis may all seem sophisticated methodologies that can improve the simulations to some degrees [30–33]. However, the recent trends toward the extended finite element method (XFEM) have well proved its efficiency and robustness for modeling of arbitrary discontinuities [34]. Since the creation of XFEM [35–38], it has successfully been applied to cohesive zone modeling. Moës and Belytschko [39] and Wells and Sluys [40] were the first to adopt XFEM to address the mixed-mode cohesive crack propagation. Zi and Belytschko [41] proposed special tip elements without branch enrichments for cohesive cracks. The effect of higher order tip enrichment functions based on the existing analytical solutions around the tip of the cohesive cracks was studied by Cox [42], while the use of energy based criteria for accurate computation of crack propagation angle was examined numerically by Dumstorff and Meschke [43]. In addition to handling crack propagation, XFEM was adopted to model dislocations in two and three dimensions as well as in thin shell carbon nanotubes [44–46]. Coupling XFEM with the atomistic scale allowed for efficient simulation of dislocations and their motion [47,48], and XFEM–DD was developed for direct dynamic dislocation simulations in electromechanical fields [49]. For a comprehensive review of development of XFEM, see [50,51]. In the present research, an extended multi-scale framework is presented to model mixed-mode crack propagation in precipitate contained crystalline materials. At the micro scale, the line dislocation dynamics code DDLab [52,53] is modified to model dislocation–precipitate interactions [54]. The macro scale cracked domain is modeled by utilizing XFEM, which allows mixed-mode crack propagation without the need for any remeshing. The two adopted methods, XFEM and the modified dislocation dynamics are consistent in the way that both are designed to avoid complex mesh and remeshing procedures. As a result, in the whole process of multi-scale analysis, the mesh remains independent of the geometry of cracks and precipitates. This simplifies solution of large systems with arbitrary distribution of precipitates on a very simple finite element mesh, with less DOFs and low computational costs. The present XFEM–DD framework is an extension of FEM– DD framework [14], in which dislocation movements are handled by the ordinary three dimensional dislocation dynamics approach and the mixed-mode crack propagation is modeled through XFEM. 2. Governing equations 2.1. Modified line dislocation dynamics formulation In the line dislocation dynamics analysis, a dislocation curve is discretized into straight segments; each defined by its two end nodes. As a result, the overall dislocation movements are determined by the movement of segment nodes [52]. The movement of each dislocation segment contributes to the total generated macroscopic plastic strain ep , [52]

ep ¼

N X li v i ðni  bi þ bi  ni Þ 2V RVE i¼1

ð1Þ

with summation over the number of segments N. li , v i , ni and bi are the segment length, magnitude of glide velocity, unit normal vector to the slip plane, and the Burgers vector of dislocation, respectively. V RVE is the volume of the finite element in which the dislocation glide occurs. Due to significant importance of dislocation–precipitate interaction in computational materials modeling [55], several methods have been developed to model precipitates in dislocation dynamics [56–59]. An efficient computational technique to model both penetrable and impenetrable precipitates has recently been presented by Keyhani et al. [54]. In order to model impenetrable precipitates by utilizing this technique, it is considered that a dislocation node which is closer than a specific distance to a precipitate is locked. Consequently, the problem of dislocation–precipitate interaction turns into the Frank–Read (F–R) mechanism, Fig. 1. The critical stress of two mechanisms must be equal to obtain the same result. Therefore, it is considered that the dislocation rounds the precipitate with a different diameter than its original size. This modeling diameter ðDm Þ is defined by

Dm ¼ L þ D  2pLb½lnðD1 =r 0 Þ

1

ð2Þ

where L is the precipitate spacing, D is the precipitate diameter, r 0 is the core radius of a dislocation and equals to the magnitude 1

of the Burgers vector [52], D1 ¼ ðD1 þ L1 Þ , and b is a constant representing the anisotropy of dislocation line tension and considered to be 0.81 [55]. A dislocation node is locked when it reaches a distance of modeling diameter to the precipitate. At each step, the local shear stress related to the local curvature at the locked node is compared with the precipitate resistance. If the local shear stress exceeds the precipitate resistance, it is released. 2.2. XFEM–DD formulation In order to model crack propagation in the presence of dislocations, the problem can be separated into two parts by applying the superposition principle, as depicted in Fig. 2 [8–11]. The first part is handled by the dislocation dynamics theory, where dislocations and their interactions with precipitates in an infinite domain are computed (Fig. 2b). In the second part, the boundary conditions are satisfied by introducing dislocation image stresses within an XFEM framework (Fig. 2c). Similar to [16–24], the universal binding law of Rose et al. [27] is adopted to describe the constitutive response of cohesive crack. A finite domain X, surrounded by the boundary C, is defined  with the external tractions t and the prescribed displacement u on part of C. Boundaries of a cohesive crack are decomposed into two parts: the cohesive region at the crack tip, denoted by Cc , and the crack surfaces which are included in C. The traction t1 is induced from the dislocations elastic field on the boundary C of the infinite domain (Fig. 2b). The stress field resulting from the image traction t1 on the boundary C of the XFEM framework along with the dislocations elastic field in the infinite domain allow for the complete account of dislocation–surface and dislocation– crack interactions. The governing boundary value equation for an isotropic elastic body under isothermal conditions and assuming the small deformation regime can be expressed as

rrþb¼0

ð3Þ

where r is the Cauchy stress tensor and b represents the body force, including the forces due to the presence of dislocations (f B ) and the movement of dislocations ðf P Þ. Essential and natural boundary conditions of the XFEM framework are

 on a part of C u¼u

ð4Þ

100

A. Keyhani et al. / Computational Materials Science 104 (2015) 98–107

Fig. 1. A dislocation line interacting with an array of precipitates. Nodes positioned closer than a specific distance to the precipitate are locked. A dislocation segment between two precipitates with the spacing L acts as an F–R source with length Lf [54].

Fig. 2. Coupled XFEM–DD boundary value problem using the superposition principle, (a) the main problem, (b) the dislocation problem in the infinite domain, and (c) the continuum scale modeled with XFEM.

r  nC ¼ t  t1 on C

ð5Þ

r  nCc ¼ tc nCc on Cc

ð6Þ

nC and nCc are the unit normal vectors to C and Cc , respectively. The traction tc defines the transferring cohesive traction at the cohesive region at the tip of a crack ðCc Þ, which is extracted from the governing interfacial law in terms of the crack opening displacement. The final weak form of the governing equations can be written as

Z

T

dðLu uÞ rdX 

X

Z

duT bdX 

Z

duT ðt  t1 ÞdC þ

C

X

Z

T

dsut tc dCc ¼ 0

Cc

ð7Þ where dðuÞ is the displacement variation, dsut ¼ duþ  du and +/ signs correspond to upper and lower edges of the crack, respectively. Lu is defined as

2

@ @x

0

6 Lu ¼ 4 0

3 ð8Þ

@ @x

ð9Þ

I¼1

where Nstd are the basic finite element shape functions, ustd repreI I are the enrichment sent the standard nodal point unknowns, uenr I DOFs, and HðxÞ is the Heaviside enrichment function,



1 if ðx  xcr Þ  n P 0 0

ð11Þ

where

f coh ¼

R R

Cc tc Nu dCc

 u dC f ext ¼ C tN R f B ¼ X SD Bu dX R f 1 ¼  C t1 Bu dC R f p ¼ X Dep Bu dX

u ¼ ½ustd

ð12Þ

uenr  Nstd H

Nu ¼ ½N

n n X X std enr uh ðxÞ ¼ Nstd Nstd I ðxÞuI þ I ðxÞHðxÞuI

HðxÞ ¼

rBu dX þ f coh ¼ f ext þ f B þ f 1 þ f P

X

std

The XFEM approximation is now adopted to accurately account for modeling of arbitrary discontinuity paths,

I¼1

Z

with

@ 7 @y 5

@ @y

field. Assuming the infinitesimal strain theory similar to [16], the final discretized form of the XFEM–DD formulation, which accounts for the full crack–dislocation coupling based on the superposition strategy of Fig. 2 and the corresponding boundary conditions, is written based on [14],

otherwise

ð10Þ

where xcr is the nearest point to a point x on the crack surface and n is the normal vector to crack surface at the point xcr . The second part of Eq. (9) represents the discontinuous part of the displacements

ð13Þ

Bu ¼ rNu D is the elastic constitutive tensor, which is assumed to be isotropic, as widely adopted in dislocation dynamics analysis of plasticity in single crystals [16–20,22]. For instance, despite the fact that the elastic response of Cu crystal is anisotropic at the micro scale, several studies have reported low sensitivity of the results with respect to the elastic anisotropy [23,24,60]. In Eq. (13), u represents the nodal displacement vector. Superscripts std and enr stand for the standard and enrichment parts of the solution, respectively, and SD is the homogenized stress due to the presence of dislocations in the finite volume of V RVE ,

SD ¼

1 V RVE

Z XRVE

rD dX

ð14Þ

101

A. Keyhani et al. / Computational Materials Science 104 (2015) 98–107

1.5 1

1 0.5

[0 1 0] μm

[0 0 1] μm

0.5 0 −0.5

0 −0.5

−1 −1 0

−1

1

−1.5 −1

[1 0 0] μm

−0.5 0

1

0.5

[0 1 0] μm

1.5 −1.5 −1.5

−1

−0.5

0

0.5

1

1.5

[1 0 0] μm Fig. 3. The initial geometry of the problem: The red lines and the blue surface present the dislocation lines (F–R sources) and the penny crack, respectively. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 4. (a) The dislocation glide plane, (b) two-dimensional XFEM–DD model, red squares and blue circles represent the heaviside enriched nodes and the domain where f B is applied. The image force due to the crack free surfaces, f 1 , is applied on the upper crack edge, (c) Typical sub triangulation integration technique, where crosses present the integration points of triangular sub elements for evaluation of rD . (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

where rD is the stress field induced from the presence of dislocations. As a dislocation line is discretized into straight segments, rD is evaluated by summation of stress fields of all segments. To

−5

0

x 10

−0.5

K (N.μm3/2)

−1

−1.5

−2

K Eq. (22.1) I

−2.5

K XFEM−DD I

KII Eq. (22.2)

−3

K XFEM−DD II

−3.5

0

500

1000

1500

2000

2500

3000

3500

r/b Fig. 5. Comparision of numerical and analytical values of shielding stress intensity factor for h ¼ p=4.

evaluate the stress field generated by a dislocation segment, the numerical approach used in the DDLab code is adopted. For further details see [61]. The well-developed sub-triangulation technique, which has been widely used in XFEM [50,51], is adopted to evaluate the integral accurately (see typical Fig. 4c). The integrand rD in Eq. (14) is evaluated on the integration points of the triangular or tetrahedral sub-elements (for two or three dimensional problems, respectively) of the finite elements which contain a dislocation. The sub-triangulation technique ensures that the integration points are not located on the position of the singular points. The last three terms of the right hand side of Eq. (11) are related to the line dislocation dynamics. They are the equivalent force vectors due to boundary conditions and movements of dislocations. f B is the body force vector due to dislocations long range interaction in an infinite homogenous domain. Therefore, f 1 is included to represent the effect of finite boundary conditions. The last term, f p , is the body force vector from the plastic strains, which is caused by the dislocation motions. Eq. (11) is then linearized using the Newton–Raphson method. The final residual equations for iteration i can then be written as i

i

i

i

i

Ri ¼ Kui þ f coh  f ext  f 1  f B  f P

ð15Þ

R

where K ¼ X BTu DBu dX. The following incremental algorithm is iterated in each step of load increment until the solution converges,

Riþ1 ¼ Ri þ J½du

iþ1

¼0

ð16Þ

102

A. Keyhani et al. / Computational Materials Science 104 (2015) 98–107 −5

0

x 10

J¼ Kþ

K (N.μm3/2 ) I

−1 −1.5

2 0 sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi13  2 1 K KI I h ¼ 2 arctan 4 @  þ 8A5 4 K II K II

−2 −2.5

Eq. (22.1) Mesh 101x101 Mesh 125x125 Mesh 151x151

−3

0

500

1000

1500

2000

2500

3000

3. Results and discussions 3500

Fig. 6. Mesh sensitivity analysis of the results.

iþ1

 ¼ J1 Ri

½uiþ1  ¼ ½ui  þ ½du

ð17Þ iþ1



ð18Þ

where J is the Jacobian matrix (derivative of the residual vector),

J¼Kþ

@f coh @f 1 @f B @f P    @u @u @u @u

ð19Þ

Variation of forces arising from dislocations depends on dislocation movements. Nevertheless, due to the fact that dislocation movements in each step are limited to acutely small distances (less than one percent of element dimension), therefore, without the loss of accuracy, the equivalent forces due to dislocations are assumed to be constant in each step. As a result, Jacobian is further simplified to

@f 1 @u

B P ¼ @f ¼ @f ¼ 0 and the @u @u

Various simulations are performed to study the crack tip plasticity and its effects on mixed-mode crack propagation by utilizing the present framework. In Section 3.1, crack shielding due to the presence of dislocations is compared with the analytical solutions. The effects of loading ratio and precipitates density on crack propagation are studied in Sections 3.2 and 3.3, respectively. Section 3.4 is dedicated to a special case where the mixed-mode crack propagation occurs while the crack is loaded in mode I. The elastic behavior of Cu crystal is assumed to be isotropic with the shear modulus of G ¼ 54:6 GPa, the Poisson’s ratio t ¼ 0:324, and the magnitude of the Burgers vector is b ¼ 0:255 nm. Except for Section 3.1, a linear decaying cohesive law with the tensile strength f t ¼ 0:02G and the fracture energy Gf ¼ 5 N=mm is adopted. A special case of penny shape crack with isotropic growth is considered in the present simulations (i.e. the crack remains penny shape after propagation). The initial diameter of crack is 1 lm and the crack is loaded perpendicular to the crack plane (mode-I loading). Simulations are continued until the crack diameter reaches to about 2 lm. The initial geometry of the problem, the penny shape crack and the surrounding F–R sources are presented in Fig. 3. The red lines present dislocations and the blue surface

δε/δt=4.407e7

δε/δt=6.414e7

1

1

0

0

−1 −1

−1 −1

0

1

−1

0

0

1

1

δε/δt=7.242e7 1

0

0

−1 −1

−1 −1 1

−1

0

−1

0

1

δε/δt=1.051e8

1

0

ð21Þ

where K I and K II are computed by the interaction integral [39].

r/b

½du

ð20Þ

Crack propagation is activated whenever the value of normal stress at the tip of discontinuity exceeds the tensile strength of materials. Then, a new segment with a predefined length is added to the existing crack length. The direction of crack propagation is governed by the maximum hoop stress criterion in terms of the mode I and II stress intensity factors,

−0.5

−3.5

@f coh @u

1

0

1

−1

0

Fig. 7. Three-dimensional presentation of F–R sources for different applied strain rates.

1

103

A. Keyhani et al. / Computational Materials Science 104 (2015) 98–107

δε/δt=4.407e7

1.5 1

1

0.5

0.5

0

0

−0.5

−0.5

−1

−1

−1.5 −1.5

−1

−0.5

0

0.5

1

1.5

−1.5 −1.5

−1

δε/δt=7.242e7

1.5

1

0.5

0.5

0

0

−0.5

−0.5

−1

−1 −1

−0.5

0

0.5

1

1.5

−1.5 −1.5

−0.5

0

0.5

1

1.5

1

1.5

δε/δt=1.051e8

1.5

1

−1.5 −1.5

δε/δt=6.414e7

1.5

−1

−0.5

0

0.5

Fig. 8. Geometry of F–R sources in the plane of crack for different applied strain rates.

δε/δt=4.407e7

δε/δt=7.242e7

x 10

−3

δε/δt=6.414e7

−3

x 10

5

5

4

4

3

3

2

2

1

1

x 10

−3

δε/δt=1.051e8

−3

x 10

5

5

4

4

3

3

2

2

1

1

Fig. 9. Contours of the effective plastic strain for different applied strain rates.

shows the initial crack. The initial dislocation density is 0:33 lm2 for Sections 3.2–3.4. In Sections 3.2 and 3.3, the F–R sources are placed symmetric with respect to the crack plane, whereas in

Section 3.4, the lower F–R sources are shifted about 0:1 lm outward of the crack direction in order to generate a mixed-mode crack propagation mode due to microstructural non-symmetry.

104

A. Keyhani et al. / Computational Materials Science 104 (2015) 98–107

1

K dis II ¼

Crack lentgh (μm)

0.95

0.9

0.85

δε/δt=4.407e7 δε/δt=6.414e7 δε/δt=7.242e7 δε/δt=1.051e8 No dislocation

0.8

0.75

0.7 5.2

5.4

5.6

5.8

6

Applied strain εzz

x 10

−3

Fig. 10. The crack length versus the applied strain for different applied strain rates.

3.1. Crack shielding due to the presence of dislocations

app dis dis K I;II ¼ K app I;II þ K I;II where K I;II and K I;II are the applied stress intensity factor and the stress intensity factor caused by the presence of dislocations, respectively. It is called dislocation shielding if the stress intensity factor induced by presence of dislocations declines

the local stress intensity factor (i.e. K dis < 0). Conversely, it is denoted dislocation anti-shielding if it intensifies the local stress intensity factor (i.e. K dis > 0). Analytical solutions for the induced stress intensity factor due to the presence of dislocations are [62,63]

3Gb h pffiffiffiffiffiffiffiffiffi sin h cos 2 2ð1  tÞ 2pr

ð22:2Þ

where G, t and b are the shear modulus, the Poisson’s ratio and the magnitude of the Burgers vector, respectively and the subscripts refer to the crack mode. r and h are the crack tip polar coordinates, as defined in Fig. 4a. The two-dimensional XFEM model of crack–dislocation interaction is depicted in Fig 4b. It is governed by the full crack–dislocation coupling of Eq. (11) based on the superposition strategy of Fig. 2 and its corresponding boundary conditions. The stress field due to the presence of dislocations in an infinite domain is modeled by f B . While the modeling domain is considered large enough to simulate an infinite domain condition, f 1 is applied to treat the boundary conditions due to crack free surfaces. Fig. 5 compares the predicted dislocation-induced stress intensity factors with the reference analytical solution for h ¼ p=4 and 400b 6 r 6 4000b ’ 1 lm (Fig. 4a), which shows a good agreement. In addition, the mesh sensitivity of the results is evaluated in Fig. 6 for three different element sizes, which clearly illustrate that the results are mesh independent. 3.2. The effect of loading rate on mode-I cohesive crack propagation

The elastic field arising from the presence of dislocations near the crack tip may intensify or decline the stress intensity factor. The local stress intensity factor is calculated by superimposing

K dis I ¼

  Gb h h pffiffiffiffiffiffiffiffiffi 2 cos h cos  sin h sin 2 2 2ð1  tÞ 2pr

ð22:1Þ

The isotropic mode-I crack propagation of the penny crack in an infinite domain for different loading rates is studied by the present multi-scale framework. As previously mentioned, the dissipative mechanism can only operate if the fracture is delayed sufficiently. Therefore, by increasing the applied strain rate, the behavior changes from ductile to brittle. In this simulation, different applied strain rates in the [0 0 1] direction (perpendicular to the crack plane) are applied. Geometry of the F–R sources for different applied strain rates are shown in Fig. 7, which indicates that the mobility of dislocations is decreased as the applied strain rate is increased. In order to present the extent of crack propagation more clearly, the results are shown in the [1 0 0][0 1 0] view in Fig. 8. The corresponding effective plastic strains for different applied strain rates are depicted in Fig. 9, which shows that the plastic region due to dislocation motions becomes smaller by increasing the

Fig. 11. Three-dimensional presentation of F–R sources interacting with the random distribution of precipitates for different precipitate densities.

A. Keyhani et al. / Computational Materials Science 104 (2015) 98–107

105

Fig. 12. Geometry of F–R sources interacting with the random distribution of precipitates for different precipitate densities in [1 0 0][0 1 0] view.

1.1 1.05

Crack lentgh (μm)

1 0.95 0.9

No dislocation No precipitate 1% 2% 3% 4%

0.85 0.8 0.75 0.7 5.2

5.4

5.6

5.8

Applied strain εzz

6

6.2 −3

x 10

Fig. 13. Variation of the crack length versus the applied strain for different precipitate densities.

applied strain rate. Therefore, the contribution of dissipation mechanisms in the total fracture energy decreases. As a result, it is expected that the crack propagates at smaller applied strain values when the applied strain rate is larger, as clearly shown in Fig. 10.

resistance and diameter. The main purpose of this simulation is to show anisotropy of the plastic domain near the crack tip in precipitate hardened materials. This anisotropy causes mixed-mode crack propagation which is the subject of the next simulation. By increasing the density of precipitates, the dislocation motion becomes more restricted. The geometry of dislocations and precipitates for different dislocation densities are presented in Figs. 11 and 12. To improve the presentation of the geometry of dislocations, precipitates are not shown on the left half of the simulation domain. Two limiting cases are possible for each applied strain rate: one is when no dislocation exists and therefore the material behavior is completely brittle, which means that the plastic dissipation does not contribute to the total fracture energy. The other limiting case is when no precipitates exist to resist dislocation motions. Expectedly, any other cases with a random or arbitrary distribution of precipitates may be assumed in between these cases. By increasing the density of precipitates, which results in higher resistance against dislocation motions, the contribution of plastic dissipation in the total fracture energy decreases significantly, and the total fracture energy is consequently reduced. Variation of the crack length versus the applied strain for different precipitate densities are presented in Fig. 13, which shows that the crack propagates at lower applied strains by increasing precipitates density.

3.3. The effect of precipitate density on mode-I cohesive crack propagation

3.4. Mixed-mode cohesive crack propagation due to microstructural non-symmetry

The effect of random distribution of precipitates on the mode-I crack propagation is studied for different precipitate densities. All precipitates are considered to be impenetrable with the diameter of D ¼ 100 nm; however, there is no limitation on the precipitate

Microstructural non-symmetry and anisotropy affect the crack propagation direction. Such non-symmetry and anisotropy may arise from different phenomena such as the random distribution and various properties of F–R sources and precipitates. This final

106

A. Keyhani et al. / Computational Materials Science 104 (2015) 98–107

Fig. 14. Three-dimensional propagated crack topologies for two considered cases: (a) mixed-mode crack propagation due to dislocation non-symmetry with respect to the crack plane, (b) mixed-mode crack propagation due to plastic anisotropy from random distribution of precipitates.

σ (MPa)

σ (MPa)

zz

zz

1500

1500

1000

1000

500

500

0

(a)

0

(b)

Fig. 15. Stress plot ðrzz Þ for two cases: (a) dislocation non-symmetry with respect to the crack plane, (b) plastic anisotropy arising from random distribution of precipitates.

propagated crack topologies for both cases are plotted in Fig. 14. It should be noted that only the normal axis to the crack plane (and not dislocation paths) is scaled by a factor of 200. The stress plots ðrzz Þ, shown in Fig. 15, illustrate that the stress is not distributed symmetrically due to non-symmetric nature of F–R sources with respect to the crack plane (Fig. 15a). This clearly affects the crack propagation direction, as shown in Fig. 16. Comparing the results in Figs. 15 and 16, shows that in case (a)

−4

15

x 10

Initial crack Case (a) Case (b)

Height (μm)

10

5

1.1

0

1.05 1

0

0.2

0.4

0.6

0.8

1

1.2

1.4

1.6

Crack length (μm) Fig. 16. The initial and propagated crack paths; case (a) for dislocation nonsymmetry with respect to the crack plane and case (b) for random distribution of precipitates in front of the crack.

Crack lentgh (μm)

−5

Mode−I: No dislocation Mode−I: Precipitate density = 2% Mixed−mode: Case (b) Mixed−mode: Case (a) Mode−I: Precipitate density = 0%

0.95 0.9 0.85 0.8

simulation is devoted to show the capability of the present framework to handle the mixed-mode crack propagation. Here, the mixed-mode crack propagation due to microstructural non-symmetry and random distribution of precipitate are studied while the crack is purely loaded in mode I. In case (a), the lower F–R sources are shifted about 0:1 lm in the outward crack direction to disturb the symmetry. In case (b), a random distribution of precipitates with the volume density of 2% is considered. The

0.75 0.7 5.2

5.4

5.6

5.8

6

Applied strain ε

zz

6.2 −3

x 10

Fig. 17. Comparison of the crack length versus the applied strain for mode-I and mixed-mode crack propagations.

A. Keyhani et al. / Computational Materials Science 104 (2015) 98–107

when the crack tip meets the first dislocation, it starts to propagate upward due to high stress ratio resulted from the dislocation glide. After a specific propagation in the direction of first dislocation glide, it is affected by the second dislocation and moves downward. This procedure repeats until the crack tip passes all dislocations. Afterwards, the crack propagates in an almost pure mode I, as shown in Fig. 16. In Fig. 15b, the stress field is not symmetric due to plastic anisotropy arising from the random distribution of precipitates. The first lower F–R source faces high density of precipitates on its glide plane, which consequently decreases the plastic strain. As a result, the crack propagates along the first higher F–R until it reaches the second lower and higher F–R sources. After passing all F–R sources, the crack propagates in mode I. Expectedly, the glide distance of dislocations in Fig. 15b is less than Fig. 15a due to precipitate resistance against the dislocation glide. In Fig. 17, the results of crack length versus the applied strain for two cases of mixed-mode crack propagations are compared with the previous results, in which the crack was forced to propagate in mode I. Even in mixed-mode crack propagation, the presence of precipitates causes the crack to propagate at lower applied strain rates. However, after passing all F–R sources, all cases converge to the same value (Fig. 15). 4. Conclusions An extended three-dimensional multi-scale framework is presented to model plasticity and fracture in the vicinity of precipitates by merging the line dislocation dynamics analysis and the extended finite element method. At the micro-scale, the modified line dislocation dynamics is applied to analyze dislocation motions and dislocation–precipitate interaction. Since the methodology which is used to model precipitate is independent of the macro analysis, it is not required to generate a mesh to define the precipitates geometry. At the macro-scale, XFEM is applied to provide mixed-mode crack propagation without remeshing. Other available studies intentionally avoided mixed-mode crack propagation by a symmetric assumption for the crack. The present multi-scale framework is more efficient since it has to no restrictions for mixed-mode crack propagation and requires less computational cost due to decreasing DOFs and avoiding the remeshing procedure. Various simulations are performed to study the crack tip plasticity and demonstrate the capability of the present framework. It is shown that the crack propagates at lower loadings by increasing the loading rates or precipitate density since these factors decrease the share of dissipation mechanism in the total fracture energy. In addition, it is shown that the crack may propagate in mixed-mode even if it is loaded purely in mode I due to microstructural non-symmetry. Considering elastic anisotropy and improving some of the limitations of the present approach can be considered for future studies. Acknowledgements The authors would like to acknowledge the technical support of High Performance Computing Lab., School of Civil Engineering, University of Tehran. Furthermore, the financial support of Iran National Science Foundation (INSF) is gratefully acknowledged. References [1] [2] [3] [4] [5]

J.R. Rice, J.S. Wang, Mater. Sci. Eng., A 107 (1989) 23–40. A.N. Gulluoglu, C.S. Hartley, Modell. Simul. Mater. Sci. Eng. 1 (1992) 17. A.N. Gulluoglu, C.S. Hartley, Modell. Simul. Mater. Sci. Eng. 1 (1993) 383. L.P. Kubin, G. Canova, Scr. Metall. Mater. 27 (1992) 957–962. T.A. Khraishi, H.M. Zbib, Comput. Mater. Sci. 24 (2002) 310–322.

107

[6] F. Roters, D. Raabe, G. Gottstein, Comput. Mater. Sci. 7 (1996) 56–62. [7] R.N. Yellakara, Z. Wang, Comput. Mater. Sci. 75 (2013) 79–85. [8] J.A. El-Awady, S. Bulent Biner, N.M. Ghoniem, J. Mech. Phys. Solids 56 (2008) 2019–2035. [9] M.C. Fivel, C.F. Robertson, G.R. Canova, L. Boulanger, Acta Mater. 46 (1998) 6183–6194. [10] A. Needleman, Acta Mater. 48 (2000) 105–124. [11] E. Van der Giessen, A. Needleman, Modell. Simul. Mater. Sci. Eng. 3 (1995) 689. [12] H. Yasin, H.M. Zbib, M.A. Khaleel, Mater. Sci. Eng., A 309–310 (2001) 294–299. [13] M. Wallin, W.A. Curtin, M. Ristinmaa, A. Needleman, J. Mech. Phys. Solids 56 (2008) 3167–3180. [14] H.M. Zbib, T. Diaz de la Rubia, Int. J. Plast 18 (2002) 1133–1163. [15] T.K. Bhandakkar, A.C. Chng, W.A. Curtin, H. Gao, J. Mech. Phys. Solids 58 (2010) 530–541. [16] H. Cleveringa, E. Van der Giessen, A. Needleman, J. Mech. Phys. Solids 48 (2000) 1133–1157. [17] H. Cleveringa, E. Van der Giessen, A. Needleman, Mater. Sci. Eng., A 317 (2001) 37–43. [18] V. Deshpande, A. Needleman, E. Van der Giessen, Acta Mater. 49 (2001) 3189– 3203. [19] V. Deshpande, A. Needleman, E. Van der Giessen, Acta Mater. 50 (2002) 831– 846. [20] V. Deshpande, A. Needleman, E. Van der Giessen, Acta Mater. 51 (2003) 1–15. [21] S. Olarnrithinun, S.S. Chakravarthy, W.A. Curtin, J. Mech. Phys. Solids 61 (2013) 1391–1406. [22] E. Van der Giessen, V. Deshpande, H. Cleveringa, A. Needleman, J. Mech. Phys. Solids 49 (2001) 2133–2153. [23] S.S. Shishvan, E. van der Giessen, Modell. Simul. Mater. Sci. Eng. 21 (2013) 065006. [24] S.S. Shishvan, E. van der Giessen, Modell. Simul. Mater. Sci. Eng. 21 (2013) 065007. [25] M.C. Fivel, C.R. Phys. 9 (2008) 427–436. [26] H.M. Zbib, J. Eng. Mater. Technol. 131 (2009). 041209-041201. [27] J.H. Rose, J. Ferrante, J.R. Smith, Phys. Rev. Lett. 47 (1981) 675–678. [28] P. Bocca, A. Carpinteri, S. Valente, Int. J. Solids Struct. 27 (1991) 1139–1153. [29] A. Carpinteri, G. Colombo, Comput. Struct. 31 (1989) 607–636. [30] T. Belytschko, D. Organ, C. Gerlach, Comput. Methods Appl. Mech. Eng. 187 (2000) 385–399. [31] E. Ooi, Z. Yang, Eng. Anal. Bound. Elem. 33 (2009) 915–929. [32] T. Rabczuk, G. Zi, Comput. Mech. 39 (2007) 743–760. [33] C.V. Verhoosel, M.A. Scott, R. de Borst, T.J. Hughes, Int. J. Numer. Meth. Eng. 87 (2011) 336–360. [34] T. Belytschko, R. Gracie, G. Ventura, Modell. Simul. Mater. Sci. Eng. 17 (2009) 043001. [35] T. Belytschko, T. Black, Int. J. Numer. Meth. Eng. 45 (1999) 601–620. [36] J.M. Melenk, I. Babuška, Comput. Methods Appl. Mech. Eng. 139 (1996) 289– 314. [37] N. Moës, J. Dolbow, T. Belytschko, Int. J. Numer. Meth. Eng. 46 (1999) 131–150. [38] T. Strouboulis, I. Babuška, K. Copps, Comput. Methods Appl. Mech. Eng. 181 (2000) 43–69. [39] N. Moës, T. Belytschko, Eng. Fract. Mech. 69 (2002) 813–833. [40] G. Wells, L. Sluys, Int. J. Numer. Meth. Eng. 50 (2001) 2667–2682. [41] G. Zi, T. Belytschko, Int. J. Numer. Meth. Eng. 57 (2003) 2221–2240. [42] J.V. Cox, Int. J. Numer. Meth. Eng. 78 (2009) 48–83. [43] P. Dumstorff, G. Meschke, Int. J. Numer. Anal. Meth. Geomech. 31 (2007) 239– 259. [44] R. Gracie, J. Oswald, T. Belytschko, J. Mech. Phys. Solids 56 (2008) 200–214. [45] J. Oswald, R. Gracie, R. Khare, T. Belytschko, Comput. Methods Appl. Mech. Eng. 198 (2009) 1872–1886. [46] T. Belytschko, R. Gracie, Int. J. Plast. 23 (2007) 1721–1738. [47] R. Gracie, T. Belytschko, Int. J. Numer. Meth. Eng. 86 (2011) 575–597. [48] R. Gracie, T. Belytschko, Int. J. Numer. Meth. Eng. 78 (2009) 354–378. [49] O. Skiba, R. Gracie, S. Potapenko, Modell. Simul. Mater. Sci. Eng. 21 (2013) 035003. [50] S. Mohammadi, Extended Finite Element Method: For Fracture Analysis of Structures, John Wiley & Sons, 2008. [51] S. Mohammadi, XFEM Fracture Analysis of Composites, John Wiley & Sons, 2012. [52] V.V. Bulatov, W. Cai, Computer Simulations of Dislocations, Oxford University Press, 2006. [53] W. Cai, J. Deng, K. Kang, A short course on DDLab and ParaDiS (2005). [54] A. Keyhani, R. Roumina, S. Mohammadi, Submitted for publication, 2015. [55] A.J. Ardell, Metall. Trans. A 16 (1985) 2131–2165. [56] C.S. Shin, M.C. Fivel, M. Verdier, K.H. Oh, Phil. Mag. 83 (2003) 3691–3704. [57] A. Takahashi, N.M. Ghoniem, J. Mech. Phys. Solids 56 (2008) 1534–1553. [58] Y. Xiang, D.J. Srolovitz, L.T. Cheng, Acta Mater. 52 (2004) 1745–1760. [59] K. Yashiro, F. Kurose, Y. Nakashima, K. Kubo, Y. Tomita, H.M. Zbib, Int. J. Plast. 22 (2006) 713–723. [60] S.S. Shishvan, S. Mohammadi, M. Rahimian, E. Van der Giessen, Int. J. Solids Struct. 48 (2011) 374–387. [61] W. Cai, A. Arsenlis, C.R. Weinberger, V.V. Bulatov, J. Mech. Phys. Solids 54 (2006) 561–587. [62] S.M. Ohr, J. Phys. Chem. Solids 48 (1987) 1007–1014. [63] J.R. Rice, R. Thomson, Phil. Mag. 29 (1974) 73–97.