Developing micro-scale crystal plasticity model based on phase field theory for modeling dislocations in heteroepitaxial structures

International Journal of Plasticity 81 (2016) 267e283

Contents lists available at ScienceDirect

International Journal of Plasticity journal homepage: www.elsevier.com/locate/ijplas

Developing micro-scale crystal plasticity model based on phase field theory for modeling dislocations in heteroepitaxial structures Liyuan Wang, Zhanli Liu*, Zhuo Zhuang** Applied Mechanics Lab., School of Aerospace Engineering, Tsinghua University, Beijing 100084, China

a r t i c l e i n f o

a b s t r a c t

Article history: Received 7 June 2015 Received in revised form 12 January 2016 Available online 3 February 2016

The capability of conventional crystal plasticity theory is extended in this paper to model a single dislocation behavior in heteroepitaxial structures. The plastic slip associated with each slip system is described by a continuous smooth field independently, which is the phase field interpolating between a slipped and an unslipped region. A dislocation is identified with location where the phase field value changes smoothly, to represent a smeared dislocation. Under a thermodynamically consistent framework that distinguishes between stored energy and dissipated energy during plastic deformation, the coupled balance equations of plastic slip evolution and quasi-static stress equilibrium are derived by using the principle of virtual power. With numerical implementation by finite element method, it is flexible to deal with material anisotropy and elastic modulus mismatch in heteroepitaxial structures at micro-scale, which are advantages of the proposed model compared to Khachaturyan-type phase field model of dislocations. Another advantage is that it is straightforward to handle the interaction and co-evolution of several types of material microstructures, such as dislocations and interfaces, within a unified continuum mechanics framework. Several examples are presented to illustrate the applicability and accuracy of the new method in modeling dislocations in complex heteroepitaxial structures. © 2016 Elsevier Ltd. All rights reserved.

Keywords: A. Dislocations C. Finite elements B. Crystal plasticity Phase field

1. Introduction Heteroepitaxial structures, e.g. epitaxial films or coreeshell nanopillars, have recently received much attention because of their important applications in engineering, such as electronics, optoelectronics and solar cells (Freund, 2000; Fu et al., 2004; Panda and Tseng, 2013). However, dislocation-free heteroepitaxial structures cannot be grown with arbitrary thickness and misfit dislocation will form above a critical thickness. The dislocation and its strong interaction with material interface not only influence the mechanical properties such as fracture toughness and strength, but also change the electrical or optical properties of heteroepitaxial crystalline materials (Bennett, 2010; Zbib et al., 2011; Liu et al., 2013; Abdolrahim et al., 2014). Understanding the underlying dynamics of dislocations in heteroepitaxial structures is crucial to provide a better insight into

* Corresponding author. Tel./fax: þ86 10 62783014. ** Corresponding author. Tel./fax: þ86 10 62783014. E-mail addresses: [email protected] (Z. Liu), [email protected] (Z. Zhuang). http://dx.doi.org/10.1016/j.ijplas.2016.01.010 0749-6419/© 2016 Elsevier Ltd. All rights reserved.

268

L. Wang et al. / International Journal of Plasticity 81 (2016) 267e283

the material design and property prediction. However, there is still lacking of an effective tool to study the dislocation behaviors in complex heteroepitaxial structures. In the past decades, discrete dislocation dynamics (DDD) that explicitly tracks dislocation lines has proven to be a useful simulation method in modeling dislocations in homogeneous bulk crystals (Amodeo and Ghoniem, 1990; Zbib et al., 1998). However, it is difficult to deal with the dislocations in heteroepitaxial structures due to the image force produced by the interfaces. Currently, most DDD methods are based on the superposition of analytical solutions of stress fields for dislocations in infinite domain. However, the existence of multi-material Green's functions is limited (Han and Ghoniem, 2005). The solution on a finite domain is the sum of infinite domain dislocation solutions and the image field corrections from the model boundary, which is usually calculated by boundary element method (El-Awady et al., 2008; Takahashi and Ghoniem, 2008) or finite element method (FEM) (Giessen and Needleman, 1995; Fivel and Canova, 1999; Weygand et al., 2001). The superposition DDD method cannot reflect the basic concept of plastic strain and is limited by relying on the existed analytical solutions. As an alternative, the coupled discrete-continuum model (DCM) was proposed, in which dislocations are associated with eigenstrain and the equilibrium field is determined directly by FEM (Lemarchand et al., 2001; Zbib and de la Rubia, 2002; Liu
et al., 2014). The critical thickness of epitaxial films for plastic relaxation et al., 2009; Gao et al., 2010; Cui et al., 2014; Vattre (Groh et al., 2003) and the dislocation behaviors in heteroepitaxial films (Cui et al., 2015a) have been studied by this method. However, the DCM still needs the analytical solution to modify the stress field when the dislocation is close to the interface (Cui et al., 2015b). In addition, how to handle the finite deformation is still a challenging issue that all the discrete based models are facing. Recently, several dislocation modeling methods completely based on continuum mechanics framework are developed. The most advantage of these methods is that the equilibrium field is solved directly. Based on field dislocation mechanics, a general crystal plasticity model was developed by Acharya (2001). The dislocation density serves as the primary internal variable and a finite element discretization of the model was presented (Roy and Acharya, 2005). It has been used to study the mechanical response of multi-layer thin films (Puri et al., 2011). Inspired by crack modeling using extended finite element method (XFEM), the displacement fields of dislocation are approximated by the sum of standard finite element part and discontinuous enriched part (Belytschko and Gracie, 2007; Gracie et al., 2007). However, the enrichment for the dislocation core is not available for anisotropic materials (Gracie et al., 2008; Oswald et al., 2009). Among all the dislocation modeling methods, the phase field method (PFM) which represents the discontinuities associated with dislocations by regularization, is receiving more and more attentions (Khachaturyan, 1983; Hu and Chen, 2001; Wang et al., 2001; Rodney et al., 2003; Koslowski et al., 2004). The advantage is that it can automatically take into account the interaction and co-evolution of dislocations and interfaces within a unified continuum mechanics framework (Wang et al., 2003; Lei et al., 2013). Currently, most phase field models of dislocations are based on the microelasticity theory of Khachaturyan and Shatalov (Khachaturyan, 1967; Khachaturyan and Shatalov, 1969; Khachaturyan, 1983). In the model, it takes ‘stress-free’ inelastic strain to represent dislocations. The elastic strain energy caused by dislocations is then expressed as a closed form function of the inelastic strain through the exact Green's function. While, the analytical Green's function solution is not available for complex structures or complex boundary conditions. For elastically inhomogeneous structures, the virtual misfit strain which is considered as additional phase fields has to be introduced and increases the number of equations to be solved (Wang et al., 2002). The dislocation behaviors in heteroepitaxial thin films have been studied by the phase field microelasticity model (Wang et al., 2003). The film and the substrate are assumed having the same elastic modulus to avoid dealing with the additional misfit strain (Wang et al., 2003). On the other hand, the equations are solved by the fast Fourier transform (FFT) method, which restricts the application to complex structures, like heteroepitaxial coreeshell nanopillars. Besides, based on Ginzburg-Landau equation, the model does not distinguish between stored energy and dissipated energy during plastic deformation (Wang et al., 2001, 2003; Lei et al., 2013). In this work, a micro-scale crystal plasticity model based on phase field theory is developed to overcome the difficulties in the Khachaturyan-type phase field model of dislocations. In contrast, the ‘stress-free’ inelastic strain is directly considered as the plastic strain based on crystal plasticity theory in the proposed model. The plastic slip associated with each slip system is described by an independent phase field variable to model a single dislocation. The elastic strain energy is expressed as a function of the elastic strain through the crystal plasticity constitutive model. Under a thermodynamically consistent framework that distinguishes between stored energy and dissipated energy during plastic deformation, the coupled balance equations of plastic slip evolution and quasi-static stress equilibrium are derived by using the principle of virtual power. Then the boundary value problem is solved directly by FEM. It can be used for complex structures or complex boundary conditions where the analytical Green's function solution is not available, which is one advantage of the proposed model. Another advantage is that the elastic modulus mismatch in heteroepitaxial structures is easy to be treated without additional complications. Besides, with numerical implementation by FEM, it is flexible to deal with finite deformation in heteroepitaxial structures at micro-scale. An outline of this paper is given as follows. The theoretical model is developed in Section 2. Three computational demonstrations are presented in Section 3, which include a screw dislocation near a free surface, a screw dislocation in an anisotropic material and an edge dislocation near a bimaterial interface. Through these examples, the accuracy of the proposed model is studied by comparing with analytical solutions. Dislocations in heteroepitaxial structures are simulated in Section 4. The conclusions are provided in Section 5.

L. Wang et al. / International Journal of Plasticity 81 (2016) 267e283

269

2. Theoretical model 2.1. Basic equations of crystal plasticity theory In this study, small strain conditions are assumed for simplicity where crystal geometry changes are neglected. The formulation extended to large deformation is more complicated and will be considered in the future work. So the total strain ε can be divided into elastic and plastic parts, εe and εp, as (Gurtin, 2002; Kuroda and Tvergaard, 2008)

ε ¼ Vs u ¼ εe þ εp

(1)

where u is displacement field and Vs is symmetric gradient operator. The plastic strain in crystal plasticity model is given by Gurtin (2002) and Kuroda and Tvergaard (2008)

εp ¼

X

Pa ga ; Pa ¼

a

1 ðsa 5ma þ ma 5sa Þ 2

(2)

where ga is a scalar variable that represents the plastic slip on slip system a, sa and ma are the unit vectors of the slip direction and the slip plane normal of slip system a, respectively. The expression of plastic slip ga associated with each slip system a is derived in Section 2.2. The crystal plasticity constitutive equation is expressed as

s ¼ C : εe ¼ C : ðε  εp Þ

(3)

where s is Cauchy stress, C is elastic modulus tensor which may be either isotropic or anisotropic. 2.2. Phase field description of plastic slip To extend the capability of conventional crystal plasticity theory to model a single dislocation, the same description for dislocations as in the phase field is adopted here. In phase field microelasticity theory, the total strain ε can be divided into elastic and inelastic parts, εe and εinelast, as (Khachaturyan, 1967; Khachaturyan and Shatalov, 1969; Khachaturyan, 1983).

ε ¼ εe þ εinelast

(4)

By comparing Eq. (4) with Eq. (1), the inelastic strain is directly considered as the plastic strain in the proposed model. By taking fa as the phase field of slip system a, the inelastic strain associated with dislocations from all slip systems is expressed as (Nabarro, 1951; Khachaturyan, 1983; Shen and Wang, 2003), i.e.

εinelast ¼

X a

dis εdis a fa ; εa ¼

ba ðsa 5ma þ ma 5sa Þ 2d

(5)

where εdis a is the eigenstrain of a dislocation loop surrounding the slipped region on slip system a, d is the interplanar distance of the slip planes and ba is the magnitude of Burgers vector. By comparing Eq. (5) with Eq. (2), the plastic slip ga associated with each slip system a can be expressed by the phase field fa as below.

ga ¼

ba f d a

(6)

Here the value of phase field fa represents the amount of plastic slip in a unit ba/d. Compared with the conventional crystal plasticity theory, the plastic slip ga associated with each slip system a is described by an independent phase field variable fa to model a single dislocation, instead of by a set of phenomenological evolution equations. As shown in Fig. 1(a), the dislocation is identified with location where the phase field value changes smoothly in the crystal plane and thus has finite dislocation core width. It represents a regularized dislocation topology which is different from the Volterra dislocation model, as shown in Fig. 1(b). The regularization is adjusted by the interplanar distance of the slip planes d and gives d / 0 for the Volterra dislocation. 2.3. Stored energy and dissipated energy By exploiting the above phase field description of plastic slip, the evolution of individual dislocations can be described. The dislocation evolution is driven by the energetic driving forces of the system. The stored energy of dislocation system consists of four parts: elastic energy Eel, crystalline energy Ecryst, gradient energy Egrad and surface energy Esurf.

270

L. Wang et al. / International Journal of Plasticity 81 (2016) 267e283

Fig. 1. (a) Phase field dislocation model; (b) Volterra dislocation model.

E ¼ Eel þ Ecryst þ Egrad þ Esurf

(7)

The elastic energy is associated with the elastic strain of crystal lattice caused by the dislocation and external loading. It is expressed as a function of the elastic strain as below.

Z Eel ¼

1 e ε : C : εe dV 2

(8)

V

It should be noted that the formulation of elastic energy here is different from the one in Khachaturyan-type phase field model, which is a closed form function of the inelastic strain through the exact Green's function (Khachaturyan, 1967; Khachaturyan and Shatalov, 1969; Khachaturyan, 1983). Different from the conventional crystal plasticity theory, two additional energy terms, including the crystalline energy Ecryst and gradient energy Egrad, are introduced here to describe the recoverable strong nonlinear deformation behavior around the dislocation core. The crystalline energy represents the potential energy of the crystal subjected to general slips and can be directly related to the generalized stacking fault energy for a particular slip plane. The slip vector associated with each slip system tends to integer multiples of the lattice vector naturally, which can be achieved by a choice of periodic function. For a full dislocation, the crystalline energy can be approximated by a sinusoidal function. The simplest expression by decoupling different slip systems was given by Wang et al. (2001),

Ecryst ¼ A

X Z a

sin2 ðpfa ÞdV

(9)

V

where A ¼ mðε0 Þ2 =2p2 is a constant (the unstable stacking fault energy), m is the shear modulus associated with the corresponding slip system and ε0 ¼ b/d is the typical shear strain. The gradient energy represents the higher-order energy associated with the higher-order derivatives of the displacement. It helps smooth and widen the dislocation core for numerical stability. It is proportional to the length of the dislocation, which can be achieved by using the term ma  Vfa. The simplest expression by decoupling different slip systems was given by Wang et al. (2005),

Egrad ¼ B

X Z a

ðma  Vfa Þ2 dV

(10)

V

where B is a positive constant. The surface effect becomes accentuated as heteroepitaxial structures become increasingly small. A form of stored energy generated by the surface steps as dislocations exit a crystal was given by Hurtado and Ortiz (2012). A diffusional surface model was also developed to simulate surface roughening in heteroepitaxial films (Wang et al., 2004). Currently, the surface energy is neglected for simplicity. The specific form of surface energy Esurf will be considered in the future work. Then, substituting Eqs. (8)e(10) into Eq. (7) and taking the derivative of Eq. (7) with respect to time, the energy storage rate can be achieved as follows.

E_ ¼

# Z " X  e _ _ s : ε_ þ Ap sinð2pfa Þfa þ 2Bma  Vfa  ma $Vfa dV V

a

(11)

L. Wang et al. / International Journal of Plasticity 81 (2016) 267e283

271

On the other side, the dislocation evolution is also an energy dissipation process. In this work, it is implemented by using a dissipation function which only depends on the rate of phase fields for simplicity,

D¼

1 X k 2 a

Z



2 f_ a dV

(12)

V

where k is dissipation coefficient. It should be pointed out that the more complex dissipation functions, e.g. also depending on the gradient of phase fields, can be used. 2.4. Principle of virtual power Regarding to the displacement field, the boundary surface of a body can be decomposed by S ¼ Su ∪St , where displacement u0 and surface traction t are prescribed on Su and St, respectively. Regarding to the phase fields, the boundary surface of a body can be decomposed by S ¼ Sfa ∪Sqa , where phase fields fa and higher-order surface force qa are prescribed on Sfa and Sqa , respectively. Clearly, there are Su ∩St ¼ 0 and Sfa ∩Sqa ¼ 0. Then the external power is expressed as follows,

Z P¼

_ f$udV þ

Z

_ þ t$udS

a

st

V

X Z

qa f_ a dS

(13)

sqa

where f is a given body force. Under a thermodynamically consistent framework, the macroscopic and microscopic balance equations are derived by using the principle of virtual power. With energy storage rate Eq. (11), dissipation function Eq. (12) and external power Eq. (13) given above, it requires

dE_ þ dD ¼ dP

(14)

for all admissible rates of the displacement field du_ and phase fields df_ a . The balance Eq. (14) makes the internal virtual power equilibrate with the external virtual power, where the former includes stored and dissipated parts. With the formulations given above, we have

# Z " X  e _ _ _ _ kfa dfa þ 2Bma  Vfa  ma $Vdfa þ Ap sinð2pfa Þdfa dV s : d_ε þ a

V

Z ¼

_ þ t$dudS

Z V

S

_ f$dudV þ

X Z a

qa df_ a dS

(15)

S

Interestingly, Eq. (15) has some similarities with the counterpart of higher-order theory of single crystal plasticity. The energy balance in the higher-order crystal plasticity can be expressed as below (Gurtin, 2002; Kuroda and Tvergaard, 2008),

# Z " Z Z X X Z  e p _ þ _ t$dudS f$dudV þ pa dg_ a dS ta dg_ a þ xa $Vdg_ a dV ¼ s : d_ε þ V

a

S

V

a

(16)

S

where g_ a is the plastic slip rate on slip system a, tpa , xa and pa are the slip resistance, higher-order force and higher-order surface force on slip system a, respectively. By comparing Eq. (15) with Eq. (16), we can find that kf_ a corresponds to the slip resistance tpa , 2Bma  Vfa  ma corresponds to the higher-order force xa and qa corresponds to the higher-order surface force pa, respectively, in a unit d/ba. While, Eq. (15) has one more item Ap sinð2pfa Þdf_ a than Eq. (16). With this item, the proposed micro-scale crystal plasticity model can be used at a smaller size scale where a single dislocation behavior is important than the higher-order theory of single crystal plasticity. To carry out a comparison study with the higher-order theory, the evolution of a dislocation loop due to its self-stress is studied. As shown in Fig. 2(a), a circular dislocation loop of radius R ¼ 30 nm in a cubic domain with the side length L ¼ 120 nm is considered. The dislocation loop is located at the center of the cubic domain and lies in the xey plane. The Burgers vector b is parallel to the x-axis with the magnitude of 0.25 nm. The material properties are E ¼ 100 GPa and n ¼ 0.34. A hexahedral mesh with element size of l ¼ 5 nm is used. The initial phase fields are set the same in our proposed micro-scale crystal plasticity model and the higher-order theory, to model the pre-existing dislocation loop. The dislocation loop is identified with location where the phase field value changes smoothly and thus has finite dislocation core width. As shown in Fig. 2(b), the initial phase field in (001) cross-section presents the initial circular dislocation loop with finite dislocation core width. While, the evolution of phase fields are different by using the micro-scale crystal plasticity model and higher-order theory. Our micro-scale crystal plasticity model captures the process of dislocation loop shrinking. As shown in Fig. 2(c),

272

L. Wang et al. / International Journal of Plasticity 81 (2016) 267e283

Fig. 2. (a) The geometry of a dislocation loop in a 3D domain. (b) The initial phase field in (001) cross-section in both the micro-scale crystal plasticity model and higher-order theory; (c) The phase field in (001) cross-section at time step ¼ 150 in micro-scale crystal plasticity model; (d) The phase field in (001) cross-section at time step ¼ 150 in higher-order theory.

the phase field in (001) cross-section at time step ¼ 150 shows the loop shrinking to an elliptical shape elongated along the x direction when using micro-scale crystal plasticity model. However, the higher-order theory doesn't have the ability to capture the dislocation evolution. Fig. 2(d) shows the distribution of phase field in (001) cross-section at time step ¼ 150 when using the higher-order theory. 2.5. Coupled balance equations By applying Gauss theorem, Eq. (15) is further expressed as

Z "

ðV$s þ fÞ$du_ þ

V

X

# Z X  _ a þ 2Bðma  VÞ2 fa  Ap sinð2pfa Þ df_ a dV  ðn$s  tÞ$dudS _  : s  k f εdis a

a

Z 

S

a

(17)

½n$ð2Bma  Vfa  ma Þ  qa df_ a dS ¼ 0

S

Since the choices of all admissible rates of the displacement field du_ and phase fields df_ a are arbitrary, the following coupled balance equations are obtained,

V$s þ f ¼ 0

(18)

L. Wang et al. / International Journal of Plasticity 81 (2016) 267e283

2 _ εdis a : s  kfa þ 2Bðma  VÞ fa  Ap sinð2pfa Þ ¼ 0

273

(19)

and subjected to the boundary conditions,

u ¼ u0

on Su

(20)

n$s ¼ t

on St

(21)

fa ¼ fa

on Sfa

(22)

n$ð2Bma  Vfa  ma Þ ¼ qa

on Sqa

(23)

where n is the outward unit normal vector of the boundary surface. In addition, Eq. (19) is supplemented with initial conditions for phase fields.

fa ðr; 0Þ ¼ fa0 ðrÞ r2V

(24)

The initial phase fields can be used to model pre-existing dislocations, which is discussed in Section 3. Compared with the Khachaturyan-type phase field models (Khachaturyan, 1983; Chen and Khachaturyan, 1991; Wang et al., 1993) in which the stress field is expressed as a function of the inelastic strain through the exact Green's function, the stress field is directly calculated by solving the boundary value problem expressed by Eqs. (18)e(24) in the proposed model. It is very convenient to deal with the complex structures or complex boundary conditions where the analytical Green's function solution is not available. Besides, the elastic modulus mismatch in heteroepitaxial structures is easy to be treated without additional complications. 2.6. Finite element discretization By exploiting the variational structure of the proposed model Eq. (15), the finite element discretization can be achieved. In the Galerkin method, the explicit representations of u(r,t) and fa(r,t) in terms of base functions and nodal variables are as follows,

ui ðr; tÞ ¼

nu X

N I ðrÞuIi ðtÞ

(25)

I¼1

fa ðr; tÞ ¼

nf X

jJa ðrÞfJa ðtÞ

(26)

J¼1

where nu and nf are the dimensions of discrete spaces, NI(r) and jJa ðrÞ are the standard shape functions, i is the spatial degreeJ of-freedom number, uIi ðtÞ and fa ðtÞ are the nodal values of ui(r,t) and fa(r,t), respectively. When Eq. (25) and Eq. (26) are substituted into the weak formulation Eq. (15) and the integration is performed, the discrete FEM formulation is achieved. The following conventions

o n n T fug ¼ fðu1 ÞT ; /; ðuna ÞT gT ; ffa g ¼ f1a ; /; faf ; NI ¼ IN I

(27)

and Voigt notation are used. Thus, the matrix BI and BJa can be defined as the following definitions,

2

I N;x 6 6 0 6 6 0 6 BI ¼ 6 I 6 N;y 6 6 0 4 I N;z

0 I N;y 0 I N;x I N;z 0

3 0 7 2 J 3 0 7 7 ja;x I 7 N;z 7 J 6 7 7; Ba ¼ 4 jJa;y 5 0 7 7 jJa;z I 7 N;y 5 I N;x

So the matrix form of discrete FEM formulation can be achieved,

(28)

274

L. Wang et al. / International Journal of Plasticity 81 (2016) 267e283

X ½Kauf ffa g ¼ ff ext g

½Kuu fug 

(29)

a

½Kbuf T fug þ

X ba     ½Kff ffa g þ ½Kbff  fb þ ½Xb  f_ b þ fQ bcryst g ¼ fQ bext g

(30)

a

where

Z

ðKuu ÞIJ ¼

ðBI ÞT ½CBJ dU0

(31)

J ðBI ÞT ½C½εdis a ja dU0

(32)

T dis J jIb ½εdis b  ½C½εa ja dU0

(33)

U0

Z

ðKauf ÞIJ ¼

U0

Z

IJ ðKba ff Þ ¼

U0

Z

ðKbff ÞIJ ¼ 2B

ðBIb ÞT ½I  mb mTb BJb dU0

(34)

U0

Z

ðXb ÞIJ ¼ k

jIb jJb dU0

(35)

U0

ðQ bcryst ÞI

0

Z ApjIb

¼

sin@2p

ðf ext ÞI ¼

NI tdG0 þ

Z

(36)

Z

NI fdU0

(37)

U0

Gt0

ðQ bext ÞI ¼

1 J J jb fb AdU0

J¼1

U0

Z

nf X

jIb qb dG0

(38)

Gqb 0

At time t, the integral evaluations in Eqs. (31)e(38) are performed once. By employing the forward Euler time integration algorithm in Eq. (30), {fb} at tþDt: ffb g0tþDt is first estimated. Then equilibrium iterations are performed to satisfy Eq. (29) at  correcting {fb} to ffb gi tþDt, while simultaneously to satisfy Eq. (30). The iteration i þ 1 is assumed to be converged when tþDt  J iþ1 J i  maxffb gtþDt  ffb gtþDt   106 . J

3. Computational demonstrations In this section, three examples are given to demonstrate the accuracy of the proposed micro-scale crystal plasticity model. To model a pre-existing dislocation, the phase field of relevant slip system is initially defined as 1 on the slip plane where the dislocation has slipped and transformed to 0 gradually across the slip plane. The transition width is at least of several grids to keep the numerical stability. It is chosen as one grid on both sides of the slip plane in the examples below. The quadratic Lagrange element is used for displacement field and the linear Lagrange element is used for phase field to keep compatible. The examples below show that they are sufficient to the model. In the simulation, the mesh size l is always greater than the real interplanar distance of slip planes d. So ε0 ¼ b/l should be chosen as the typical shear strain instead of ε0 ¼ b/d to maintain the accurate strain field. The dimensionless parameters A* ¼ 2 A=mðε0 Þ which controls the crystalline energy, B* ¼ B=mðε0 Þ2 l2 which controls the gradient energy and Dt * ¼ mðε0 Þ2 Dt=k which decides the stable time step are defined. A* ¼ 1/2p2, B* ¼ 0.07 and Dt* ¼ 0.1 provide reasonable values for the dislocation systems in this paper. 3.1. Dislocation near a free surface A screw dislocation near a free surface of a semi-infinite domain is first considered. The free surface is located at y ¼ 0, and the semi-infinite domain occupies the domain y > 0, as shown in Fig. 3(a). The infinitely long screw dislocation of Burgers vector b is along the positive z-axis and located at a distance of L ¼ 100b from the free surface. The slip plane is parallel to the

L. Wang et al. / International Journal of Plasticity 81 (2016) 267e283

275

Fig. 3. (a) A screw dislocation near a free surface, including the numerical simulation subdomain ABCD. (b) The scaled stress syz/(mb/l) in the grey region labeled in (a).

free surface. The analytical solution of stress fields can be obtained by superposing the infinite domain analytical solution of the same screw dislocation and that of an image screw dislocation (Hirth and Lothe, 1982),

sxz syz

" # mb ðy  LÞ ðy þ LÞ ¼  2p x2 þ ðy  LÞ2 x2 þ ðy þ LÞ2 " # mb x x ¼  2p x2 þ ðy  LÞ2 x2 þ ðy þ LÞ2

(39)

sxy ¼ sxx ¼ syy ¼ szz ¼ 0 where m is shear modulus. The problem is solved on a subdomain ABCD with dimensions 5000b  5000b  10b, as plotted in Fig. 3(a), which can approximately be regarded as a semi-infinite domain. The periodic boundary condition is applied on the surfaces z ¼ 0 and z ¼ 10b to simulate the infinite length along z-direction. The displacements in all three directions on the surface y ¼ 5000b are constrained to restrict the rigid body motion. On the other three boundary surfaces, the free boundary condition is applied. A hexahedral mesh with element size of l ¼ 10b is used. The elastic modulus is E ¼ 26 GPa and Poisson's ratio is n ¼ 0.3. The scaled stress field syz/(mb/l) obtained by this method in an 80l  40l portion which is labeled by the grey region in Fig. 3(a) is plotted in Fig. 3(b). To make quantitative comparisons of syz/(mb/l) between the numerical and analytical solutions, the results along the line in the slip plane y ¼ 10l and along the line that is of one-grid distance from the slip plane y ¼ 9l are plotted in Fig. 4(a) and (b), respectively. It can be found that the numerical solutions agree well with the analytical solutions outside the dislocation core. The linear elastic analytical solutions are singular within the dislocation core, while the numerical solutions are not. 3.2. Dislocation in an anisotropic material In order to verify the accuracy of the method to anisotropic materials, a screw dislocation in an infinite anisotropic material is considered, as shown in Fig. 5(a). The infinitely long screw dislocation of Burgers vector b is along the positive z-axis and

Fig. 4. Quantitative comparisons of syz/(mb/l) between the numerical and analytical solutions; (a) along the line in the slip plane y ¼ 10l; (b) along the line that is of one-grid distance from the slip plane y ¼ 9l.

276

L. Wang et al. / International Journal of Plasticity 81 (2016) 267e283

Fig. 5. (a) A screw dislocation in an anisotropic material, including the numerical simulation subdomain ABCD. (b) The scaled stress syz/(mb/l) in the grey region labeled in (a).

located in the middle of the domain. The slip plane is parallel to x-axis. The stress field has an analytical solution (Hirth and Lothe, 1982). The same geometry, mesh and boundary conditions are used as those in Section 3.1. The material selected is facecentered cubic anisotropic copper and elastic constants are c11 ¼168.4 GPa, c12 ¼ 121.4 GPa and c44 ¼ 75.4 GPa when the coordinate axes coincide with the cubic axes, which are taken from Hirth and Lothe (1982). Since the x, y and z axes are along ½121, [111] and ½101 directions here, the elastic constant matrix needs to be transformed. The stresses are normalized by mb/l, where m ¼ 54.6 GPa is Voigt average elastic constant (Hirth and Lothe, 1982). The scaled shear stress syz/(mb/l) obtained by this method in a 160l  80l portion which is labeled by the grey region in Fig. 5(a) is plotted in Fig. 5(b). Compared with Fig. 3(b), the contour plot is not symmetric due to elastic anisotropy. To make quantitative comparisons of syz/(mb/l) between the numerical and analytical solutions, the results along the line in the slip plane y ¼ 0 and along the line that is of one-grid distance from the slip plane y ¼ l are plotted in Fig. 6(a) and (b), respectively. The numerical solutions agree well with the analytical solutions outside the dislocation core as in the isotropic case.

3.3. Dislocation near a bimaterial interface As shown in Fig. 7(a), an edge dislocation in the vicinity of a bimaterial interface is considered. The bimaterial interface is located at x ¼ 0. The subdomain occupied x > 0 and the subdomain occupied x < 0 are both semi-infinite. The infinitely long edge dislocation of Burgers vector b is along the positive z-axis and located at a distance of L ¼ 95b from the interface. The slip plane is perpendicular to the bimaterial interface. The analytical solution was first given by Head (1953), further clarification was made by Lubarda (1997). It can be regarded as a plane strain problem and solved on a subdomain ABCD with dimensions 10,000b10,000b, as shown in Fig. 7(a). The displacements in all three directions on the bottom line y ¼ 5000b are constrained to restrict the rigid body motion. On the other boundary lines, the free boundary condition is applied. A structured quadrilateral mesh with element size of l ¼ 10b is used. In the subdomain x > 0, elastic modulus E1 ¼ 31 GPa and Poisson's ratio n1 ¼ 0.276; in the subdomain x < 0, elastic modulus E2 ¼ 94.7 GPa and Poisson's ratio n2 ¼ 0.276.

Fig. 6. Quantitative comparisons of syz/(mb/l) between the numerical and analytical solutions; (a) along the line in the slip plane y ¼ 0; (b) along the line that is of one-grid distance from the slip plane y ¼ l.

L. Wang et al. / International Journal of Plasticity 81 (2016) 267e283

277

Fig. 7. (a) An edge dislocation near a bimaterial interface, including the numerical simulation subdomain ABCD. (b) The scaled stress sxy/(m1b/l) in the grey region labeled in (a).

The scaled shear stress sxy/(m1b/l) obtained by this method in a 30l  30l portion which is labeled by the grey region in Fig. 7(a) is plotted in Fig. 7(b). m1 is the shear modulus of the subdomain x > 0. It can be found that sxy/(m1b/l) is continuous across the bimaterial interface as expected. To make quantitative comparisons of sxy/(m1b/l) between the numerical and analytical solutions, the results along the line in the slip plane y ¼ 0 and along the line that is of two-grid distance from the slip plane y ¼ 2l are plotted in Fig. 8(a) and (b), respectively. The numerical solutions agree well with the analytical solutions outside the dislocation core. The linear elastic analytical solutions are singular within the dislocation core, while bounded solutions are predicted in the calculation. 4. Applications in heteroepitaxial structures In this section, two examples are given to demonstrate the applicability of the micro-scale crystal plasticity model for modeling dislocations in typical heteroepitaxial structures. 4.1. Critical shell thickness of heteroepitaxial coreeshell nanopillars Heteroepitaxial coreeshell nanopillars have recently received much attention because of their important applications in electronics, optoelectronics and solar cells. The coreeshell interface is coherent and free of dislocations when the shell is initially deposited. As the shell layer thickens, the stored strain energy due to the lattice misfit strain increases. When the shell is above a critical thickness, a stable misfit dislocation will form since the strain energy due to the lattice misfit strain exceeds

Fig. 8. Quantitative comparisons of sxy/(m1b/l) between the numerical and analytical solutions; (a) along the line in the slip plane y ¼ 0; (b) along the line that is of two-grid distance from the slip plane y ¼ 2l.

278

L. Wang et al. / International Journal of Plasticity 81 (2016) 267e283

the dislocation formation energy. The critical shell thickness of Ge/Si coreeshell nanopillars is predicted by employing the micro-scale crystal plasticity model developed in Section 2. In addition to energetic considerations, it is also important to know the kinetic processes that how the misfit dislocation comes into existence to predict the critical shell thickness. The glide-relaxation mechanism is pursued here for Ge/Si coreeshell nanopillars, which is consistent with the experimental studies (Goldthorpe et al., 2008; Dayeh et al., 2013). For Ge/Si coreeshell nanopillars which have diamond crystal structures, aSi/2〈110〉 perfect dislocations nucleate from the surface of the shell and slip on the inclined {111} slip planes to deposit elliptical dislocation loops at the coreeshell interface. The slip system associated with the smallest critical shell thickness is the most likely activated in actuality. In the calculation, [111] growth direction is considered. For convenience, a coordinate system with x, y and z axes along ½101, ½121 and [111] directions respectively is defined, as shown in Fig. 9(a). The lattice constants of Si and Ge crystals are aSi ¼ 0.5433 nm and aGe ¼ 0.5660 nm, respectively (Schwarz, 1999). In a Ge/Si coreeshell nanopillar, the lattice misfit strain of Si with respect to Ge is εm ¼ (aGeaSi)/aSi. To calculate the elastic field due to the lattice misfit strain, the initial strain field εinitial¼εmI is introduced in the shell to match the lattice of the shell to the core. So the total strain field is given by ε ¼ Vsuþεinitial. Then the elastic field can be obtained by solving the balance equation (18) and boundary conditions (20)e(21). In this section, the anisotropic elastic constants are used. c11 ¼128.4 GPa, c12 ¼ 48.2 GPa, c44 ¼ 66.7 GPa correspond to the core material Ge and c11 ¼166.2 GPa, c12 ¼ 64.4 GPa, c44 ¼ 79.8 GPa correspond to the shell material Si, when the coordinate axes coincide with the cubic axes (Freund and Suresh, 2003). The elastic field due to the lattice misfit strain is first calculated for the coreeshell nanopillar which is shown in Fig. 9(a). The core radius is 15 nm and the shell thickness is 5 nm. The height of the coreeshell nanopillar is set to 240 nm. An FEM mesh generated by the tetrahedral elements is used. To identify the most likely nucleation sites, the resolved shear stress (RSS) computed by the proposed method is plotted in Fig. 9(b), for a dislocation with b ¼ aSi/2[110] Burgers vector and ð111Þ slip plane. The RSS with a positive value promotes the nucleation and propagation of the dislocation. While, the RSS with a zero or negative value will suppress the nucleation and propagation of the dislocation. So the dislocation is expected to nucleate from the location where the RSS is with the maximum positive value and to evolve in the shell by inspecting Fig. 9(b). To determine the critical shell thickness, the reduction in stored energy of coreeshell nanopillar DE from the dislocation-free state to the misfit dislocation state is considered. The misfit dislocation tends to form when the following criterion is satisfied,

DE ¼ Em  E0  0

(40)

where Em is the stored energy of coreeshell nanopillar with the presence of misfit dislocation loop at the coreeshell interface and E0 is the stored energy due to lattice misfit strain only. The critical shell thickness is then evaluated when the reduction in stored energy DE equals to zero. For the Ge core radius is 15 nm, the reduction in stored energy DE with the variation of the Si shell thickness from 2 nm to 5 nm is plotted in Fig. 10(a). It shows that the results are convergent with an average mesh size of l ¼ 1.5 nm and the critical shell thickness is 3.5 nm when the core radius is 15 nm. It is consistent with the experimental result for Ge/Si coreeshell nanopillars, which is ~3 nm (Dayeh et al., 2013). And the result predicted by the proposed method is closer to the experimental result than the analytical prediction made by Chu et al. (2013), which is 3.7 nm under isotropic assumption of materials properties.

Fig. 9. (a) The geometry and mesh of coreeshell nanopillar with [111] growth direction, the core radius is 15 nm and the shell thickness is 5 nm. (b) Resolved shear stress (RSS) in a ð111Þ plane crossing the nanopillar in (a) for the [110]ð111Þ slip system.

L. Wang et al. / International Journal of Plasticity 81 (2016) 267e283

(b) Mesh size: 1.5 nm Mesh size: 3 nm

1.0 0.5 0.0 -0.5 -1.0 -1.5 -2.0

Core radius=15 nm 2.0

2.5

3.0 3.5 4.0 4.5 Shell thickness (nm)

5.0

Critical shell thickness (nm)

-16

Reduction of energy (10 J)

(a) 1.5

279

Simulation results Experimental data Analytical predictions

8

6

4

2

5

10 15 20 Core radius (nm)

25

Fig. 10. (a) The reduction in stored energy DE from the dislocation-free state to the misfit dislocation state, for the Ge core radius is 15 nm and the Si shell thickness changes from 2 nm to 5 nm. (b) The dependence of the critical shell thickness on the core radius predicted by the proposed method, the experimental and analytical critical shell thicknesses are also plotted for comparison.

The predicted critical shell thicknesses with respect to different core radiuses for Ge/Si coreeshell nanopillars are shown in Fig. 10(b). It suggests that the variation of the critical shell thickness depends on the core radius, and the reduction in core radius can increase the critical shell thickness. When the core radius is sufficiently small, misfit dislocation will not nucleate any more. While, the critical shell thickness tends to a constant as the core radius increases. These trends are in good agreement with the theoretical results (Chu et al., 2013), as shown in Fig. 10(b). 4.2. Dislocations in heteroepitaxial thin films Like heteroepitaxial coreeshell nanopillars, the dislocation-free film cannot be grown with arbitrary thickness on substrate and misfit dislocation will form above a critical thickness. Here, the proposed micro-scale crystal plasticity model is employed to simulate dislocation behaviors in Si1-xGex heteroepitaxial thin film grown on Si substrate. In particular, the operation of a dislocation source that can generate dislocation loops is considered. The same lattice constants of Si and Ge crystals are used as those in Section 4.1. Then the stress-free lattice constant of Si 1xGex can be well approximated by the Vegard's law, a SiGe ¼ aSi(1x) þ a Gex, where x represents the atomic fraction of Ge in the Si1-x Gex . In the simulation, x ¼ 0.24 is assumed. So the internal stresses will be generated in the heteroepitaxial thin films due to the lattice misfit and the lattice misfit strain is εm ¼ (aSiGea Si)/aSi ¼ 0.0418x. If the film thickness is larger than the critical thickness, this internal stresses can be relieved by the nucleation and evolution of dislocations. Here we use the thermal expansion process to achieve the internal stresses caused by the misfit strain (Cui et al., 2015a). Both Si and Ge are assumed isotropic, with Young's modulus ESi ¼ 130 GPa (Shen, 2008) and EGe ¼ 102.5 GPa (Neuberger, 1971), respectively, as well as the same Poisson ratio nSi¼nGe ¼ 0.28. The Young's modulus and Poisson ratio of Si1-xGex are approximated by the Vegard's law, ESiGe ¼ ESi(1x) þ EGex and nSiGe¼nSi(1x) þ nGex, respectively. Then the equivalent thermal expansion coefficients of the film and substrate are given as follows (Freund and Suresh, 2003),

af ¼ εm

hs Ms hf Mf þ hs Ms

hf Mf as ¼ εm hf Mf þ hs Ms

(41)

where hf and hs are the thicknesses of the film and substrate, Mf ¼ ESiGe/(1nSiGe) and Ms ¼ ESi/(1nSi) are the biaxial modulus of the film and substrate, respectively. Then the internal stresses caused by the misfit strain are equal to those caused by the thermal expansion in one unit. In the simulation, a side length of 3 mm is used for both the film and substrate. The thicknesses of the film and substrate are set to hf ¼ 0.3 mm and hs ¼ 0.6 mm, respectively. The out-of-plane displacement on the bottom surface of substrate is constrained and a lateral periodic boundary condition is applied to exclude the lateral bend. A diamond cubic epitaxial film system with (100) interface is considered. Supposing there is a source in (111) slip plane that can generate dislocation loops of 0.15 mm with Burgers vector b of magnitude 0.236 nm along ½101 in the middle of the film, it naturally multiplies under the effect of internal stress caused by the misfit strain and image force from the free surface, without any a priori constraint. An FEM mesh generated by the tetrahedral elements is used with an average element size of l ¼ 60b, as shown in Fig. 11. The evolution of a dislocation source in (111) cross-section of the epitaxial film at various time moments is illustrated in Fig. 12(a)-(d). As the dislocation loop expands and encounters the filmesubstrate interface and the free surface, it forms two threading arms and deposits a misfit dislocation segment at the filmesubstrate interface. The threading arms near the free

280

L. Wang et al. / International Journal of Plasticity 81 (2016) 267e283

Fig. 11. The geometry and mesh of the epitaxial film with (100) interface.

surface try to stay perpendicular to the free surface because of the image force produced by the free surface, which is accurately captured by the proposed model. The misfit segment is prevented from entering the substrate due to the internal stress of opposite sign in the substrate. The slip-steps on the free surface of the film produced by the passage and exit of the threading arms are described in Fig. 13. The deformation is enlarged 300 times and the scaled displacement nephogram of z direction w/b is also plotted in Fig. 13. As a comparison, the evolution of the same dislocation source in a periodic multilayer is considered. The same model geometry is used as the one of epitaxial film, except that the periodic boundary condition is applied in the vertical direction. The simulated evolution of the dislocation source in (111) cross-section of the multilayer at various time moments is illustrated in Fig. 12(e)e(h). In the multilayer, there is no free surface and the dislocation evolution is not affected by the image force from the free surface. As the dislocation loop expands and encounters the filmesubstrate interfaces, it forms two threading arms and deposits two misfit dislocation segments at the filmesubstrate interfaces. The misfit segments are also prevented from entering the substrates due to the internal stress of opposite sign in the substrates. The threading arms bow out symmetrically, which is different from the situation in epitaxial film. The results are consistent with previous studies (Wang et al., 2003; Cui et al., 2015a). By comparing the two cases, it can be found that the threading arms in the epitaxial film move faster than those in the multilayer, as illustrated in Fig. 12. It can be explained from an energetic standpoint. The stored elastic energy is the same in the two systems initially. While, the decrement of elastic energy is larger in the epitaxial film which indicates that the velocity of dislocation slip is larger in the epitaxial film, as shown in Fig. 14(a). As shown in Fig. 14(b), the crystalline energy and gradient energy are both higher in the multilayer at the same time step, which is caused by the formation of an extra misfit dislocation segment in the multilayer. It can be noted that the stress in the film can be relieved by dislocation

Fig. 12. The simulated evolution of a dislocation source in (111) cross-section of the epitaxial film and the multilayer at various time moments, as well as Mises stress: (a) and (e) Time step ¼ 0; (b) and (f) Time step ¼ 200; (c) and (g) Time step ¼ 1130; (d) and (h) Time step ¼ 1550.

L. Wang et al. / International Journal of Plasticity 81 (2016) 267e283

281

Fig. 13. The scaled displacement nephogram of z direction w/b of epitaxial film at time step ¼ 1550 and the deformation is enlarged 300 times.

slip in both cases. While, the decreased amount of stress is smaller in the epitaxial film, illustrating that the degradation of device performance is smaller when use epitaxial film than multilayer. However, the stress at the filmesubstrate interface is increased when misfit dislocations pile up at the interface. Such locations may be the preferred nucleation sites for cracks.

5. Conclusions Under a thermodynamically consistent framework, we have developed a micro-scale crystal plasticity model based on phase field theory. Compared with the widely used Khachaturyan-type phase field model of dislocations, the proposed model has some great advantages in the analysis of dislocations in complex heteroepitaxial structures. It can be used for complex structures or complex boundary conditions where the analytical Green's function solution is not available. And the elastic modulus mismatch in heteroepitaxial structures is easy to be treated without additional complications. Besides, with numerical implementation by FEM, it is flexible to deal with finite deformation in heteroepitaxial structures at micro-scale. Compared with conventional crystal plasticity model, the proposed model can be used at a smaller size scale where a single dislocation behavior is important. The accuracy of the proposed model is studied by comparing with the analytical solutions in three problems: a screw dislocation near a free surface, a screw dislocation in infinite anisotropic material and an edge dislocation near a bimaterial interface. Since FEM is used to solve the boundary value problem, the mesh should be sufficiently refined around the dislocation core. The numerical solutions agree well with the analytical solutions outside the dislocation core. The linear elastic analytical solutions are singular within the dislocation core, while bounded solutions are predicted in the simulation. The critical shell thickness of heteroepitaxial coreeshell nanopillars is predicted by employing the proposed model. The result is consistent with the experimental result and more accurate than the analytical prediction. The results suggest that the variation of the critical shell thickness depends on the core radius, and the reduction in core radius can increase the critical shell thickness. The model also captures the fundamental physical trends that misfit dislocation will not nucleate when the core radius is sufficiently small and the critical shell thickness tends to a constant as the core radius increases. The evolution of a dislocation source in heteroepitaxial thin films is also simulated by employing the proposed model. It shows that the image forces produced by the free surface and filmesubstrate interface can be properly estimated. At the same time, the structure deformation is also calculated and the slip-steps on the free surface of film produced by the passage and exit of dislocations are well captured. At last, the reason for why the threading arms in the epitaxial film move faster than those in the multilayer is explained by comparing the energy of the two systems.

(a) 3.25 Eel Eel_multilayer

6.00 -14

3.23 3.22

Ecryst Egrad Ecryst_multilayer Egrad_multilayer

7.00 Energy (10 J)

-11

Energy (10 J)

3.24

(b) 8.00

5.00 4.00 3.00 2.00 1.00

3.21

0

200 400 600 800 1000 12001400 1600 Time step

0.00

0

200 400 600 800 1000 12001400 1600 Time step

Fig. 14. (a) The elastic energy of two systems. (b) The crystalline energy and gradient energy of two systems.

282

L. Wang et al. / International Journal of Plasticity 81 (2016) 267e283

Next, the micro-scale crystal plasticity model developed in this paper can be improved on in terms of large deformation and interface energy. Acknowledgments This work is supported by the National Natural Science Foundation of China under Grant Nos. 11132006 and 11302115, National Basic Research Program of China under Grant No. 2015CB351900. References Abdolrahim, N., Zbib, H.M., Bahr, D.F., 2014. Multiscale modeling and simulation of deformation in nanoscale metallic multilayer systems. Int. J. Plast. 52, 33e50. Acharya, A., 2001. A model of crystal plasticity based on the theory of continuously distributed dislocations. J. Mech. Phys. Solids 49, 761e784. Amodeo, R.J., Ghoniem, N.M., 1990. Dislocation dynamics. I. A proposed methodology for deformation micromechanics. Phys. Rev. B 41, 6958e6967. Belytschko, T., Gracie, R., 2007. On XFEM applications to dislocations and interfaces. Int. J. Plast. 23, 1721e1738. Bennett, S.E., 2010. Dislocations and their reduction in GaN. Mater. Sci. Technol. 26, 1017e1028. Chen, L.Q., Khachaturyan, A.G., 1991. Computer simulation of structural transformations during precipitation of an ordered intermetallic phase. Acta Metall. Mater 39, 2533e2551. Chu, H.J., Zhou, C.Z., Wang, J., Beyerlein, I.J., 2013. An analytical model for the critical shell thickness in core/shell nanowires based on crystallographic slip. J. Mech. Phys. Solids 61, 2147e2160. Cui, Y.N., Lin, P., Liu, Z.L., Zhuang, Z., 2014. Theoretical and numerical investigations of single arm dislocation source controlled plastic flow in FCC micropillars. Int. J. Plast. 55, 279e292. Cui, Y.N., Liu, Z.L., Zhuang, Z., 2015a. Quantitative investigations on dislocation based discrete-continuous model of crystal plasticity at submicron scale. Int. J. Plast. 69, 54e72. Cui, Y.N., Liu, Z.L., Zhuang, Z., 2015b. Theoretical and numerical investigations on confined plasticity in micropillars. J. Mech. Phys. Solids 76, 127e143. Dayeh, S.A., Tang, W., Boioli, F., Kavanagh, K.L., Zheng, H., Wang, J., Mack, N.H., Swadener, G., Huang, J.Y., Miglio, L., Tu, K.-N., Picraux, S.T., 2013. Direct measurement of coherency limits for strain relaxation in heteroepitaxial core/shell nanowires. Nano Lett. 13, 1869e1876. El-Awady, J.A., Bulent Biner, S., Ghoniem, N.M., 2008. A self-consistent boundary element, parametric dislocation dynamics formulation of plastic flow in finite volumes. J. Mech. Phys. Solids 56, 2019e2035. Fivel, M.C., Canova, G.R., 1999. Developing rigorous boundary conditions to simulations of discrete dislocation dynamics. Modell. Simul. Mater. Sci. Eng. 7, 753. Freund, L.B., 2000. The mechanics of electronic materials. Int. J. Solids Struct. 37, 185e196. Freund, L.B., Suresh, S., 2003. Thin Film Materials. Cambridge University Press, Cambridge. Fu, Y., Du, H., Huang, W., Zhang, S., Hu, M., 2004. TiNi-based thin films in MEMS applications: a review. Sens. Actuators, A 112, 395e408. Gao, Y., Liu, Z.L., You, X.C., Zhuang, Z., 2010. A hybrid multiscale computational framework of crystal plasticity at submicron scales. Comp. Mater. Sci. 49, 672e681. Giessen, E.V.D., Needleman, A., 1995. Discrete dislocation plasticity: a simple planar model. Modell. Simul. Mater. Sci. Eng. 3, 689e735. Goldthorpe, I.A., Marshall, A.F., McIntyre, P.C., 2008. Synthesis and strain relaxation of Ge-core/Si-shell nanowire arrays. Nano Lett. 8, 4081e4086. Gracie, R., Oswald, J., Belytschko, T., 2008. On a new extended finite element method for dislocations: core enrichment and nonlinear formulation. J. Mech. Phys. Solids 56, 200e214. Gracie, R., Ventura, G., Belytschko, T., 2007. A new fast finite element method for dislocations based on interior discontinuities. Int. J. Numer. Methods. Eng. 69, 423e441. Groh, S., Devincre, B., Kubin, L.P., Roos, A., Feyel, F., Chaboche, J.L., 2003. Dislocations and elastic anisotropy in heteroepitaxial metallic thin films. Philos. Mag. Lett. 83, 303e313. Gurtin, M.E., 2002. A gradient theory of single-crystal viscoplasticity that accounts for geometrically necessary dislocations. J. Mech. Phys. Solids 50, 5e32. Han, X., Ghoniem, N.M., 2005. Stress field and interaction forces of dislocations in anisotropic multilayer thin films. Philos. Mag. 85, 1205e1225. Head, A.K., 1953. Edge dislocations in inhomogeneous Media. Proc. Phys. Soc. B 66, 793e801. Hirth, J.P., Lothe, J., 1982. Theory of Dislocations. Wiley, New York. Hu, S.Y., Chen, L.Q., 2001. Solute segregation and coherent nucleation and growth near a dislocation - a phase-field model integrating defect and phase microstructures. Acta Mater 49, 463e472. Hurtado, D.E., Ortiz, M., 2012. Surface effects and the size-dependent hardening and strengthening of nickel micropillars. J. Mech. Phys. Solids 60, 1432e1446. Khachaturyan, A.G., 1967. Some questions concerning theory of phase transformations in solids. Sov. Phys. Solid State 8, 2163. Khachaturyan, A.G., 1983. Theory of Structural Transformations in Solids. John Wiley and Sons, New York. Khachaturyan, A.G., Shatalov, G.A., 1969. Elastic-interaction potential of defects in a crystal. Sov. Phys. Solid State 11, 118. Koslowski, M., LeSar, R., Thomson, R., 2004. Dislocation structures and the deformation of materials. Phys. Rev. Lett. 93. Kuroda, M., Tvergaard, V., 2008. On the formulations of higher-order strain gradient crystal plasticity models. J. Mech. Phys. Solids 56, 1591e1608. Lei, L., Marin, J.L., Koslowski, M., 2013. Phase-field modeling of defect nucleation and propagation in domains with material inhomogeneities. Modell. Simul. Mater. Sci. Eng. 21. Lemarchand, C., Devincre, B., Kubin, L.P., 2001. Homogenization method for a discrete-continuum simulation of dislocation dynamics. J. Mech. Phys. Solids 49, 1969e1982. Liu, Y., Chen, Y., Yu, K.Y., Wang, H., Chen, J., Zhang, X., 2013. Stacking fault and partial dislocation dominated strengthening mechanisms in highly textured Cu/Co multilayers. Int. J. Plast. 49, 152e163. Liu, Z.L., Liu, X.M., Zhuang, Z., You, X.C., 2009. A multi-scale computational model of crystal plasticity at submicron-to-nanometer scales. Int. J. Plast. 25, 1436e1455. Lubarda, V.A., 1997. Energy analysis of dislocation arrays near bimaterial interfaces. Int. J. Solids Struct. 34, 1053e1073. Nabarro, F.R.N., 1951. The synthesis of elastic dislocation fields. Philos. Mag. 42, 1224e1231. Neuberger, M., 1971. Volume 5_ Group IV Semiconducting Materials- Handbook of Electronic Materials. Plenum Press, London. Oswald, J., Gracie, R., Khare, R., Belytschko, T., 2009. An extended finite element method for dislocations in complex geometries: thin films and nanotubes. Comput. Methods Appl. Mech. Eng. 198, 1872e1886. Panda, D., Tseng, T.-Y., 2013. Growth, dielectric properties, and memory device applications of ZrO2 thin films. Thin Solid Films 531, 1e20. Puri, S., Das, A., Acharya, A., 2011. Mechanical response of multicrystalline thin films in mesoscale field dislocation mechanics. J. Mech. Phys. Solids 59, 2400e2417. Rodney, D., Le Bouar, Y., Finel, A., 2003. Phase field methods and dislocations. Acta Mater 51, 17e30. Roy, A., Acharya, A., 2005. Finite element approximation of field dislocation mechanics. J. Mech. Phys. Solids 53, 143e170. Schwarz, K.W., 1999. Simulation of dislocations on the mesoscopic scale. II. Application to strained-layer relaxation. J. Appl. Phys. 85, 120e129.

L. Wang et al. / International Journal of Plasticity 81 (2016) 267e283

283

Shen, C., Wang, Y., 2003. Phase field model of dislocation networks. Acta Mater 51, 2595e2610. Shen, Y.L., 2008. Externally constrained plastic flow in miniaturized metallic structures: a continuum-based approach to thin films, lines, and joints. Prog. Mater Sci. 53, 838e891. Takahashi, A., Ghoniem, N., 2008. A computational method for dislocationeprecipitate interaction. J. Mech. Phys. Solids 56, 1534e1553.
, A., Devincre, B., Feyel, F., Gatti, R., Groh, S., Jamond, O., Roos, A., 2014. Modelling crystal plasticity by 3D dislocation dynamics and the finite element Vattre method: the discrete-continuous model revisited. J. Mech. Phys. Solids 63, 491e505. Wang, Y.U., Chen, L.Q., Khachaturyan, A.G., 1993. Kinetics of strain-induced morphological transformation in cubic alloys with a miscibility gap. Acta Metall. Mater 41, 279e296. Wang, Y.U., Jin, Y.M., Cuitino, A.M., Khachaturyan, A.G., 2001. Nanoscale phase field microelasticity theory of dislocations: model and 3D simulations. Acta Mater. 49, 1847e1857. Wang, Y.U., Jin, Y.M., Khachaturyan, A.G., 2002. Phase field microelasticity theory and modeling of elastically and structurally inhomogeneous solid. J. Appl. Phys. 92, 1351. Wang, Y.U., Jin, Y.M., Khachaturyan, A.G., 2003. Phase field microelasticity modeling of dislocation dynamics near free surface and in heteroepitaxial thin films. Acta Mater 51, 4209e4223. Wang, Y.U., Jin, Y.M., Khachaturyan, A.G., 2004. Phase field microelasticity modeling of surface instability of heteroepitaxial thin films. Acta Mater 52, 81e92. Wang, Y.U., Jin, Y.M., Khachaturyan, A.G., 2005. Mesoscale modelling of mobile crystal defectsddislocations, cracks and surface roughening: phase field microelasticity approach. Philos. Mag. 85, 261e277. Weygand, D., Friedman, L.H., van der Giessen, E., Needleman, A., 2001. Discrete dislocation modeling in three-dimensional confined volumes. Mater. Sci. Eng. A 309, 420e424. Zbib, H.M., de la Rubia, T.D., 2002. A multiscale model of plasticity. Int. J. Plast. 18, 1133e1163. Zbib, H.M., Overman, C.T., Akasheh, F., Bahr, D., 2011. Analysis of plastic deformation in nanoscale metallic multilayers with coherent and incoherent interfaces. Int. J. Plast. 27, 1618e1639. Zbib, H.M., Rhee, M., Hirth, J.P., 1998. On plastic deformation and the dynamics of 3D dislocations. Int. J. Mech. Sci. 40, 113e127.