Coupled DDD–FEM modeling on the mechanical behavior of microlayered metallic multilayer film at elevated temperature

Journal of the Mechanics and Physics of Solids 85 (2015) 74–97

Contents lists available at ScienceDirect

Journal of the Mechanics and Physics of Solids journal homepage: www.elsevier.com/locate/jmps

Coupled DDD–FEM modeling on the mechanical behavior of microlayered metallic multilayer film at elevated temperature Minsheng Huang a,b, Zhenhuan Li a,b,n a b

Department of Mechanics, Huazhong University of Science and Technology, Wuhan 430074, China Hubei Key Laboratory of Engineering Structural Analysis and Safety Assessment, Luoyu Road 1037, Wuhan 430074, China

a r t i c l e in f o

abstract

Article history: Received 23 December 2014 Received in revised form 15 July 2015 Accepted 15 September 2015 Available online 24 September 2015

To investigate the mechanical behavior of the microlayered metallic thin films (MMMFs) at elevated temperature, an enhanced discrete-continuous model (DCM), which couples rather than superposes the two-dimensional climb/glide-enabled discrete dislocation dynamics (2D-DDD) with the linearly elastic finite element method (FEM), is developed in this study. In the present coupling scheme, two especial treatments are made. One is to solve how the plastic strain captured by the DDD module is transferred properly to the FEM module as an eigen-strain; the other is to answer how the stress field computationally obtained by the FEM module is transferred accurately to the DDD module to drive those discrete dislocations moving correctly. With these two especial treatments, the interactions between adjacent dislocations and between dislocation pile-ups and interphase boundaries (IBs), which are crucial to the strengthening effect in MMMFs, are carefully taken into account. After verified by comparing the computationally predicted results with the theoretical solutions for a dislocation residing in a homogeneous material and nearby a bi-material interface, this 2D-DDD/FEM coupling scheme is used to model the tensile mechanical behaviors of MMMFs at elevated temperature. The strengthening mechanism of MMMFs and the layer thickness effect are studied in detail, with special attentions to the influence of dislocation climb on them. & 2015 Elsevier Ltd. All rights reserved.

Keywords: DDD–FEM coupling scheme Metallic multilayer film High temperature Layer thickness effect Strengthening mechanism Dislocation climb

1. Introduction With rapid progresses in the magnetron sputtering deposition technology, by which two or more metals can be alternately deposited layer by layer with the controllable layer thickness ranging from several to hundreds of nanometers, various microlayered metallic multilayer films (MMMFs) can be synthesized rapidly and expediently. High density interphase boundaries (IBs) endow the MMMFs a series of excellent mechanical properties, including extraordinary strength (Misra et al., 2002; 2005), good morphological stability even at high temperatures (Mara et al., 2004), improved fatigue resistance (Wang et al., 2006) and enhanced irradiation damage resistance (Höchbauer et al., 2005; Misra et al., 2007). Due to these outstanding properties, the MMMFs are increasingly used as multifunctional and structural materials in some rising fields, for example, integrated circuit chip, coating, micro/nano-sensors and actuators, structural or functional modules for the Micro–Electro-Mechanical Systems (MEMS) etc. (Was and Foecke, 1996; Misra et al., 2002; Ghoniem and Han, 2005). However, in materials science, there exists a well-known dilemma, i.e. the enhancement of material strength is usually at n

Corresponding author at: Department of Mechanics, Huazhong University of Science and Technology, Wuhan 430074, China. Fax: þ 86 27 87543501. E-mail address: [email protected] (Z. Li).

http://dx.doi.org/10.1016/j.jmps.2015.09.007 0022-5096/& 2015 Elsevier Ltd. All rights reserved.

M. Huang, Z. Li / J. Mech. Phys. Solids 85 (2015) 74–97

75

the expense of its ductility. How to achieve a marked improvement of material strength without sacrificing much of its ductility is a highly interesting topic for material scientists. A satisfactory solution to this dilemma depends on deep understanding of strengthening and toughening mechanisms of materials. During the past two decades, a considerable amount of research has been done on the strengthening mechanisms of MMMFs; however, most of these works ignored the influence of temperature on the strengthening mechanisms of MMMFs although they usually work at high temperature. How high temperature-induced dislocation climb influences the strengthening mechanism of materials has come into focus recently (Gao et al., 2011, 2013; Davoudi et al., 2012; Ayas et al., 2012, 2014; Ahmed and Hartmaier, 2013; Danas and Deshpande, 2013; Huang et al. 2014a, b; Yang et al., 2015) but few on the MMMFs. A series of experiments repeatedly indicated that the strength of MMMFs is much higher than that predicted by the often-used mixture rule (Keralavarma and Benzerga, 2007), achieving 1/3–1/2 the theoretical strength of constituent materials of MMMFs (Akasheh et al., 2007a, b). This ultra-high strength is believed to stem from the interaction between the lattice dislocations and the inter-phase boundaries (IBs) (Chu et al., 2013). It has been experimentally observed that the dislocation–IB interaction mechanism in the MMMFs alters with changing the layer thickness, resulting in different dependences of the MMMF strength on the individual layer thickness h (Misra et al., 2002, 2005; Zbib et al., 2011; Abdolrahim et al., 2014a). With decreasing layer thickness, the strength of MMMFs first increases gradually, then reaches a maximum plateau and finally levels off or even decreases slightly (Misra et al., 2002, 2005). This strength vs. layer thickness relation of MMMFs can be roughly divided into three regimes where the plastic deformation is governed by three different mechanisms, respectively (Misra et al., 2002, 2005; Akasheh et al., 2007a; Keralavarma and Benzerga, 2007; Zbib et al., 2011; Zhu et al. 2014). In regime I, where the layer thickness h of MMMFs is above dozens of nanometers, the dislocation pile-ups at the IBs can easily be formed in these sufficiently thick metallic layers and govern the plastic deformation of MMMFs, so the strength dependence of MMMFs on the layer thickness h can be depicted by the classical Hall–Petch relation as σ ∝ h−n with n ≈ 1/2 similar to that in polycrystals. In regime II, where h is in the range of dozens of nanometers to several nanometers, long dislocation pile-ups at IBs are difficult to form due to layer thickness restriction but the slip of single hairpinshaped dislocation loop parallel to IBs is still permitted within such confined layer, thus the strength dependence of MMMFs on the layer thickness h can be depicted by the confined layer slip (CLS) model as σ ∝ ln (αh) /h (Matthews and Blakeslee, 1974). In regime III, where the layer thickness h is at several nanometers, the strength reaches a plateau value and in some cases drops to a certain extent, indicating a considerable modification of the strengthening mechanism. At such small layer thickness closing to the width of dislocation core, dislocations can cross IBs and propagate from one layer to the other, thus the strength of MMMFs is controlled by the interfacial barrier strength (IBS) rather than by the layer thickness (Li and Anderson, 2005). Obviously, the dependence of MMMFs strength on the layer thickness is closely related to the strong constraint of abundant IBs to gliding lattice dislocations. Besides, a variety of structural and material factors also play key roles in the strength vs. layer thickness relation of MMMFs (Keralavarma and Benzerga, 2007), including the elastic constant mismatch-induced dislocation image stress, the lattice mismatch-induced coherency stress and interfacial dislocation network, the slip system mismatch-induced cross slip resistance etc. In order to describe quantitatively the layer thickness effect on the strength of MMMFs, several analytical models have been developed and employed, including the dislocation-based interface shear model (Chu et al., 2013), dislocation glide across layers model (Chu and Barnett, 1995), CLS model (Matthews and Blakeslee, 1974), Hall–Petch model (Misra et al., 2005) etc.. Although most of these analytical models can depict the strength vs. layer thickness relation in different regimes mentioned above, they are all based on the ad hoc assumption of the governing strengthening mechanism in respective regimes (Keralavarma and Benzerga, 2007), and the image stress on lattice dislocations due to elastic constant mismatch between two phases is not considered explicitly. In addition, other factors influencing the hardening behaviors of MMMFs, such as dislocation–dislocation interactions and intersections (Misra and Kung, 2001), kinetic effects (Dodson and Tsao, 1987), Peirels friction (Matthews and Crawford, 1970) and slip step formation at surfaces and interfaces (Hoagland et al., 2004), are not included in these models. These simplifications generally render the strength of MMMFs to be underestimated by above analytical models (Akasheh et al., 2007a, b; Zbib et al., 2011). In order to evade these simplifications, the numerical simulation should be adopted. Since atomic scale modeling can consider above factors directly and provide atomic details of deformation mechanisms, many molecular dynamics (MD) simulations have been performed in recent years to study the nature of dislocations, interfaces and their interactions in the MMMFs (Hoagland et al., 2002, 2006; Mastorakos et al., 2009; Shao and Medyanik, 2010; Chu et al., 2013; Shao et al., 2013; Abdolrahim et al., 2014a, b; Zhu et al., 2014; Wang et al., 2014). Although these MD simulations can capture the atomic details of deformation mechanisms, they are computationally very expensive, which makes them unsuitable for modeling of MMMFs with individual layer thickness at the submicron scale or above. In addition, as well-known, MD also cannot deal with practical temporal scale on the order of microsecond or larger (Zbib et al., 2011). In order to model plastic behavior of materials at larger spatial scale and longer temporal scale, discrete dislocation dynamics (DDD), which omits the atomic details of dislocations, has been developed and extensively used (Van der Giessen and Needleman, 1995; Li et al., 2009; Zbib et al., 1998, 2002, 2011; Huang and Li, 2013). By means of a series of short/long range physical rules, which can be obtained from MD, DDD can physically capture the dislocation–dislocation and dislocation– interface/surface interactions thus becomes a popular and effective tool extensively used to capture the size-dependent plastic behavior of micron-sized materials and its underlying dislocation mechanisms. In the last two decades, a series of DDD simulations have been performed to study the layer thickness effect on the mechanical behavior of MMMFs and its underlying strengthening mechanism. By employing 3D-DDD, the interactions between lattice threading dislocation and

76

M. Huang, Z. Li / J. Mech. Phys. Solids 85 (2015) 74–97

interfacial dislocation and their strengthening effect were modeled carefully (Schwarz and Tersoff, 1996; Schwarz,1999; Akasheh et al., 2007a,b; Zbib et al., 2011; Bellou et al., 2011). Although these 3D-DDD modelings show good potential to describe the deformation behavior of MMMFs, enormous computational cost of dealing with the image stress induced by the elastic mismatch between different metallic layers severely limits their scope of investigation (Keralavarma and Benzerga, 2007), and the solution is numerically intractable especially when the influence of anisotropic properties of metallic layers on dislocations is taken into account (Han and Ghoniem, 2005; Ghoniem and Han, 2005). Compared with 3D-DDD, 2D-DDD modeling is computationally more tractable and cheaper. For this, 2D-DDD modeling has been performed to study the strengthening mechanisms and the layer thickness-dependent strength of MMMFs, with special focus on elastic modulus mismatch (e.g. Koehler barrier), single-dislocation bow-out within layers, interface bypassing strength and dislocation pile-ups at IBs (i.e. the Hall–Petch mechanism) (Keralavarma and Benzerga, 2007). In these 2D-DDD simulations, the crystal anisotropic effect on the deformation field of dislocation is not taken into account, although it might affect the microscopic plastic details produced by each individual dislocation (Shishvan et al., 2011; Ghoniem and Han, 2005). In addition, available 2D-DDD modelings usually employ the superposition framework (Van der Giessen and Needleman, 1995), where the computational costs of polarization stress and long-range interactions between N dislocations are very expensive. In order to lighten the computational burden, unrealistically high loading rate (such as ε̇ = 1000/s ) usually is applied. To overcome these limitations, a feasible approach is to replace the classical 2D-DDD superposition framework by a new 2D-DDD/FEM coupling scheme. In all theoretical models and computational schemes mentioned above, the constraint of IBs to lattice dislocations, which is the most important factor inducing the layer thickness strengthening in MMMFs, is paid high attention to. Due to the strong constraint of IBs to gliding lattice dislocations, dislocations pile-ups can be formed at IBs especially in the MMMFs with layer thickness at the submicron scale or above. With temperature increasing to 0.3Tm and above, where Tm is the melting temperature of the metals, the diffusion of vacancies becomes intensive, and thus the dislocation climb induced by interaction between dislocations and vacancies becomes increasingly important (Ayas et al., 2012, 2014; Ahmed and Hartmaier, 2013; Huang et al. 2014a,b; Hafez Haghighat et al. 2013; Yang et al., 2015). With the aid of vacancies diffusioninduced climb, dislocations can move out of their original slip planes and thereby prevent the build-up of intensive pile-ups at IBs, which may have a marked influence on the strengthening mechanism and the layer thickness effect (Lee et al., 2011; Greer and De Hosson, 2011; Ayas et al., 2012, 2014; Ahmed and Hartmaier, 2013; Huang et al., 2014a,b). To the best of our knowledge, for the MMMFs usually working at elevated temperature due to their high-temperature morphological stability (Mara et al., 2004), how dislocation climb influences the strengthening mechanisms and their layer thickness effect is still open. A thorough study on this topic is of academic and practical importance for the extended use of MMMFs in hightemperature applications. Motivated by above knowledge, a 2D-DDD/FEM coupling scheme based on the hybrid discrete-continuous model (DCM) (Lemarchand et al., 2001; Liu et al., 2009; Vattré et al., 2014) is developed here, with dislocation climb and anisotropic properties of metallic layers being specially included. Different from the often-used 2D-DDD superposition scheme (Van der Giessen and Needleman, 1995; Keralavarma and Benzerga, 2007), no calculation of polarization stress is needed and the computational burden on the long-range interaction between dislocations is also significantly reduced, which make it suitable for modeling of dislocation climb in the MMMFs with difference layer thicknesses in the range from submicron to micron scale. The two major concerns of this study are: to develop a 2D-DDD/FEM coupling scheme suitable for modeling of the strengthening behavior of MMMFs and to reveal the influences of dislocation climb on the strengthening mechanisms and the layer thickness effect of the MMMFs working at elevated temperature. This paper is organized as follows: Section 2 L/2 B

A Cu

Al

Al

Cu

y x slip planes

1200

h1

C

L

h2

D L/2

Fig. 1. Two-dimensional Al/Cu multilayer model.

M. Huang, Z. Li / J. Mech. Phys. Solids 85 (2015) 74–97

77

introduces the climb-enabled and anisotropy-included 2D-DDD/FEM coupling scheme and the computational model of MMMFs. Section 3 presents the results and discussions. Finally, Section 4 ends the paper with some main conclusions.

2. Computational model and 2D-DDD/FEM coupling scheme 2.1. Computational model for MMMFs To capture the layer thickness effect on the strength of MMMFs at elevated temperature and its underlying dislocation mechanism, a 2D representative cell (RC) model, as shown in Fig. 1, was built to perform the DDD modeling under monotonic tension. For simplicity and without loss of generality, the MMMF comprising submicron/micron-sized copper and aluminum layers is considered here. Different from the free-standing single bilayer model, the present RC contains two stacks of Cu/Al bilayers and the periodic displacement boundary condition is applied in the thickness direction (i.e. the xdirection) as follows:

UxA − UxB = U xC − UxD

(1a)

As indicated in Fig. 1, a strain-controlled monotonic tension loading is applied to the RC in the y -direction with given strain rates ε̇ as follows:

UyA − U Cy = UyB − UyD = {UxA,

UxB,

U xC ,

UxD }

1 2

∫ ϵ̇L dt

{UyA,

UyB,

U yC ,

(1b) UyD }

where and are the displacements of representative points A, B, C and D in the x and y directions, respectively. In addition, the thicknesses h1 of Al layer and h2 of Cu layer are selected larger than 100 nm, above which the dislocation pile-up becomes dominant strengthening mechanism. The modulation period h1 + h2 is taken to be one-third of the total length L of the RC. In the present Cu/Al MMMFs, two potential slip systems with an intersection angle of 760° to the x -axis were considered in each layer, as illustrated in Fig. 1. This is an approximation of the three active slip systems for the plane strain state of FCC crystals with relative intersection angle {0°, ± 54.736°} (Rice, 1987). The 0° slip system is neglected since its Schmid factor is 0 (Keralavarma and Benzerga 2007). According to Shishvan et al. (2011), the x -axis and the y -axis are along the [101] and [010] directions, respectively. In addition, the Cu/Al inter-phase boundaries (IBs) are assumed to be fully coherent and no slip system mismatch between two adjacent layers is introduced as shown schematically in Fig. 1, and therefore no initial stress field needs be considered in the present simulation. For simplicity, all IBs are assumed to be dislocation-impenetrable (i.e. rigid). Both Cu and Al layers are assumed to be single crystal without any grain boundary (GB). The elastic behavior of single crystal is usually linearly anisotropic and can be depicted by Hooke’s law σij = Cijkl εij with Cijkl the components of the stiffness tensor. For FCC Cu and Al, there are only three independent elastic constants, namely c11, c12 and c44 . If the basic lattice vectors of the FCC crystal are adopted as the reference coordinate axes, the crystal anisotropy can be reflected by three nonzero elements of stiffness tensor, i.e. Ciiii = c11, Ciijj = c12 and Cijij = c44 with no summation over the repeated i , j indices. Based on the transformation rule given in Shishvan et al. (2011), the elastic tensor in any other coordinate systems (included the selected one shown in Fig. 1) can also be obtained straightforwardly. These independent elastic constants for Cu Cu Cu Cu Al Al Al {c11 , c12 , c44 } and for Al {c11 , c12 , c44 }, respectively, are listed in Table 1 (Shishvan et al., 2011; Hu et al., 2010). With these elastic constants, not only the elastic anisotropic mechanical behaviors of single crystals Al and Cu but also their influence on the dislocation stress field can be considered to a certain extent as discussed in below sections. It should be pointed out that, although these elastic constants are employed in the present dislocation simulations at elevated temperature T¼800 K, they correspond to values at room temperature actually. This approximate treatment for the elastic constants at T¼800 K is also adopted in the DDD simulations performed by Davoudi et al. (2012), Keralavarma et al. (2012) and Ayas et al. (2014). The main emphasis of this study is on how the dislocation climb influences the strengthening mechanism and the layer thickness effect, and thus this approximation should not change the main results and tendencies significantly since the elastic constants only influence the dislocation stress field to a certain extent. 2.2. 2D-DDD/FEM coupling scheme In the last two decades, tremendous progress has been made on the DDD method, which provided an opportunity to Table 1 Elastic constants for FCC Cu and Al layers.

Cu Al

c11 (GPa)

c12 (GPa)

c44 (GPa)

168.2 108.2

121.4 61.3

75.4 28.5

78

M. Huang, Z. Li / J. Mech. Phys. Solids 85 (2015) 74–97

capture the plastic deformation behavior of materials at the submicron-to-micron scale. The DDD frameworks available in the literatures can be roughly divided into two categories, i.e. the superposition scheme and the coupling scheme. In order to solve the boundary condition problem for finite solids, Van der Giessen and Needleman (1995) innovatively proposed a 2DDDD framework by superposing the discrete dislocation dynamics (DDD) and finite element method (FEM), where the DDD result represents the discrete solution for all dislocations in an infinite unbounded crystal and the FEM result is the complementary continuous elastic field which ensures the external and internal boundary conditions of finite solid to be satisfied accurately. By incorporating some additional constitutive rules for the 3D short-range dislocation interactions into this superposition framework, Benzerga et al. (2004) developed a so called 2.5D DDD model. In addition, Shishvan et al. (2011) improved this framework by using the displacement and stress field solutions of dislocations in anisotropic crystal to consider the anisotropic effect. Just recently, Davoudi et al. (2012), Keralavarma et al. (2012), Danas and Deshpande (2013), Ayas et al. (2014) and Huang et al. (2014a, b) further introduced the dislocation climb into this only glide-enabled superposition framework. Although the superposition framework has been extensively employed to model the plastic deformation behavior of materials at the submicron-to-micron scale and achieved great success, it still has some limitations. In the superposition scheme, the long-range interactions between N dislocations are calculated directly, which is an O (N2) problem and thereby the computational burden is very high, especially when the influence of elastic anisotropy on the dislocation field is taken into account (Shishvan et al., 2011). This greatly limits the spatial and temporal scales the superposition scheme can cope with, and thus an ultrahigh loading rate usually is applied in order to save the computational time. In addition, for the present MMMFs, an extra calculation for the so-called polarization stress induced by the difference in elastic properties between adjacent layers must be performed to take the dislocation image stress into account as suggested by Keralavarma and Benzerga (2007). Since dislocations may appear in each metallic layer, the polarization stress must be calculated in every layer of the MMMFs, which also heavily increases the computational burden and greatly reduces the computation efficiency of the superposition scheme. In order to overcome the limitations of the superposition scheme, Lemarchand et al. (2001) proposed a 3D discrete-continuous method (DCM), where the plastic strain ε p in the continuous FEM is calculated by means of the discrete DDD model rather than the classical continuous plastic constitutive relation. In this DDD–FEM coupling scheme, the long-range interaction between dislocations is calculated by the FEM rather than by direct summation of the stress fields induced by all individual dislocation segments as done in the superposition scheme, which evades the computationally expensive O (N2) problem successfully. Liu et al. (2009) further extended this 3D-DCM to solve the finite boundary condition problem, where an initial stress field was introduced to consider the pre-existing stationary dislocations and a Burgers vector density distribution function was adopted to distribute the plastic strain ε p obtained by DDD to the Gaussian integration points of FEM. However, the price to pay the above advantages of DCM is that the singularity of dislocation stress field is smeared out and homogenized in the elements around the dislocation core to a certain extent. If two dislocations are close enough or reside within a mesh element, their elastic interaction is seriously underestimated by the DCM due to the stress homogenization and the continuous interpolation in the FEM. For this, the simulated volume must be meshed as fine as possible, which drastically limits the dimension of computational models. More disadvantageously, for the present MMMFs with abundant IBs where the dislocation pile-ups exist in great numbers, the homogenization of singular dislocation stress field in this DCM inevitably renders the stress induced by the dislocation pile-ups against the IBs to be seriously underestimated, which influences not only the dislocation climb but also the interactions between dislocations in the pile-ups and the interactions between dislocations and IBs significantly. With these in mind, the pervious DCM is improved here by explicitly considering the dislocation–dislocation and dislocation–IB interactions within a small cut-off distance following Lemarchand et al.’s (2001) advice. Similar to the previous DCM, the improved DCM still includes the DDD and FEM modules. In the DDD module, the dynamic evolution of discrete dislocations is

FEM

DD Al

Cu

Al

Cu

Al

Cu

Al

p

Fig. 2. The coupling scheme between the DDD module and the FEM module.

Cu

M. Huang, Z. Li / J. Mech. Phys. Solids 85 (2015) 74–97

79

modeled and the plastic strain ε p induced by the glide and climb of dislocations is directly computed. In the FEM module, the boundary value problem is solved, where the plastic strain ε p is informed directly by the DDD module rather than calculated by means of the phenomenal continuous plastic constitutive laws as done by traditional elastoplastic FEM. In other words, the traditional continuum plastic constitutive calculation is no longer needed in this DDD–FEM coupling scheme. The coupling between these two modules is schematically illustrated in Fig. 2, which involves two key steps. First, the plastic strain ε p , which is computationally obtained by the DDD module, is treated as the eigen-strain and transferred to the FEM module (i.e. DDD→ FEM). Second, the stress σ at the Gaussian integration points (GIPs) of the FE mesh, which is computationally obtained by the FEM module, is interpolated to the discrete dislocations or dislocation sources in the DDD module where the P-K force should be computed (i.e. FEM→ DDD). 2.2.1. Climb-enabled DDD and transfer scheme for plastic strain from DDD to FEM Since a 2D plane-strain small deformation problem is modeled for the MMMFs described in Section 2.1, only edge dislocations are considered here. At the temperature above 0.3Tm , both the climb and glide of these edge dislocations govern the plastic deformation of the MMMFs. As well known, the glide of edge dislocation is confined on its slip plane but the climb is along the normal direction of its slip plane. The in-plane (i.e. glide) and out-of-plane (i.e. climb) motions of edge dislocations are controlled by the glide and climb forces (i.e. the Peach–Koehler force), respectively, which can be expressed as

f g(k) = m(k)⋅Σ(k)⋅b (k) and f c(k) = − s(k)⋅Σ(k)⋅b (k)

(2)

where b(k ) is the Burgers vector of the kth dislocation, and {m(k ), s(k ) } are the unit vectors normal and tangent to its slip plane, respectively. In addition, Σ(k ) is the stress tensor at the position of the kth dislocation, which includes the contributions from the long-range stress fields of all the other dislocations and externally applied load. This tensor can be achieved from the FEM module, of which the computational formulation is provided in the next section. The glide velocity vg(k ) of the kth dislocation is related to the glide force f g(k ) by a simple drag law as follows (Van der Giessen and Needleman, 1995):

vg(k) = f g(k) /Bg

(3)

where Bg is the linear drag coefficient. Different from the dislocation glide, the dislocation climb stems from the interaction between dislocation and vacancies, and thereby it is controlled by the diffusion of vacancies. Therefore, the relation between the climb velocity vc(k ) and the climb force f c(k ) is more complicated. Keralavarma et al. (2012) solved the problems of dislocation climb and vacancy production/diffusion concurrently. Ayas et al. (2012) and Danas and Deshpande (2013) incorporated the temperature, the vacancy volume, the equilibrium vacancy concentration, and average dislocation spacing into the climb drag coefficient Bc and adopted a simple linear drag-type relation similar to Eq. (3). In addition, Davoudi et al. (2012) related the dislocation climb velocity vc(k ) to the vacancy self-diffusion energy ΔEsd and the climb force f c(k ) as follows:

vc(k) =

⎛ f (k) ΔV */b ⎞ ⎤ ⎛ ΔE ⎞ ⎡ 2π D 0 ⎟⎟−1⎥ exp ⎜ − sd ⎟ ⎢ exp ⎜⎜ c ⎝ kB T ⎠ ⎢⎣ b ln (R/b) kB T ⎝ ⎠ ⎥⎦

(4)

where b is the magnitude of the Burgers vector, ΔV * the vacancy formation volume, kB the Boltzmann constant, T the absolute temperature, D0 the pre-exponential diffusion constant, and the argument R/b of the logarithm term can be approximately taken as 2π (Hirth and Lothe, 1968; Mordehai et al., 2008; Davoudi et al., 2012; Huang et al., 2014a, b). In fact, as justified by Ayas et al. (2012), the linear climb drag relation (Danas and Deshpande, 2013), which is a simplification of Eq. (4) with the reasonable assumption fc ΔV * ≪ bkT , is applicable to both the sink and diffusion limited cases. Therefore, although Eq. (4) does not consider directly the coupling of dislocation climb and vacancy diffusion as that by Keralavarma et al. (2012), it is still employed in the present simulations for simplicity. A detailed consideration of the interaction between vacancy diffusion and dislocation climb will be made in future work. In the 2D-DDD modeling by van der Giessen and Needleman (1995), each dislocation source is assumed as a point, and all dislocation sources in the material with a given density are randomly distributed on the slip planes. Once the magnitude of the shear stress at the source exceeds the source strength τnuc during a period of time tns = 10 ns, a dislocation dipole nucleates from it. The source strength τnuc is assumed to follow a Gaussian distribution with the a mean strength τ¯nuc and a standard deviation δτnuc . To consider the effect of Frank–Read source length on the strength of bi-layer metallic thin films, Keralavarma and Benzerga (2007) expressed the dislocation source strength as

τnuc = μb/S

(5)

in which μ is the shear modulus, b the magnitude of the Burgers vector and S the length of the dislocation source. It deserves to be specially noted that, for the considered MMMFs, the length S of Frank–Read sources may has the same order of magnitude as that of the layer thickness, so the improved nucleation model by Keralavarma and Benzerga (2007) is also employed here. The other constitutive law for dislocation annihilation is similar to those for glide-only 2D-DDD simulation, as detailed in Van der Giessen and Needleman (1995). Compared with the dislocation glide, the dislocation climb is very slow. Therefore, the time step Δtclimb for dislocation climb modeling must be chosen carefully, which is 100 times larger than

80

M. Huang, Z. Li / J. Mech. Phys. Solids 85 (2015) 74–97

Fig. 3. The scheme for distributing the plastic strain induced by dislocation motion to the Gaussian integration points in the FEM module.

that for glide. For more details of the climb-enabled 2D-DDD, please refer to Huang et al. (2014a, b). In the DDD modeling, the plastic strain ε ip induced by the glide and climb of dislocation i can be calculated as follows:

εip =

Ai ( ni ⊗ bi + bi ⊗ ni ) 2V

(6)

where Ai is the area swept by dislocation glide or climb, ni the unit vector normal to the slip plane, bi the Burgers vector of the ith dislocation and V is the reference volume. The computationally obtained plastic strain ε ip by the DDD module should be transferred to the Gaussian integration points of the FEM module as the eigen-strains. In the work by Zbib and Diaz de la Rubia (2002), the plastic strain ε ip induced by the ith dislocation was homogeneously transferred to a single element where the ith dislocation resides. Obviously, this simple transfer method can induce a strong discontinuity of plastic strain between adjacent elements and thus influences the solution accuracy of the whole boundary condition problem (Liu et al., 2009). For this, Liu et al. (2009) and Vattré et al. (2014) suggested independently two different methods, by which the plastic strain ε ip is non-uniformly transferred to some Gaussian integration points. This can relieve the incorrect disturbance on the FE results induced by the homogeneous distribution of plastic strain ε ip in a single element (Liu et al., 2009). Considering for this, an improved plastic-strain transfer strategy, which is similar to but distinctly different from that of Liu et al. (2009), is adopted here. The improved transfer strategy for the plastic strain ε ip is illustrated in Fig. 3. When the dislocation i nucleated from the source point Sp glides to the point dp , a plastic strain ε ip is produced and should be transferred/distributed to the domain Ωi as shown in Fig. 3. This domain Ωi comprises the elements with the minimum distance from their center to the slip path Sp dp less than a given value a . This parameter a is analogous to the dislocation core spreading radius (Cai et al., 2006), which can vary with different element size, integration scheme, and the precision of considered problems (Liu et al., 2009). It should be specially noted that, for the considered MMMF with abundant impenetrable inter-phase boundaries (IBs) for dislocations, a careful treatment on the domain Ωj is required. As shown in Fig. 3, when a gliding dislocation j comes up against an IB, the domain Ωj cannot cross the IB and thus is strictly confined in the metallic layer where the dislocation j resides. Within the domain Ωi , the plastic strain εkp distributed to the Gaussian integration point k is assumed as follows:

ϵ kp =

p ωk ∫Ωi ϵ dV ω A = k i ( ni ⊗ bi + bi ⊗ ni ) ∑k ωk Vk ω 2Vk

(7)

where Vk is the volume assigned to the Gaussian integration point k . For the considered 2D case in which a unit length is taken in the direction perpendicular to paper, Ai equals to the length lsdp of glide path Sp dp and Vk equals to one quarter of the element area for the four-integration-points element adopted. In addition, ωk is a plastic strain distribution coefficient on the Gaussian integration point k , which is selected to be the same as the Burgers vector spreading function suggested by Cai et al. (2006)

M. Huang, Z. Li / J. Mech. Phys. Solids 85 (2015) 74–97

ωk =

1 (rk2 + a2)3.5

81

(8)

where rk is the distance between the Gaussian integration point k and the slip plane on which the dislocation i glides. ω = ∑k ωk is the sum of distribution coefficients for all Gaussian integration points in the domain Ωi . It deserves to be specially noted that if the low-order element with four nodes is employed, the distribution coefficients for the four Gaussian integration points (GIPs) 1–4 in a given element shown in Fig. 3 should be the same, although {r1, r2, r3, r4 } are different from each other. Otherwise, a wrong dislocation deformation field might be computationally obtained for the high gradient of the distributed plastic strain in one element, especially when the slip plane intersects obliquely with the element. For this, an effective distance reff = (r1 + r2 + r3 + r4 ) /4 is defined and the distribution coefficients for all the GIPs 1–4 in any given ele2 ment are identically set as ω1 − 4 = 1/(reff + a2)3.5 for the low-order element. For high-order element with eight nodes, no such effective distance and special treatment are needed. It should be mentioned that, in all the plastic strain distribution schemes by Zbib and Diaz de la Rubia (2002), Liu et al. (2009) and Vattré et al. (2014), the incremental form of Eq. (7) was adopted in each time step Δt as follows:

⎧ ϵ p (t ) = ϵ p (t −Δt ) + Δϵ p (Δt ) k k k ⎪ ⎪ ⎨ ⎪ Δϵ p (Δt ) = ωk V Δϵ p = ωk ΔAi ( ni ⊗ bi + bi ⊗ ni ) ⎪ ω Vk i ω 2Vk ⎩ k

(9)

This means that the plastic strain ϵkp (t −Δt ) at the time t − Δt should be saved. In addition, for the considered 2D cases, ΔAi t t −Δt = v Δt of the dislocation i during the time step Δt . When the moving distance v Δt equals to the moving distance lsd − lsd i i p p of the dislocation i is very small, the domain ΩiΔt and the number of GIPs within it should also be very small. As a result, the plastic strains distributed to small number of GIPs in this small domain ΩiΔt are seriously localized and thus an artificial nonuniform plastic strain distribution along the slip path Sp dp of the dislocation i may be induced, which inevitably decreases the accuracy of the following FE calculation. With these in mind, the incremental form of plastic strain distribution (i.e. Eq. (9)) is not employed here. Instead, at the time t , the plastic strain ε ip induced by the gliding dislocation i is totally redistributed by Eq. (7) to all GIPs within the domain Ωi along its slip path from the source point Sp to the current point dp as shown in Fig. 3, where ϵkp (t −Δt ) at the time t − Δt does not need be saved. Nevertheless, if the dislocation i climbs to a new slip plane, its plastic strain distribution before climb still should be saved. Once the dislocation i climbs successfully from one slip plane to another, the new location of the dislocation i is treated as a new source point Sp , which can again be used for plastic strain distribution when the dislocation i glides continually on its new slip plane. By means of these treatments, the distribution domain Ωi and the number of GIPs in it increase gradually with continuous motion of the dislocation i , which leads to a smoother and more appropriate plastic strain distribution along its motion path. On the other hand, when the incremental form of plastic strain distribution is employed as Eq. (9), a special treatment should be performed near the inter-phase boundaries as that by Vattré et al. (2014) to ensure the plastic strain distribution near the phase boundaries is similar with that away from them. However, for the total plastic strain re-distribution scheme employed here, no special treatment is needed when the dislocations approach the phase boundaries since all the GIPs near or away from these boundaries are completely equivalent for the total plastic strain distribution scheme used. As a result, the plastic strain distribution schemes in domain Ωi and Ωj as shown in Fig. 3 are totally the same. It should be specially pointed out that, although the total plastic strain distribution scheme has above merits, it might not suitable to solve the finite deformation problems with large lattice and slip plane rotation. 2.2.2. FEM governing equation and transfer scheme for the stress from FEM to DDD As shown in Fig. 2, after the plastic strain ε p induced by the motion of all dislocations is distributed to appropriate GIPs as the eigen-strain, the standard FEM can be employed to solve the boundary value problem of MMMFs governed by the equation as follows:

[K]{U} = { f a} + { f p}

(10)

where [K] = ∫ [B]T [Ce][B] dV is the elastic stiffness matrix with the independent components of the elastic tensor Ce V listed in Table 1, {U} is the nodal displacement vector, {f a} = ∫ t a [N] dS is the applied force vector, St {f p} = ∫ [B]T [Ce] ε pdV = ∑k Wk [B]T [Ce] εkp is the force vector induced by the plastic eigen-strain ε p transferred from the DDD V module, [B] = grad [N] is the geometrical matrix with N being the shape function vector, Wk is the weight function for Gaussian integral and εkp is the plastic strain distributed to the GIP k . Obviously, Eq. (10) can be readily solved in a unified continuum mechanics framework by the linear finite element method (LFEM). N In the DDD superposition framework, the stress Σp at the point p is calculated as Σp = ^ σFE + ∑i = 1 σ˜i (x p), where ^ σFE is N obtained by employing FEM to solve the image problem in which the fields are smooth, and ∑i = 1 σ˜i (x p) is the analytical solution for the stress at x p induced by all N dislocations. According to this, for a dislocation at x p , the computational cost on its interaction with other dislocations is O (N ). Hence, for all N dislocations at various x p , the total calculation cost is O (N2) for this DDD superposition scheme. In the classical DCM framework, the stress σGIP at the GIPs is calculated as

82

M. Huang, Z. Li / J. Mech. Phys. Solids 85 (2015) 74–97

σGIP = C e : (ϵ − ϵ p),

(11) p

σGIP at

σ NOD

where the plastic strain ε at the GIPs can be achieved by Eq. (7). Once the stress the GIPs and thus at the FE nodes are obtained, the total stress Σp at the point p in an element can be readily solved by the FE interpolation as NODE = 8 , with NI being the shape function of node I . Obviously, different from the DDD superposition Σp = σ FE NI σ NOD p = ∑I = 1 I scheme, the classical DCM scheme avoids the direct calculation of long-range interactions between N dislocations thus can raise the calculation efficiency greatly. However, the classical DCM scheme also has some limitations. Since the stress Σp at the point p is solved by the FE NODE = 8 , the stress Σp near the dislocation core with high level and ultra-high gradient is interpolation as Σp = σ FE NI σ NOD p = ∑I = 1 I inevitably smeared out or homogenized to a certain degree (Lemarchand et al., 2001; Liu et al., 2009). Consequentially, if the distance between the dislocation i and the considered point p is small enough, especially when they reside in the same element, the stress Σp induced by the dislocation i is inevitably underestimated seriously by the FE-interpolation stress σ FE p . This serious underestimation is not permitted for the present modeling on the MMMFs with abundant IBs, since it may induce incorrect configurations of dislocation pile-ups against IBs and thus underestimates seriously the dislocation climb velocity according to Eq. (2). Very refined meshes must be employed in order to improve the precision of the classical DCM. This will increase the calculation cost heavily, especially for those problems with larger sizes, which renders the classical DCM to lose its advantage over the DDD superposition scheme. Furthermore, even if very refined FE mesh is adopted, the classical DCM still cannot appropriately capture the interaction between dislocations residing in the same element (Lemarchand et al., 2001). For these reasons, an enhanced DCM was developed in this contribution. This enhanced DCM algorithm not only overcomes the above limitations of the classical DCM but also has higher computational efficiency than the DDD superposition scheme, and thus can be used to model the MMMFs with abundant IBs and dislocation pile-ups successfully. In our enhanced DCM algorithm, all dislocations around the considered point p are divided into two parts as shown in Fig. 4. The first part includes those blue dislocations residing in the circle cp with its center at the point p and the radius of rcp = 2a (a is the Burger vector spreading radius), and the other part includes those green dislocations residing outside the circle cp . Then the stress Σp at the point p can be calculated by

Σ p = σ FE p +

⎛ 2dpi ⎞ ⎟⋅ ⎝ a ⎠

∑ σTE (x i, x p) − ∑ σTE (x i, x cpi) H ⎜ i ∈ cp

i ∈ cp

(12)

When the stress Σp at the point p is accurately obtained by Eq. (12), it can be transferred to the DDD module to calculate the Peach–Koehler (PK) forces on the point p by means of Eq. (2). If there is a dislocation source at the point p, Eq. (12) can be employed directly; however, if there is a dislocation j at this point, Eq. (12) should be modified by substituting (i ∈ cp ) with (i ∈ cp and i ≠ j ). Then the dislocation nucleation, climb and glide can be determinate by this PK force in the DDD module. Different from the classical DCM, the present DCM algorithm introduces two additional terms (i.e. the second and third terms in the right side of Eq. (12). For those blue dislocations i within the circle cp (i.e. i ∈ cp ), the stress at the point x p

Cu

Al

rcp =2a C

X pi

p

Xp

Circle C p

dpi

dislocation i

Circle Ci

Fig. 4. Schematic diagram for calculating the stress field at point p induced by the dislocation i .

M. Huang, Z. Li / J. Mech. Phys. Solids 85 (2015) 74–97

83

induced by them cannot be accurately captured by the first FE interpolation term σ FE p thus some amendments to it must be made as denoted by the second and third terms in the right side of Eq. (12). However, for those green dislocations i outside the circle cp (i.e. i ∉ cp ), the stress at the point x p induced by them can be described by σ FE p accurately enough, and thus two amendment terms in Eq. (12) should automatically disappear. In order to meet these requirements, the analytical stress solution σ TE (x i, x p) is specially introduced to depict directly the stress at the point p induced by those blue dislocations i within the cut-off circle cp (i.e.i ∈ cp ). For the micro-layered metallic thin films (MMTFs) considered, the analytical formula of

σ TE (x i, x p) developed by Lubarda (1997) for the dislocation near a bi-material interface was adopted to take the influence of inter-phase boundaries on the dislocation stress field into account. To avoid repetition of consideration, the contribution to the first FE-interpolation term σ FE p by those blue dislocations i within the cut-off circle cp (i.e.i ∈ cp ), which can be denoted as

σ FE p

i ∈ cp

, must be removed. For this, the third term −σ FE p

i ∈ cp

= − ∑i ∈ cp σTE (x i, x cpi ) H (2dpi /a) is introduced into Eq. (12) for those

blue dislocations i ∈ cp , where σ TE (x i, x cpi ) is the analytical stress at the point x cpi (x cpi = x i + rci (x p − x i/ x p − x i )) induced by the dislocation i at x i , as shown in Fig. 4. Obviously, the point x cpi is the radial projection of the point p on the blue circle ci with its center at point x i and radius rci = 2a . Since the point x cpi is away from the dislocation i enough, the FE-interpolation stress

σ FE p

i ∈ cp

at x cpi is basically equal to the analytical solution ∑i ∈ cp σTE (x i, x cpi ) (Lubarda, 1997) at x cpi . However, at the point p be-

tween x cpi and x i , σ FE p

i ∈ cp

varies with x p changing along the path from x i to x cpi , it cannot be accurately expressed only by the

analytical solution σ TE (x i, x cpi ) at x cpi . For this, a function H (Ria ), which can capture accurately enough the variation of σ FE p

i ∈ cp

with x p , is specially introduced into the third term as

⎛ 2dpi ⎞ 25 Ria H ( Ria ) = H ⎜ , ⎟= 4 (Ria + 1)2 ⎝ a ⎠

(13)

where dpi = x p − x i is the distance between the point p and dislocation i , and the parameter Ria is defined as 2dpi /a . When x p = x cpi , Ria = 2dpi /a = 4 and H (Ria ) = 1. Why the function H (Ria ) in the third term is selected as Eq. (13) is explained in Appendix A. With the introduction of last two terms in Eq. (12), the present algorithm becomes not as sensitive to the FE mesh as the classical DCM scheme. Even for coarser mesh, the stress field near the dislocation core can be accurately calculated by the second theoretical terms in Eq. (12), and the stress field outside the circle with a cut-off radius rci = 2a also can be appropriately achieved only by the first FE-interpolation term. Our numerical results show that the computationally obtained stress–strain curves of MMTFs by the present DCM algorithm are basically the same for two different mesh sizes, which is not given out here owing to the paper length limitation. It should be specially mentioned that, since the cut-off radius is suggested as rci = 4l as shown in the Appendix A, the cut-off radius and the dislocation number M in the cut-off circle inevitably increase when coarser meshes are employed. Although coarser mesh can decrease the calculation cost in FEM, it increases the computational cost on the interaction between dislocations in the cut-off circle, thus a proper mesh size should be chosen. If the mesh size in the present DCM is the same as that used in the classcial DDD superposition scheme, their computational precision should be very close. On the other hand, for each point p, the computation cost of present DCM algorithm is O (M ) scale, with M being the number of dislocations in the circle cp . If there are N dislocations in the whole system, the stress should be calculated on these N points, and thus the total computation cost is O (N⋅M ) scale. It deserves to be noted that the calculation cost of O (N⋅M ) basically increases linearly with increasing the total dislocation number N , and thus the computational efficiency is almost similar to the fast multi-pole expansion scheme by Keralavarma (2011). In addition, the so-called dislocation polarization stress induced by elastic property mismatch between two phases (Van der Giessen and Needleman, 1995; Keralavarma and Benzerga, 2007) needs not to be calculated directly in the present DCM scheme, which also can reduce the computational burden substantially. Since M in the small cut-off circle cp is much less than the total number N , the present DCM algorithm is more efficient than the DDD superposition scheme with the total computation cost being O (N2) scale, especially for general system with larger size and more dislocations N . Obviously, higher computational efficiency of the present DCM algorithm is attributed to the cut-off circle cp introduced. For the dislocation located at the center of the circle cp , only the interactions between it and those M dislocations within the circle cp need be considered by the second and third terms of Eq. (12); however, the interactions between it and those dislocations outside the circle cp can be accurately obtained only by the first FE-interpolation term of Eq. (2) as done in the classical DCM scheme. Hence, although the computational cost of present enhanced DCM slightly increase compared with the classical DCM method due to the introduction of last two terms in Eq. (12), it is still much smaller than that the DDD superposition scheme. To sum up, the present enhanced DCM has several merits. Firstly, it has similar high computational efficiency to the classical DCM and high computational precision to the DDD superposition scheme. Since this enhanced DCM algorithm can dramatically reduce the computation burdens, it can be used to perform the simulations with loading rate (ϵ = 50/s or even smaller) significantly lower than that in the classical DDD methods. Most importantly, the limitation of classical DCM is

84

M. Huang, Z. Li / J. Mech. Phys. Solids 85 (2015) 74–97

effectively overcome by the present enhanced DCM algorithm. The significant underestimation of the ultra-high stress near the dislocation core (where the analytical stress should be singular but usually is smeared out to zero by FE interpolation) is successfully avoided. The interaction between neighboring dislocations can be considered accurately as done in the DDD superposition framework even they resides in the same element. Hence, there is no special requirement for FE mesh refinement in the present DCM algorithm, even when the MMTFs with abundant IBs are simulated. Secondly, since the plastic strain ε p is transferred to the FEM module where the anisotropic elastic constants are assigned to different metallic layers, the influences of anisotropic elasticity and crystal orientation on the stress field of dislocations can be depicted by the first term σpFE to a certain extent, although no explicit complex anisotropic dislocation stress formulation is introduced as that by Shishvan et al. (2011). In addition, the image stress on the dislocation induced by IBs can also be captured by the first term σpFE or the second term in Eq. (12) directly. Finally, by solving Eq. (10) with FEM, the displacement fields both by the discrete dislocations and by the external loadings can be readily achieved. No analytical dislocation displacement field is needed in the present enhanced DCM algorithm. Provided the plastic strain ε p by the dislocation glide and climb is distributed along its moving path correctly, no special correction on dislocation displacement field should be appended as done by Ayas et al. (2012) and Huang et al. (2014) to ensure the right dislocation displacement discontinuity when the dislocation climbs out of its slip plane. This makes the present algorithm is easy to implement. For above these merits, the present DDD–FEM coupling scheme is appropriate to model the micro-scale plastic behavior of the submicron-to-micron sized materials with abundant interfaces (such as grain boundaries and phase boundaries) or passivated/coated surfaces. The detailed validation of the present amended DCM is provided in Appendix A for both singlephase crystals and materials with internal interfaces.

3. 2D DDD/FEM coupling modeling on mechanical behaviors of MMMFs The microlayered Cu/Al multilayer thin film with the layer thickness h ≥ 250 nm is considered, as shown in Fig. 1, with a special focus on how the dislocation climb influences the strengthening mechanism and the layer thickness effect at elevated temperature. For the FCC Al crystal, the density of dislocation sources ρsource generally falls into the range of 30–300 μm  2 (Keralavarma and Benzerga, 2007) and ⍴source ¼30–60 μm  2 is the typical value at room temperature. At the high temperature T = 800 K considered here, the dislocation source densities in both Al and Cu layers are assumed to 120 μm−2. Similar to Keralavarma and Benzerga (2007), the drag coefficients in Eq. (3) are selected as B = 10−4Pa s for both two materials. Uniform FE meshes with element size of l = 0.025 μm are used in all modeling of MMMFs with different layer thicknesses. As validated in Appendix A, this sized mesh can capture the dislocation-induced stress fields accurately enough. Eight-nodes and four-Gaussian-points plane-strain iso-parametric square elements are employed. The elastic properties for FE calculations are provided in Table 1, the elastic constants associated with the dislocation stress field and dislocation source strength are μAl = 26 GPa and vAl = 0.36 for Al and μCu = 41 GPa and vCu = 0.33 for Cu, respectively. The material parameters in Eq. (5) for dislocation climb are listed in Table 2, which are representative of Al (Davoudi et al., 2012) and Cu (Butrymowicz et al., 1973), respectively. Compared with the existing climb-enabled DD simulations on materials with multiple internal slip-blocking boundaries, such as passivated surfaces (Ayas et al., 2012) and grain boundaries of thin films (Davoudi et al., 2012), the present DDD calculations on the MMMFs have the following differences. First, unless otherwise specified, the strain rate applied at both ends of the MMMFs is fixed at ϵ = 50/s, which is much lower than 1000/s by Ayas et al. (2012, 2014) and 4000/s by Davoudi et al. (2012). Second, the influence of modulus mismatch on the dislocation climb effect is first investigated. Third, different from the works by Ayas et al. (2012) and Davoudi et al. (2012) where the influence of layer thickness on the dislocation source strength was neglected, the present paper considers it by the length of dislocation sources S in Eq. (4), which is randomly generated in the range of Cl and Ch . In the work by Keralavarma and Benzerga (2007), Cl is fixed at 32b and Ch is selected as one-third of the layer thickness h. Considering the layer thickness in this study is much larger than that in Keralavarma and Benzerga (2007), we suggest that Ch is expressed as follows:

Ch = min (h/3, βμb/τ¯nuc )

(14)

where τ¯nuc is the mean strength of dislocation sources in bulk materials and the factor β = 1 scales the source length. In the Cu Al is selected as 50 MPa (Cleveringa et al., 2000) for Al and τ¯nuc present simulations, τ¯nuc = 50 MPa × μCu /μAl for Cu. With this consideration, the influence of layer thickness on the dislocation source strength and thus on the size-dependent responses of MMMFs can be investigated. Table 2 Climb parameters for Cu and Al layers.

Al Cu

Burgers vector, b (nm)

ΔEsd (eV)

D0 (cm2/s)

0.25 0.2557

1.28 2.036

0.1185 0.18

M. Huang, Z. Li / J. Mech. Phys. Solids 85 (2015) 74–97

85

Dashed: without climb considered 706

Red: 0.25

m

Blue: 0.5

m

Green:1.0

m

Stress

(MPa)

solid: with climb considered

T=800K 483

305 244 203 187

pn

pc

Strain

p

Fig. 5. The computationally obtained stress–strain curves of MMMFs with different layer thicknesses h = {0.25, 0.5, 1.0} μm , without and with dislocation climb considered for comparison.

3.1. Influence of dislocation climb on the stress–strain responses of MMMFs To investigate the effect of dislocation climb on the stress–strain behavior of MMMFs, the computationally obtained σ ~ε curves with and without dislocation climb considered are plotted together in Fig. 5 for three different layer thicknesses h = {0.25, 0.5, 1.0} μm . Similar distributions of dislocation sources are assigned for two cases with and without climb considered. It is clearly shown in Fig. 5 that, when the dislocation climb is taken into account, both the flow stress σ and hardening rate ∂σ /∂ε are significantly reduced for all the layer thicknesses h considered, especially when the layer thickness h is lower. For comparison, the stress at the plastic strain εp = 0.2% is defined as the nominal yield stress σ0.2 (i.e. the stress value of the intersection point between the stress–strain curves and the εp = 0.2% offset line). It can be found from Fig. 5 that, when the dislocation climb is considered, σ0.2 is reduced by 223 MPa from 706 MPa to 483 MPa for the MMMF with smaller layer thickness h = 0.25 μm but only by 16 MPa from 203 MPa to 187 MPa for the MMMF with larger layer thickness h = 1.0 μm . At least two conclusions can be drawn from this figure that: (1) with consideration of dislocation climb, the Hall–Petch-like layer thickness effect on the strength seems still valid for the considered MMMFs with layer thickness larger

T=800K

h=0.5

m

Red: Without climb considered

Dislocation Number (N)

Blue: With climb considered No symbol: Total dislocation number Square: Dislocation number in Al layers Circle: Dislocation number in Cu layers

Strain Fig. 6. The dislocation number in the whole MMMF and in individual layers with the layer thickness h = 0.5 μm , without and with dislocation climb considered for comparison.

M. Huang, Z. Li / J. Mech. Phys. Solids 85 (2015) 74–97

m -2

86

With climb considered

Dislocation density

h=0.25

T=800K

m

h=0.375 m h=0.5

m

h=1.0

m

Strain Fig. 7. The dislocation density in the MMMFs with different layer thicknesses h = {0.25, 0.375, 0.5, 1.0} μm , with the consideration of dislocation climb.

than 250 nm , i.e. the thicker the layer thickness h is, the lower the nominal yield stress σ0.2 is; (2) the dislocation climb can significantly influence the Hall–Petch-like layer thickness effect, which will be further analyzed in the following section. As we all known, the dynamic evolution of dislocations in the MMMFs is responsible for their stress–strain responses. Considering for this, the evolution of total dislocation number Ntot in the MMMFs with the applied strain ε is plotted in Fig. 6 for the layer thickness h = 0.5 μm , where Ntot = NAl + NCu with NAl and NCu being dislocation numbers in Al layers and Cu layers, respectively. For comparison, the curves with and without consideration of dislocation climb are plotted together. It can be clearly seen that, at a given applied strain ε , the dislocation numbers Ntot , NAl and NCu with consideration of dislocation climb are markedly lower than those without consideration of it. However, at a given strain (e.g. ε = 0.011) shown in Fig. 5, the corresponding plastic strain εpc with dislocation climb considered is much larger than εpn without it. This means that, when the dislocation climb becomes active at elevated temperature, less dislocation can induce larger plastic strain. This is because, when the climb is activated, the dislocations in the pile-ups against Al/Cu interphase boundaries can climb away from their original slip plane to a new slip plane and glide forward, and thereby they can produce larger plastic strain εp . Further, according to the simple relation σ = E (ε − εp ), the larger plastic strain εp by dislocation climb can reduce the flow stress level σ and hardening rate as shown in Fig. 5; as a result, the number of dislocations successfully nucleated from their sources is reduced under this decreased stress level. In addition, it is also shown in Fig. 6 that, when the dislocation climb is neglected, the dislocation number in the Al layers is basically the same with that in the Cu layers. However, with the dislocation climb taken into account, the dislocation number in the Al layers becomes much less than that in the Cu layers.

a

b

Dashed : without climb considered

Dashed: without climb considered Solid: with climb considered 1011MPa

(MPa)

Solid: with climb considered Red : 4000/s

Blue: Ni/Al system Red: Cu/Al system

(MPa)

Blue: 1000/s Green: 50/s

Stress

Stress

670MPa 626MPa

439MPa

h=0.5 m h=0.25 m

p

Strain

p

Strain

Fig. 8. The influence of (a) loading rate and (b) modulus mismatch on the dislocation-climb effect.

M. Huang, Z. Li / J. Mech. Phys. Solids 85 (2015) 74–97

87

This indicates that the dislocation climb can result in inhomogeneous dislocation distribution across different layers, which can aggravate heterogeneous plastic deformation in the MMMFs. To further understand the layer thickness effect in the MMMFs as shown in Fig. 5, the total dislocation density ρdis in the MMMFs is plotted in Fig. 7 as a function of the applied total strain ε for four layer thicknesses h = {0.25, 0.375, 0.5, 1.0} μm , with the dislocation climb considered. It is clearly shown that, with the layer thickness h decreasing from 1 μm to 0.375 μm , ρdis gradually increases although this increasing is not very remarkable. However, with the layer thickness h further decreasing to 0.25 μm, ρdis does not increase instead decrease to a certain extent. A similar result also can be obtained when no dislocation climb is taken into account. However, as shown in Fig. 5, the yield stress and hardening rate of MMMFs monotonically increases with the decrease of layer thickness h, so the strengthening mechanisms of the MMMFs and the layer thickness effect might be dominated by the dislocation pile-up against IBs and the confined dislocation–source length by limited layer thickness, rather than by the mechanism associated with the dislocation density ρdis , especially when the layer thickness is smaller. 3.2. Influence of loading rate and modulus mismatch on the dislocation climb effect Due to the merits of present discrete-continuous model, a strain rate ϵ = 50/s much lower than those in the previous climb-enabled DD simulations (ϵ = 1000/s and 4000/s, Ayas et al., 2012, 2014; Davoudi et al., 2012) can be applied on the MMMFs with multiple internal phase boundaries. To investigate the influence of loading rate on the dislocation climb effect, the stress–strain responses of MMMFs with layer thickness h = 0.5 μm at three different loading strain rates ϵ = {50, 1000, 4000}/s are plotted in Fig. 8(a) for both cases with and without dislocation climb considered. For comparison, the same initial dislocation-source density and spatial distribution are assigned for all the cases considered. It can be clearly seen from Fig. 8(a) that, when only dislocation glide is considered, very weak strain loading rate effect is captured on the stress–strain curves. This is true for this case because there are only two parameters associated with the time-dependent evolution of dislocations in the only glide-enabled DDD simulations. One is the drag coefficient B which is related to the gliding velocity of dislocations, and the other is the Frank–Read source nucleation time tnuc . It has been validated by our simulations that changing B in the range of {10−3, 10−4 , 10−5} Pa s for general metals will not induce a remarkable loading rate effect for these glide-only cases. On the other hand, as long as the Frank–Read sources have sufficient time to nucleate at all the loading rates considered, the simulated stress–strain responses should be irrespective of ε̇. However, when the dislocation climb is taken into account, the flow stress at any given strain significantly decreases with decreasing the loading rate, showing remarkable loading rate effect. Hence, we can say that different rate sensitivities should mainly originate from the dislocation climb rather than the drag coefficient B. In other word, the drag coefficient B will not lead to different tendency of rate sensitivity. A possible explain is that the dislocation climb is much slower than dislocation glide, so the dislocation can climb more fully at low loading rate than that at high loading rate. As a result, the loading rate effect on the stress–strain response is captured when dislocation climb is considered. This is also consistent with the experimentally observed significant effect of deformation rate on the compressive behavior of indium micropillars, where dislocation climb was shown to come into effect (Lee et al., 2011). Moreover, as shown in Fig. 8(a), the difference in nominal yield stress σ0.2 between with and without the dislocation climb being considered is significant for ϵ = 50/s, but it becomes insignificant for both ϵ = {1000, 4000}/s. In other words, if an over-high loading rate is applied, the dislocation climb effect on the stress–strain response of MMMFs becomes not so remarkable, thus a lower loading rate ϵ = 50/s is employed here to capture the dislocation climb effect effectively. As mentioned in Section 2, the influences of the modulus mismatch between different layers and the dislocation image stress due to interphase boundaries can be considered in the present DCM model, although no special calculation of the socalled polarization stress is performed as Keralavarma and Benzerga (2007). To reveal clearly the influence of modulus mismatch on the dislocation climb effect, a Ni/Al MMMFs with higher relative stiffness ENi/EAl = 3 (νNi = 0.31) is specially considered here together with the Cu/Al MMMFs with ECu/EAl ≈ 1.5. For the Ni at 800 K, the vacancy self-diffusion energy ΔEsd is selected as 2.88 eV and the pre-exponential diffusion constant is set as 0.92 cm2/s (Maier et al., 1976). Fig. 8(b) plots the stress–strain curves of two MMMFs captured by the present DCM method with and without dislocation climb being considered. For simplification and without loss of generality, only the layer thickness of h = 0.25 μm with the same initial dislocation-source density and distribution is considered here. From Fig. 8(b), it can be clearly seen that, with dislocation climb considered, the nominal yield stress σ0.2 only decreases about 187 MPa from 626 MPa to 439 MPa for the Cu/Al system, while it decreases about 341 MPa from 1011 MPa to 670 MPa for the Ni/Al system. Even if the stress decrease Cu / Al Ni / Al climb no climb is normalized by the effective modulus Eeff = (E1 + E2 ) /2 (Eeff Δσ0.2 = σ0.2 − σ0.2 ≈ 90 GPa and Eeff ≈ 140 GPa , Keralavarma and Benzerga, 2007), the value of Δσ0.2/Eeff is much larger for Ni/Al system with ENi/EAl = 3 than that for Cu/Al system with ECu/EAl ≈ 1.5. It seems that the modulus mismatch between layers can enhance the dislocation climb effect. In addition, the simulations also show that, similar to those in the Cu/Al system, almost all the dislocation climbs occur at the Al layers in the Ni/Al system. However, the dislocation climb in the Ni/Al system is more significant than that in the Cu/Al one. This might be explained as follows. In Ni/Al MMMFs, the dislocation pile-ups against the inter-phase boundaries experience larger repulsive image stress due to the higher modulus mismatch, leading to more significant dislocation climb and thus larger Δσ0.2/Eeff . As a total result, the dislocation climb effect is more remarkable in the MMMFs with higher modulus mismatch between layers.

88

M. Huang, Z. Li / J. Mech. Phys. Solids 85 (2015) 74–97

a

Al

Cu

Al

b

Cu

Al

Cu

Cu

Al

T=800K

Y

Y

T=800K

Without

With

climb

climb

considered

considered

X

X

Fig. 9. The contour of equivalent slip quantity Γ in the MMMF with the layer thickness h = 0.5 μm at applied strain ε = 0.012 (a) without and (b) with dislocation climb considered.

3.3. Influence of dislocation climb on the microscopic deformation fields in MMMFs The overall stress–strain behavior of MMMFs and the layer thickness effect on it are closely related to both the dynamic evolution of dislocations and the microscopic deformation fields in individual layers. To study the influence of dislocation climb on the microscopic deformation field in the MMMFs, a parameter Γ named as the microscopic equivalent plastic slip is defined as follows (Balint et al., 2006):

Γ = γ 1 + γ 2 with γ i = mi⋅˜ϵ p⋅ni (with i no sum)

(15)

i

where γ is the slip quantity for the i th slip system, mi and ni the tangential and normal vectors of this slip system, respectively, and ε˜ p the plastic strain tensor induced by the discrete dislocations. The contours of the microscopic plastic slip Γ in the MMMF with layer thickness h = 0.5 μm are plotted in Fig. 9 for both cases with and without dislocation climb considered at the applied strain ε = 0.012. It should be pointed out that, although the periodic displacement boundary conditions are applied as shown in Fig. 1, the lateral outer boundaries where the periodic boundary conditions are exerted are exactly the rigid layer interfaces, and thus the inherent jumps of slip/stress contours between different layer interfaces make the total slip/stress field discontinuous (i.e. not periodic like displacement) even across the periodic boundaries. Comparison between Fig. 9(a) and (b) indicates that, when no dislocation climb is considered as shown in Fig. 9(a), the

a

Al

Cu

Al

b

Cu

Al

Cu

Cu

Al

T=800K T=800K yy yy

Pa

Y

Y

Pa

Without climb

With

considered

climb considered

X

X

Fig. 10. The contour of tensile stress σyy in the MMMF with the layer thickness h ¼ 0.5 μm at the applied strain ε = 0.012 (a) without and (b) with dislocation climb considered for comparison.

M. Huang, Z. Li / J. Mech. Phys. Solids 85 (2015) 74–97

89

plastic slip Γ is lower and its distribution is basically similar in different metallic layers; however, if the dislocation climb is taken into account, as shown in Fig. 9(b), the plastic slip Γ is much larger and seriously localized into some slip bands, and the plastic slip in different metallic layers is non-uniform with averaged Γ in the Al layers much higher than that in the Cu layers. The main reason is that the dislocation climb in the Al layers is more significant than that in the Cu layers. Since those dislocations which climb successfully away from the pile-ups in the Al layers have a longer free glide path and can efficiently relieve the dislocation pile-ups, the average plastic slip in the Al layers are larger than that in the Cu layers as shown in Fig. 9 (b), although the number of nucleated dislocations in the Al layers is much less than that in the Cu layer as shown in Fig. 6. This validates again that, with the assistance of dislocation climb, the mobility of dislocations is greatly increased, and thus less dislocations can produce larger plastic deformation. Besides the plastic slip, the microscopic stress fields in individual layers of MMMFs also can be influenced by the dislocation climb. In order to clarify this influence, the contours of the normal stress σyy in the MMMF with layer thickness h = 0.25 μm are presented in Fig. 10 at the applied strain ε = 0.012 for both cases with and without dislocation climb included. It is clearly shown that for both cases, intense stress concentration can be observed near the Al/Cu IBs parallel to the tensile axis, which is believed to be induced by the dislocation pile-ups against the Al/Cu interfaces. A careful comparison between Fig. 10(a) and (b) indicates that, when the dislocation climb is neglected, the stress level and its concentration in the Al layers are basically similar to those in the Cu layers; however, when the dislocation climb is taken into account, the stress in the whole MMMF are greatly reduced especially in the Al layers, and significant stress concentration can only be observed at the Cu layer side of the Al/Cu interface. Obviously, the dislocation climb leads to more non-uniform stress distribution in different metallic layers of MMMFs. This is because that, at the temperature T ¼800 K considered, the dislocation climb in the Al layers is more significant than that in the Cu layers. On the one hand, the dislocation climb can induce larger plastic strain in the Al layers as shown in Fig. 9(b), which reduces significantly the stress in the Al layers according to the simple relation σ = E (ε − ε p). On the other hand, the dislocation climb renders the number of dislocations nucleated in the Al layers to decrease as shown in Fig. 6, which can further reduce the number of dislocations composing the pile-ups against IBs in the Al-layer side. As a result, the stress concentration induced by dislocation pile-ups at the Al-layer side of the Al/Cu IBs is greatly decreased as shown in Fig. 10(b). To understand the influence of dislocation climb on the microscopic deformation field more deeply, the dislocation slip traces in the MMMF with layer thickness h¼0.5 μm are presented in Fig. 11(a) and (b) at the applied strain ε = 0.012 for both cases with and without dislocation climb considered, respectively. Since the dislocation slip traces are defined as the trajectories of dislocation glide and climb, Fig. 11 is not a dislocation configuration at one given moment but the superposition of dislocation configurations at different loading steps. It is evident from Fig. 11(a) that, when the dislocation climb is neglected, the numbers of active slip bands in different metallic layers are basically the same. This is the reason why the plastic deformation is almost similar in different metallic layers of the MMMF without consideration of dislocation climb, as shown in Fig. 9(a). However, when the dislocation climb is included, the number of activated slip bands decreases markedly, especially in the Al layers where dislocation climb occurs more significantly. Those climbing dislocations can spread toward both sides of their original slip planes and effectively relieve the dislocation pile-ups against the Al/Cu interfaces. Continuous dislocation climb makes the slip bands become wider and wider, thus longer free gliding path is provided for dislocations, resulting in

Fig. 11. The trajectories of dislocation motion in the MMMF with the layer thickness of h ¼0.5 μm at the applied strain ε = 0.012 (a) without and (b) with climb considered for comparison.

90

M. Huang, Z. Li / J. Mech. Phys. Solids 85 (2015) 74–97

Without climb considered

(MPa)

Simulation results:

With climb considered

Without climb considered

Fitting curves:

Yielding stress

With climb considered

-0.85

60+160*h

30+160*h

-0.65

h

( m)

Fig. 12. The nominal yield stress σ 0.2 of the MMMFs vs. the layer thickness h for two cases with and without dislocation climb considered.

higher plastic deformation in these wider slip bands. This is distinctly different from the case without the consideration of dislocation climb, where dislocations only glide on their original slip plane as shown in Fig. 11(a) and their motion is strongly constrained by the pile-ups against the IBs. Due to lower dislocation mobility, only less plastic deformation can be produced as shown in Fig. 9(a), although more dislocations are successfully nucleated on more activated slip planes as indicated in Fig. 6. In addition, even in the case with dislocation climb considered, since dislocation climb in the Al layers is stronger than that in the Cu layers, distinct differences in the width and number of dislocation slip traces exist between in the Al layer and in the Cu layer as shown in Fig. 11(b), further resulting in significant differences in the plastic strain and the stress distribution among them as shown in Figs. 9(b) and 10(b), respectively. In summary, by increasing the dislocation mobility and relieving the pile-ups (especially those in the Al layers) against IBs, the dislocation climb can markedly change the number of successfully nucleated dislocations, the width and number of dislocation slip traces, and thus can render more non-uniformly distributed microscopic plastic strain and stress across the Al and Cu layers, which can further influence the strengthening behavior of MMMFs and its layer thickness effect significantly. 3.4. Influence of dislocation climb on the layer thickness effect of MMMFs As mentioned above, due to the dislocation pile-ups against IBs, the strength of MMMFs strongly depends on the layer thickness h, showing strong layer thickness effect. At elevated temperature, some dislocations can climb away from their pile-ups and achieve a longer free gliding path. Obviously, the dislocation climb as an important relaxation mechanism to dislocation pile-ups against IBs can significantly influence the layer thickness effect of MMMFs. To reveal the influence of dislocation climb on the layer thickness effect, the variation of the nominal yield stress σ0.2 with the layer thickness h is plotted in Fig. 12 for both cases with and without climb considered. For each layer thickness, at least three different random distributions of dislocation sources are considered in the present simulations. As reported by Misra et al. (2005), when the layer thickness h ≥ 0.1 μm , the Hall–Petch model can describe the dependence of the yield stress σ0.2 on the layer thickness h as σ0.2 = σ0 + kh−n , with k being the Hall–Petch slope, n the layer-thickness exponent and σ0 the yield stress of MMMF with sufficiently large h. From Fig. 12, it can be clearly seen that the computationally obtained dependence of yield stress σ0.2 on the thickness layer h can also be well described by the Hall–Petch model for both two cases with and without consideration of dislocation climb, but there are distinct difference in the model parameters {σ0, n, k} listed in Table 3, with R2 being the fitting correlation coefficient. It can be found from Table 3 that the layer-thickness exponent n is higher than 0.5 for both two cases with and without Table 3 The fitting parameters of Hall–Petch relations for MMMFs. Parameters

σ 0 (MPa)

k (MPa)

n

R2

No dislocation climb Dislocation climb

60 30

160 160

0.85 0.65

0.95 0.96

M. Huang, Z. Li / J. Mech. Phys. Solids 85 (2015) 74–97

91

Nc

m

T=800K With climb considered

h=0.25

m

h=0.50

m

h=0.75

m

h=1.00

m

cr

Strain Fig. 13. The number of climb events happened per unit area vs. the applied strain for the MMMFs with different layer thickness.

consideration of dislocation climb, showing stronger layer thickness effect than that predicted by classic Hall–Petch relation. This is because that, the length of dislocation sources in the MMMFs is confined by the layer thickness h, which is considered by Eq. (14) in the present DDD simulations. When the thickness layer h ≤ 0.375 μm , the dislocation-source length in the MMMFs is determined by the one-third of layer thickness (h/3) rather than by the average source length in bulk materials μb/τ¯nuc according to Eq. (14), which can enhance significantly the layer thickness effect since h/3 < (μb /τ¯nuc ), especially for the MMMFs with smaller layer thickness (i.e. h = 0.25 μm ). Obviously, besides the dislocation pile-up mechanism, the confined dislocation-source length mechanism induced by the layer thickness constraint also contributes to the layer thickness effect of MMMFs. A similar mechanism with n = 1 is also induced in the fine-grain polycrystals by the 3D dislocation curvature confined by grain boundaries (Ohashi et al., 2007). On the other hand, the layer-thickness exponent n in the cases with consideration of dislocation climb is smaller than that without it. This indicates that the dislocation climb at elevated temperature (800 K) weakens the layer thickness effect on the yield strength of MMMFs. This weakening effect is believed to come from the relaxation of pile-ups against the IBs induced by the dislocation climb. Moreover, the yield stress σ0 of the MMMF with sufficiently large h is lower for the case with consideration of dislocation climb than that without it, which clearly shows that the dislocation climb can significantly decrease the yield stress even for the MMMFs with sufficiently large layer thickness. In the present DDD algorithm with dislocation climb considered, a dislocation can climb successfully from one slip plane to another if its climb distance is larger than the magnitude of a Burgers vector. By assuming the climb distance is basically the same for every climb event, the number of climb events Nc occurred in unit area can be used to describe the extent (or degree) of dislocation climb to a certain extent. For this reason, the climb events Nc per unit area (i.e. density of climb events) is plotted in Fig. 13 as a function of the applied strain ε for the MMMFs with various layer thicknesses h = {0.25, 0.5, 0.75, 1.0} μm . It is clearly shown from this figure that, when the applied strain ε > εcr ≈ 0.006, Nc increases monotonically with decreasing the layer thickness h. This means that, when the applied strain is large enough, the dislocation climb event takes place more frequently in the MMMFs with smaller layer thickness, which results in more marked weakening to the layer thickness effect on the strength of MMMFs with h ≤ 0.5 μm as shown in Fig. 12.

4. Conclusions An extensive hybrid discrete-continuous model (DCM), which couples the 2D climb/glide-enabled DDD to the linearly elastic FEM, is developed. In this 2D-DDD/FEM coupling scheme, the dynamic evolutions of discrete dislocations are captured and the plastic strain induced by moving dislocations is achieved by the 2D climb/glide-enabled DDD module, while the stress field induced by both the external loading and the discrete dislocations are solved by the FEM module. The coupling scheme between the DDD module and FEM module includes two key steps: one is how the plastic strain obtained by the DDD module is properly transferred to the FEM module as an eigen-strain; the other is how the stress obtained by the FEM module is accurately transferred to the DDD module to drive those discrete dislocations to evolve. In order to solve these two key problems, some especial treatments are made in this study. A careful comparison between the results obtained by the present numerical scheme and the solutions given by the classical dislocation theory indicates that the present

92

M. Huang, Z. Li / J. Mech. Phys. Solids 85 (2015) 74–97

coupling scheme is valid and the algorithm is accurate enough. By using this coupling scheme, the tensile behaviors of microlayered Al/Cu multilayer films are modeled, and the layer thickness effect on the yield strength and its underlying dislocation mechanisms are studied carefully, with a special attention to the influences of dislocation climb on them. Some interesting results are summarized as follows: 1. With dislocation climb considered, the plastic flow stress and hardening rate of MMMFs are significantly decreased; meanwhile, the dislocation density is also reduced significantly. This means that less dislocation can accommodate larger plastic strain with the aid of climb. Since climbed dislocations are not confined on their original slip planes any more, they achieve longer free motion path and thus can produce larger plastic deformation. Moreover, the climb of dislocations in the Al layers is more significant than that in the Cu layers, resulting in more heterogeneous distributions of the dislocation density, the plastic strain, the microscopic stress in different metallic layers, which have substantial influences on the strengthening behavior of MMMFs. 2. With dislocation climb considered, remarkable strain rate dependency of the stress–strain curves is captured. The lower the loading rate is, the more significant the dislocation climb effect on the stress–strain responses is. In addition, the modulus mismatch between different layers tends to enhance the dislocation climb effect on the initial yield stress of MMMFs. 3. There are two mechanisms governing the strengthening behavior of MMMFs. One is the dislocation pile-ups against IBs, and the other is the confined dislocation-source length by the layer thickness. Due to their cooperative contributions to the strengthening behavior of MMMFs, the layer thickness effects on the yield strength of MMMFs are enhanced for both two cases with and without the consideration of dislocation climb, with the layer-thickness exponent n being larger than 0.5 predicted by the classic Hall–-Petch model where only the dislocation pile-up mechanism is considered. Except these two mechanisms, the strengthening mechanism related to the dislocation density seems not to govern strongly the responses of MMMFs with the layer thickness at the submicron scale. 4. Dislocation climb weakens the layer thickness effect on the yield strength of MMMFs. When some dislocations climb successfully away from their pile-ups against IBs, the dislocation pile-ups are relieved and those climbing dislocations achieve longer free motion path. The layer thickness effect induced by dislocation pile-up mechanism is significantly reduced by the dislocation climb at elevated temperature. Due to careful treatments on the interactions between adjacent dislocations and between dislocation and interface, the present 2D-DDD/FEM coupling scheme can capture accurately enough the size effect induced by dislocation pile-ups against various interfaces (such as IBs and GBs) and coated/passivated surfaces; therefore, it can be used to study the size effect due to the microstructural constraint in materials (Arzt, 1998). Even so, an extension of this 2D coupling scheme to 3D is of both academic and practical significance to reveal deeply the micro-constraint induced size effect and its underlying strengthening mechanisms. In addition, the dislocation may cut through the interface by climb at high temperature. Detailed consideration of dislocation cutting through or absorbed/climbing in the coherent IBs, which are also important relaxation mechanisms against the dislocation pile up strengthening mechanism, might be very important to model the high-temperature responses of MMMFs accurately. Considering the focuses of this paper are on the DDD–FEM coupling scheme and

Path B

Slip direction

Edge dislocation i

d pi

Point P

Path A

Fig. A1. Computational model for an edge dislocation at the center of a finite aluminum crystal.

M. Huang, Z. Li / J. Mech. Phys. Solids 85 (2015) 74–97

93

the layer thickness effect induced by competition between dislocation pile up and climb, we neglect the dislocation cutting through the IBs in the present study. The above two limitations will be our aims in the future study.

Acknowledgments The authors would like to express their thanks for the financial supports from NSFC (11272130 and 11272128), Natural Science Foundation of Hubei Province (2013CFA130), and Fundamental Research Funds for the Central Universities (2014TS134). The pertinent comments and good advices given by anonymous reviewers are very helpful to improve this paper.

Appendix A. Validation of the 2D-DDD/FEM coupling scheme To validate the present 2D-DDD/FEM coupling scheme, the stress field of an edge dislocation in a finite aluminum crystal is solved. As shown in Fig. A1, the dislocation is assumed to glide into the crystal from the left side on the horizontal slip plane and arrive at the center of this crystal. The square elements with eight nodes and four GIPs are adopted to discrete this finite aluminum crystal with the dimensions of 2 μm × 2 μm . The model is meshed by 80 × 80 elements, with the element size l = 100b = 0.025 μm . The Burgers vector distribution radius in Eq. (8) is selected as a = 2l = 200b (Liu et al., 2009). In fact, our simulations also suggested that the distribution radius should be selected as two times of the element size as a = 2l . To eliminate the singularity of the FE stiffness matrix, the displacement constraint condition shown in Fig. A1 is applied. In addition, two paths (i.e. the horizontal path A and the vertical path B) are defined, on which the computational stress will be compared with the analytical solution by elastic dislocation theory. The dislocation-induced shear stress σ12 on the path A is plotted in Fig. A2(a) as a function of the normalized coordinate Rpi = 2dpi /a , with dpi being the distance between the dislocation i and the considering point p on path A. The blue solid curve is the result obtained by the pure FE interpolation (i.e. the first term σpFE in Eq. (12)), while the green dashed curve is the solution by classical dislocation theory. It can be clearly seen that, when Rpi ≥ 4 (i.e. dpi ≥ 2a ), two curves are basically overlapped, which means the FE interpolation result is accurate enough on those points away from the dislocation with dpi ≥ 2a . This is the reason why the radii of the cut-off circles cp and ci defined in Eq. (12) are selected as 2a . However, it also can be clearly seen that when Rpi < 4 (i.e. dpi < 2a ), the FE interpolation result σpFE seriously underestimates the actual stress Σp on the point p. Especially, on the core of dislocation i (i.e. dpi = 0), the singular stress is smeared out to 0 by the FEM interpolation. Considering for this, the analytical solution (i.e. the second term of Eq. (12)) is introduced to describe the stress at the point p induced by those dislocations i in the circle cp with dpi < 2a , while the third term shown as the red circles in Fig. A2(a) is specially introduced to remove the contribution of FE-interpolation stress induced merely by those dislocations i within cp (i.e. σ FE ) to avoid repetition of consideration. p i ∈ cp In order to validate Eq. (12), the shear stress σ12 along the path A and the normal stress σ11 along the path B (see Fig. A1) are plotted in Fig. A2 (b) as a function of Rpi , where the theoretical solutions by the classical dislocation theory and the numerical results predicted by Eq. (12) are plotted together for comparison. It is shown that the curves predicted by Eq. (12) are well consistent with the theoretical curves along the paths A and B. Obviously, the introduction of the last two terms in Eq. (12) is very necessary to describe accurately the dislocation stress field and thus the interaction between two close dislocations.

a

b Green dashed: Theoratical value Theoretical value along path A

(MPa)

Blue solid: FE interpolation

(MPa)

Red Circle: The third term in Eq. (12)

By Eq. (12) along path A Theoretical value along path B

Stress

By Eq. (12) along path B

Rpi =1

R pi =4

2d pi /a

2dpi /a

Fig. A2. Comparison between the numerical results and the theoretical solutions for the stress field induced by the dislocation in Fig. A1 (a) without and (b) with the last two terms in Eq. (12).

94

M. Huang, Z. Li / J. Mech. Phys. Solids 85 (2015) 74–97

m

Pa

Y

Y

m

Pa

X

m

X

m

Pa

Y

Y

m

m

Pa

X

X

m

m

Pa

Y

Y

m

m

Pa

X

m

X

m

Fig. A3. The stress contours induced by the dislocation i in Fig. A1. The results in the left row are obtained by the present DCM algorithm, while those in the right row are achieved by the classical dislocation theory.

M. Huang, Z. Li / J. Mech. Phys. Solids 85 (2015) 74–97

a

95

b Pa

Y

Y

m

m

Pa

0

60

X

X

m

m

Fig. A4. (a) Contours of σ11 induced by the dislocation gliding on 60° oblique slip plane of anistropic Al single crystal obtained (a) by present algorithm and (b) by theoretical solution (Shishvan et al., 2011).

a

b

Material 1

Material 2

(MPa)

Material 2

Material 1

Solid: d= 0.0125

m

Dashed: d=0.22

m

d

Slip direction O

x

Red: Theoretical results by Lubarda (1997) Blue: Present numerical results by Eq. (12)

Phase boundary

x

c

d

Pa

2

1

Y

Y

m

m

m)

Pa

2

1

(

X

m

X

m

Fig. A5. (a) Computational model of an edge dislocation near a bi-material interface. (b) The variation of σ12 along the slip plane for two given distances d = {0.0125, 0.22} μm between dislocation and IB. The contours of σ12 obtained (c) by the present algorithm and (d) by Lubarda (1997).

96

M. Huang, Z. Li / J. Mech. Phys. Solids 85 (2015) 74–97

To further validate the present algorithm, the contours of stress fields induced by the dislocation i in Fig. A1 are plotted in Fig. A3, where the results obtained by the present DCM method are given by three figures on the left row (i.e. a1–c1) while the theoretical solutions achieved by the classical dislocation theory are given by three figures on the right row (i.e. a2–c2). It can be seen that the present computational contours agree well with the analytical solutions. It should be pointed out that the stress field contours shown in Fig. A3(a2)–(c2) are the theoretical solutions of the dislocation residing in an infinite homogeneous isotropic material. However, the stress distributions shown in Fig. A3(a1)–(c1) is calculated by the present DCM scheme and Eq. (12), where the influences of the external boundary and the elastic anisotropy on the stress field are included naturally. For these reasons, slight difference between them can be found, especially near the outer boundaries. For the plane strain problem considered, two slip systems are inclined to the element mesh. In addition, the single crystal considered in the MMMFs is elastically anisotropic. It is necessary to verify the present algorithm for dislocations on the oblique slip planes in the anisotropic materials. Considering for this, Fig. A4 gives out the contours of σ11 induced by the dislocation gliding on the 60° inclined slip plane from bottom side of the single crystal Al, with its anisotropic elastic constants listed in Table 1. For comparison, both the numerical contour by present DCM algorithm and the theoretical solution of anisotropic materials (Shishvan et al., 2011) are provided together. It can be clearly seen that, even for the dislocations gliding on the oblique slip planes of anisotropic materials, the present numerical result agrees well with the theoretical one, with slight difference on the outer boundary due to infinite boundary condition considered for the theoretical result. For the MMMFs considered in this study, there exist abundant IBs where dislocation pile-ups can be easily built up. A question, whether the present DDD–FEM coupling scheme can accurately capture the stress field induced by the dislocation near these IBs or not, should be answered first before performing the modeling of MMMFs. For this, a square bimaterial model containing an edge dislocation near a bimaterial interface is considered as shown in Fig. A5(a). This bimaterial model has the same size and mesh as those in Fig. A1. The edge dislocation gliding on the horizontal slip plane from the left side is assumed at the point x = − d . To enhance the influence of bimaterial interface on the dislocation stress field, the elastic properties of these two materials are selected as {E1 = 121.41 GPa, v1 = 0.34} for material 1 and {E2 = 0.1E1, v2 = 0.3} for material 2, respectively, similar to those selected by Gracie and Belytschko (2007). The shear stress σ12 obtained by the present DCM algorithm and by the theoretical solution of infinite bi-materials (Lubarda, 1997) is plotted along the slip plane in Fig. A5(b) for two given distances d = {0.0125, 0.22} μm between dislocation and IB. It can be seen that, for d = 0.22 μm , the computational results obtained by the present DCM scheme are almost overlapped with the theoretical ones both near and away from the dislocation core where the stress singularity appears. This is also true even for the dislocation quite near the interphase boundary with d = 0.0125 μm . In addition, both the numerical and theoretical contours of σ12 for d = 0.22 μm are given out in Fig. A5(c) and (d) for comparison. It shows that the present numerical result has a good agreement with the theoretical solution, both of which can capture the stress jump across the bimaterial interface satisfactorily. Again, the slight difference between them may come from their different boundary conditions. It should be specially noted that, the Lubarda’s (1997) theoretical solution is introduced in Eq. (12) in a cut-off circle with its center at dislocation core and radius rcp = 2a = 0.1 μm for each dislocation, whether it is near or away from the bi-material interface. Further, for the MMMFs with multiple internal interfaces, only the nearest interface to the considered dislocation is taken into account by Lubarda’s (1997) solution in Eq. (12). Above careful comparisons and verifications show that the present 2D-DDD/FEM coupling scheme is accurate enough to calculate the stress field at the point p induced by the dislocation i even though they are very close to each other, and it also can accurately capture the influence of interphase boundaries on the dislocation stress fields. These two points are crucial to model correctly the configuration of pile-ups against the IBs.

References Abdolrahim, N., Zbib, H.M., Bahr, D.F., 2014a. Multiscale modeling and simulation of deformation in nanoscale metallic multilayer systems. Int. J. Plast. 52, 33–50. Abdolrahim, N., Mastorakos, I.N., Shao, S., Bahr, D.F., Zbib, H.M., 2014b. The effect of interfacial imperfections on plastic deformation in nanoscale metallic multilayer composites. Comp. Mater. Sci. 86, 118–123. Akasheh, F., Zbib, H.M., Hirth, J.P., Hoagland, R.G., Misra, A., 2007a. Dislocation dynamics analysis of dislocation intersections in nanoscale multilayer metallic composites. J. Appl. Phys. 101, 084314. Akasheh, F., Zbib, H.M., Hirth, J.P., Hoagland, R.G., Misra, A., 2007b. Interactions between glide dislocations and parallel interfacial dislocations in nanoscale strained layers. J. Appl. Phys. 102, 034324. Ahmed, N., Hartmaier, A., 2013. A two-dimensional dislocation dynamics model of the plastic deformation of polycrystalline metals. J. Mech. Phys. Solids 58 (12), 2054–2064. Arzt, E., 1998. Size effects in materials due to microstructural and dimensional constraints: a comparative review. Acta Mater. 46, 5611–5626. Ayas, C., Deshpande, V.S., Geers, M.G.D., 2012. Tensile response of passivated films with climb-assisted dislocation glide. J. Mech. Phys. Solids 60, 1626–1643. Ayas, C., van Dommelen, J.A.W., Deshpande, V.S., 2014. Climb-enabled discrete dislocation plasticity. J. Mech. Phys. Solids 62, 113–136. Balint, D.S., Deshpande, V.S., Needleman, A., Van der Giessen, E., 2006b. Discrete dislocation plasticity analysis of the wedge indentation of films. J. Mech. Phys. Solids 54, 2281–2303. Bellou, A., Overman, C.T., Zbib, H.M., Bahra, D.F., Misra, A., 2011. Strength and strain hardening behavior of Cu-based bilayers and trilayers. Scr. Mater. 64, 641–644. Benzerga, A.A., Brechet, Y., Needleman, A., Van der Giessen, E., 2004. Incorporating three-dimensional mechanisms into two-dimensional dislocation dynamics. Model. Simul. Mater. Sci. Eng. 12, 159–196. Butrymowicz, D.B., Manning, J.R., Read, M.E., 1973. Diffusion in copper and copper alloys. J. Phys. Chem. Ref. Data 2, 643–655. Cai, W., Arsenlis, A., et al., 2006. A non-singular continuum theory of dislocations. J. Mech. Phys. Solids 54, 561–587. Chu, X., Barnett, S.A., 1995. Model of superlattice yield stress and hardness enhancements. J. Appl. Phys. 77, 4403–4411. Chu, H.J., Wang, J., Beyerlein, I.J., Pan, E., 2013. Dislocation models of interfacial shearing induced by an approaching lattice dislocation. Int. J. Plast. 41, 1–13. Cleveringa, H.H.M., Van der Giessen, E., et al., 2000. A discrete dislocation analysis of mode I crack growth. J. Mech. Phys. Solids 48, 1133–1157.

M. Huang, Z. Li / J. Mech. Phys. Solids 85 (2015) 74–97

97

Danas, K., Deshpande, V.S., 2013. Plane-strain discrete dislocation plasticity with climb-assisted glide motion of dislocations. Model. Simul. Mater. Sci. Eng. 21, 045008. Davoudi, K.M., Nicola, L., Vlassak, J.J., 2012. Dislocation climb in two dimensional discrete dislocation dynamics. J. Appl. Phys. 111, 103522. Dodson, B.W., Tsao, J.Y., 1987. Relaxation of strained-layer semiconductor structures by plastic flow. Appl. Phys. Lett. 51, 1325–1327. Gao, Y., Zhuang, Z., Liu, Z.L., You, X.C., Zhao, X.C., Zhang, Z.H., 2011. Investigations of pipe-diffusion-based dislocation climb by discrete dislocation dynamics. Int. J. Plast. 27, 1055–1071. Gao, Y., Zhuang, Z., You, X.C., 2013. Study of dislocation climb model based on coupling the vacancy diffusion theory with 3d discrete dislocation dynamics. Int. J. Multiscale Comput. Eng. 11, 59–69. Ghoniem, N.M., Han, X., 2005. Dislocation motion in anisotropic multilayer materials. Philos. Mag. 85, 2809–2830. Gracie, R., Belytschko, T., 2007. On XFEM applications to dislocations and interface. Int. J. Plastic. 23, 1721–1738. Greer, J., De Hosson, J.T.M., 2011. Plasticity in small-sized metallic systems: intrinsic versus extrinsic size effect. Prog. Mater. Sci. 56, 654–724. Hafez Haghighat, S.M., Eggeler, G., Raabe, D., 2013. Effect of climb on dislocation mechanisms and creep rates in γ′-strengthened Ni base superalloy single crystals: a discrete dislocation dynamics study. Acta Mater. 61, 3709–3723. Han, X., Ghoniem, N.M., 2005. Stress field and interaction forces of dislocations in anisotropic multilayer thin films. Philos. Mag. 85, 1205–1225. Hirth, J.P., Lothe, J., 1968. Theory of Dislocations. McGraw-Hill, New York. Höchbauer, T., Misra, A., Hattar, K., Hoagland, R.G., 2005. Influence of interfaces on the storage of ion-implanted He in microlayered multilayered metallic composites. J. Appl. Phys. 98, 1–7. Hoagland, R.G., Mitchell, T.E., Hirth, J.P., Kung, H., 2002. On the strengthening effects of interfaces in multilayer fcc metallic composites. Philos. Mag. A 82, 643–664. Hoagland, R.G., Hirth, J.P., Misra, A., 2006. On the role of weak interfaces in blocking slip in nanoscale layered composites. Philos. Mag. 86, 3537–3558. Hoagland, R.G., Kurtz, R.J., Henager Jr., C.H., 2004. Slip resistance of interfaces and the strength of metallic multilayer composites. Scr. Mater. 50, 775–779. Huang, M.S., Li, Z.H., 2013. The key role of dislocation dissociation in the plastic behaviour of single crystal nickel-based superalloy with low stacking fault energy: three-dimensional discrete dislocation dynamics modelling. J. Mech. Phys. Solids 61, 2454–2472. Huang, M.S., Tong, J., Li, Z.H., 2014a. A study of fatigue crack tip characteristics using discrete dislocation dynamics. Int. J. Plastic. 54, 229–246. Huang, M.S., Li, Z.H., Tong, J., 2014b. The influence of dislocation climb on the mechanical behaviour of polycrystals and grain size effect at elevated temperature. Int. J. Plast. 61, 112–127. Hu, S.Y., Henager Jr., C.H., Chen, L.Q., 2010. Simulations of stress-induced twinning and de-twinning: a phase field model. Acta Mater. 58, 6554–6564. Keralavarma, S.M., Benzerga, A.A., 2007. A discrete dislocation analysis of strengthening in bilayer thin films. Model. Simul. Mater. Sci. Eng. 15, S239–S254. Keralavarma, S.M., 2011. A Contribution to the Modeling of Metal Plasticity and Fracture: From Continuum to Discrete Descriptions. Texas A&M University, UMI Dissertations Publishing 3500111. Keralavarma, S.M., Cagin, T., Arsenlis, A., Benzerga, A.A., 2012. Power law creep from discrete dislocation dynamics. Phys. Rev. Lett. 109, 265504. Lemarchand, C., Devincre, B., Kubin, L.P., 2001. Homogenization method for a discrete-continuum simulation of dislocation dynamics. J. Mech. Phys. Solids 49, 1969–1982. Lee, G., Kim, J.-Y., Burek, M.J., Greer, J.R., Tsui, T.Y., 2011. Plastic deformation of indium nanostructures. Mater. Sci. Eng. A 528, 6112–6120. Li, Q., Anderson, P.M., 2005. Dislocation-based modeling of the mechanical behavior of epitaxial metallic multilayer thin films. Acta Mater. 53, 1121–1134. Li, Z.H., Hou, C.T., Huang, M.S., Ouyang, C.J., 2009. Strengthening mechanism in micro-polycrystals with penetrable grain boundaries by discrete dislocation dynamics simulation and Hall–Petch effect. Comput. Mater. Sci. 46, 1124–1134. 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, 1436–1455. Lubarda, V.A., 1997. Energy analysis of dislocation arrays near bimaterial interfaces. Int. J. Solids Struct. 34, 1053–1073. Maier, K., Mehrer, H., Lessmann, E., Schule, W., 1976. Self-diffusion in nickel at low temperatures. Phys. Stat. Sol. (b) 78, 689–698. Mara, N., Sergueeva, A., Misrab, A., Mukherjee, A.K., 2004. Structure and high-temperature mechanical behavior relationship in nano-scaled microlayered multilayered materials. Scr. Mater. 50, 803–806. Mastorakos, I.N., Zbib, H.M., Bahr, D.F., 2009. Deformation mechanisms and strength in nanoscale multilayer metallic composites with coherent and incoherent interfaces. Appl. Phys. Lett. 94, 173114. Matthews, J.W., Crawford, J.L., 1970. Accommodation of misfit between single crystal films of nickel and copper. Thin Solid Films 5, 187–198. Matthews, J.W., Blakeslee, A.E., 1974. Defects in epitaxial multilayers. I. Misfit dislocations. J. Cryst. Growth 27, 118–125. Misra, A., Kung, H., 2001. Deformation behavior of nanostructured metallic multilayers. Adv. Eng. Mater. 3, 217–222. Misra, A., Hirth, J.P., Kung, H., 2002. Single-dislocation-based strengthening mechanisms in nanoscale metallic multilayers. Philos. Mag. 82, 2935–2951. Misra, A., Hirth, J.P., Hoagland, R.G., 2005. Length-scale-dependent deformation mechanisms in incoherent metallic microlayered multilayered composites. Acta Mater. 53, 4817–4824. Misra, A., Demkowicz, M.J., Zhang, X., Hoagland, R.G., 2007. The radiation damage tolerance of ultra-high strength nanolayered composites. JOM 59, 62–65. Mordehai, D., Clouet, E., Fivel, M., Verdier, M., 2008. Introducing dislocation climb by bulk diffusion in discrete dislocation dynamics. Philos. Mag. 88 (6), 899–925. Ohashi, T., Kawamukai, M., Zbib, H., 2007. A multiscale approach for modeling scale-dependent yield stress in polycrystalline metals. Int. J. Plast. 23, 897–914. Rice, J.R., 1987. Tensile crack tip fields in elastic-ideally plastic crystals. Mech. Mater. 6, 317–335. Schwarz, K.W., Tersoff, J., 1996. Interaction of threading and misfit dislocations in a strained epitaxial layer. Appl. Phys. Lett. 69, 1220–1222. Schwarz, K.W., 1999. Simulation of dislocations on the mesoscopic scale. II. Application to strained-layer relaxation. J. Appl. Phys. 85, 120–129. Shao, S., Medyanik, S.N., 2010. Interaction of dislocations with incoherent interfaces in nanoscale FCC–BCC metallic bi-layers. Model. Simul. Mater. Sci. Eng. 18, 055010 . Shao, S., Zbib, H.M., Mastorakos, I., Bahr, D.F., 2013. Effect of interfaces in the work hardening of nanoscale multilayer metallic composites during nanoindentation: a molecular dynamics investigation. J. Eng. Mater. Technol. Trans. ASME 135, 021001. Shishvan, S.S., Mohammadi, S., Rahimian, M., Van der Giessen, E., 2011. Plane-strain discrete dislocation plasticity incorporating anisotropic elasticity. Int. J. Solids Struct. 48, 374. Van der Giessen, E., Needleman, A., 1995. Discrete dislocation plasticity: a simple planar model. Model. Simul. Mater. Sci. Eng. 3, 689. Vattré, 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 method: the Discrete-Continuous Model revisited. J. Mech. Phys. Solids 63, 491–505. Wang, Y.C., Misra, A., Hoagland, R.G., 2006. Fatigue properties of nanoscale Cu/Nb multilayers. Scr. Mater. 54, 1593–1598. Wang, J., Zhang, R.F., Zhou, C.Z., Beyerlein, I.J., Misra, A., 2014. Interface dislocation patterns and dislocation nucleation in face-centered-cubic and body-centeredcubic bicrystal interfaces. Int. J. Plast. 53, 40–55. Was, G.S., Foecke, T., 1996. Deformation and fracture in microlaminates. Thin Solid Films 286, 1–31. Yang, H., Huang, M.S., Li, Z.H., 2015. The influence of vacancies diffusion-induced dislocation climb on the creep and plasticity behaviours of nickel-based single crystal superalloy. Comp. Mater. Sci. 99, 348–360. Zbib, H.M., Rhee, M., Hirth, J.P., 1998. On plastic deformation and the dynamics of 3D dislocations. Int. J. Mech. Sci. 40, 113–127. Zbib, H.M., de la Rubia, T.D., 2002. A multiscale model of plasticity. Int. J. Plast. 18, 1133–1163. Zbib, H.M., Overman, C., Akasheh, F., Bahr, D.F., 2011. Analysis of plastic deformation in nanoscale metallic multilayers with coherent and incoherent interfaces. Int. J. Plast. 27, 1618–1638. Zhu, Y.X., Li, Z.H., Huang, M.S., 2014. The size effect and plastic deformation mechanism transition in the nanolayered polycrystalline metallic multilayers. J. Appl. Phys. 115, 233508.