Interface-Dominated Plasticity and Kink Bands in Metallic Nanolaminates

Abstract

The theoretical and computational framework of finite deformation mesoscale field dislocation mechanics (MFDM) is used to understand the salient aspects of kink-band formation in Cu-Nb nano-metallic laminates (NMLs). A conceptually minimal, plane-strain idealization of the three-dimensional geometry, including crystalline orientation, of additively manufactured NML is used to model NMLs. Importantly, the natural jump/interface condition of MFDM imposing continuity of (certain components) of plastic strain rates across interfaces allows theory-driven ‘communication’ of plastic flow across the laminate boundaries in our finite element implementation. Kink bands under layer parallel compression of NMLs in accord with experimental observations arise in our numerical simulations. The possible mechanisms for the formation and orientation of kink bands are discussed, within the scope of our idealized framework. We also report results corresponding to various parametric studies that provide preliminary insights and clear questions for future work on understanding the intricate underlying mechanisms for the formation of kink bands.

1. Introduction

A fundamental challenge in the development of efficient engineering materials is to design a material that has both high strength and high toughness, as most materials have exclusively either of these properties. Kinking in nano-metallic laminates is a deformation mechanism that may circumvent this challenge, as it leads to enhanced plasticity (ductility) at an almost constant level of stress after yielding [1]. Kink bands, in the present context of a material with an intrinsic layered structure, refer to a thin band of localized deformation that is transverse—and most commonly not normal—to the intrinsic layering, as shown in Figure 1a. Shear band formation has been the predominantly observed form of localization in NMLs, as reported in [2,3]. In these studies, NMLs are compressed perpendicular to the layer direction in micropillar configurations, meaning the compression loading direction is parallel to the layer normal direction. While this is a somewhat physically intuitive mode of deformation localization through rotation of the layers towards the loading axis and slipping along softer layers, kink band formation in NMLs for sufficiently small layer thickness under layer parallel compression has also been demonstrated [4]. In layer parallel compression of NML, the compression loading direction is perpendicular to the layer normal direction, as shown in Figure 1a. The formation of such kink bands and its dependence on the size of layers is quite non-intuitive, and it is the study of this phenomenon that forms the subject of this paper; the modeling of the development of shear bands in NMLs in micropillar configurations under layer perpendicular compression is a subject worthy of study on its own merits, a challenge that has not been successfully solved as yet [5].

In general, the formation of shear or kink bands, due to strain localization, is detrimental as it produces sites of narrow planes under intense shear or lattice rotations, which act as a source for crack initiation and lead to eventual material failure. For instance, it was reported in [6] that for grain sizes smaller than 300 nm in iron (Fe), shear bands are formed leading to strain-softening and failure at low values of plastic strain. Similarly, graphite-fiber-reinforced epoxy composites fail under compression loading due to the formation of kink bands [7]. However, for layer parallel compression of NMLs with sufficiently small laminate widths, a large amount of plastic deformation is observed before failure without much reduction in the strength of the material in the experimental work of Nizolek et al. [4] in Cu-Nb nanolaminates. It was reported there that the samples with layer thickness greater than 250 nm deform homogeneously, while the deformation behavior of samples with thicknesses 65 or 30 nm was inhomogeneous, and enhanced plastic deformation was observed after yield at a more-or-less constant load [4]. Recently, Zhang et al. [8] investigated the kink band formation process in Cu/Nb NMLs via in situ micropillar compression along the layer-parallel direction in a scanning electron microscope (SEM). It was reported that after the onset of kinking, dislocations active on the slip systems nearly parallel to the interface in NMLs are responsible for the band evolution. Cahn [9] showed the formation of kink bands in aluminum single crystals, and it was suggested that these bands are composed of dipolar arrays of edge dislocations aligned perpendicular to the primary slip direction.

Figure 1. (a) Schematic of an idealized kink band in layer parallel compression of NML, (b) geometry of the kink band at $|e_N|=25\%$ engineering strain in the experimental work of compression of Cu/Nb NMLs. (Figure in (b) reprinted from [8] with permission from Elsevier.).

Figure 1. (a) Schematic of an idealized kink band in layer parallel compression of NML, (b) geometry of the kink band at $|e_N|=25\%$ engineering strain in the experimental work of compression of Cu/Nb NMLs. (Figure in (b) reprinted from [8] with permission from Elsevier.).

With respect to modeling, Zecevic et al. [5] used Gurtin’s [10] strain gradient crystal plasticity model augmented with various choices of the defect energy function within a large strain elasto-viscoplastic Fast Fourier transform-based computational technique to model the compression of Cu/Nb NMLs. Arora and Acharya [11] modeled a 2-D single crystal under simple shearing using the mesoscale field dislocation mechanics (MFDM) framework in a finite deformation setting. In a special single slip study, kink bands were observed in the form of dipolar dislocation walls perpendicular to the single slip plane [11]. Kink band formation in a polycrystalline Mg alloy under compression was computationally studied using a strain-gradient crystal plasticity framework [12]. Asaro and Rice [13] described the onset of a shear and kink band from a bifurcation analysis of the classical, finite deformation, single slip constitutive behavior. Forest [14] used small deformation generalized crystal plasticity (Cosserat theory of crystal plasticity) to separate bifurcation modes (slip and kink bands), i.e., kink bands occur later than slip bands, and this is unlike the prediction of classical plasticity.

In this work, we model the compression of composite nanolaminates within a plane strain idealization, using the theoretical and computational framework of MFDM. The interface conditions in MFDM impose continuity of certain components of the total plastic strain rate, which models continuous plastic flow across the laminate interfaces. We consider a ‘minimal’ model setup in considering 3 co-planar slip systems for each material in the bi-metallic laminate structure and show that kink bands are formed in qualitative accord with the experimental observations of [4,8]. We also report on other studies within this idealized setup and discuss whether kink bands form or not for different cases.

The outline of the paper is as follows: after presenting some notational details used in this paper in the following paragraph, Section 2 provides a brief review of the theory of MFDM. For more details on the theoretical and computational framework of MFDM, the interested reader is referred to [11,15,16,17]. Section 3 discusses the various successes of the MFDM framework in correspondence with experimental studies in many instances. Section 4 shows the numerical results of kink-band formation in composite nanolaminates within our idealization using the computational framework of MFDM. We also discuss various parametric studies of our model characterizing this complicated problem of kink band formation in NMLs. Section 5 presents some concluding remarks and the scope of future work for modeling compression of NMLs using our framework.

Vectors and tensors are represented by boldface lower and upper-case letters, respectively. The action of a second-order tensor A on a vector b is denoted by $\mathit{A}\mathit{b}$. The inner product of two vectors is denoted by $\mathit{a}·\mathit{b}$, while the inner product of two second-order tensors is denoted by $\mathit{A}:\mathit{B}$. A rectangular Cartesian coordinate system is invoked for the ambient space and all vector and tensor components are expressed with respect to the basis of this coordinate system. The partial derivative of the quantity $(·)$ w.r.t. the $x_i$ coordinate of this coordinate system is denoted by ${(·)}_{,i}$, and $e_i$ denotes the unit vector in the $x_i$ direction. The time derivative of a quantity $(·)$ is denoted by $\dot{(·)}$. Einstein’s summation convention is always implied unless mentioned otherwise. The symbols $\mathrm{grad}$, $\mathrm{div}$, and $\mathrm{curl}$ denote the gradient, divergence, and curl on the current configuration, respectively. For a second order tensor A, vectors v, a, and c, and a spatially constant vector field b, the operations of $\mathrm{div}$, $\mathrm{curl}$, and the cross-product of a tensor $(\times )$ with a vector are defined as follows:

$$\begin{array}{cc}\hfill \left(\mathrm{div}\mathit{A}\right)·\mathit{b}=\mathrm{div}\left({\mathit{A}}^T\mathit{b}\right),& \forall \mathit{b}\hfill \\ \hfill \mathit{b}·\left(\mathrm{curl}\mathit{A}\right)\mathit{c}=\left[\mathrm{curl}\left({\mathit{A}}^T\mathit{b}\right)\right]·\mathit{c},& \forall \mathit{b},\mathit{c}\hfill \\ \hfill \mathit{c}·\left(\mathit{A}\times \mathit{v}\right)\mathit{a}=\left[\left({\mathit{A}}^T\mathit{c}\right)\times \mathit{v}\right]·\mathit{a},& \forall \mathit{a},\mathit{c}.\hfill \end{array}$$

In rectangular Cartesian coordinates, these are denoted by

$$\begin{array}{c}\hfill {\left(\mathrm{div}\mathit{A}\right)}_i=A_{ij,j},{\left(\mathrm{curl}\mathit{A}\right)}_{ri}={\epsilon }_{ijk}{A}_{rk,j},{\left(\mathit{A}\times v\right)}_{ri}={\epsilon }_{ijk}{A}_{rj}{v}_k,\end{array}$$

where ${\epsilon }_{ijk}$ are the components of the third order alternating tensor X.

When we refer to different slip systems with respect to their slip vector orientation, we mean that those orientations were in the initial unstretched lattice, and do not represent the orientations of the corresponding slip vectors in the deformed lattice, unless mentioned otherwise.

2. Review of MFDM Theory

The governing equations for MFDM [15] are as follows:

$$\begin{array}{c}\mathrm{curl}\mathit{W}=-\mathit{\alpha },\end{array}$$

(3a)

$$\begin{array}{c}\dot{\mathit{W}}+\mathit{W}\mathit{L}=\mathit{\alpha }\times \mathit{V},\end{array}$$

(3b)

$$\begin{array}{c}\mathrm{div}\left[\mathit{T}\left(\mathit{W}\right)\right]=\{\begin{array}{c}\mathbf{0},\mathrm{quasistatic}\hfill \\ \rho \dot{v},\mathrm{dynamic}.\hfill \end{array}\end{array}$$

(3c)

In Equation (3), $\mathit{W}$ is the inverse of the elastic distortion tensor (${\mathit{F}}^e$), $\mathit{\alpha }$ denotes the dislocation density tensor (Nye’s tensor), v denotes the material velocity, $\mathit{L}=\mathrm{grad}\mathit{v}$ is the velocity gradient, T is the (symmetric) Cauchy stress tensor, and V denotes the dislocation velocity field with respect to the material. Equation (3a) is the incompatibility equation for the inverse elastic distortion tensor, while Equation (3b) is derived from the statement of conservation of the Burgers vector content of an arbitrary area patch in the current configuration [18]. Equation (3c) is the linear momentum balance for both quasistatic and dynamic cases.

Equation (3b) is a system of coupled Hamilton-Jacobi equations for which, to our knowledge, no rigorous mathematical theory or numerical schemes exist in the literature. To computationally approximate the governing equations, a Stokes-Helmholtz decomposition of the field W is introduced, expressing it as a sum of a gradient $\left(\mathrm{grad}\mathit{f}\right)$ and a divergence-free part ($\mathit{\chi }$) [11]. This decomposition replaces the tensor field W with a tensor and a vector field in our computations, but it is relatively easier to develop numerical schemes for this larger set of governing equations [15].

The extended set of governing equations for MFDM [15] are shown below:

$$\begin{array}{c}\mathring{\mathit{\alpha }}\equiv \left(\mathrm{div}\left(\mathit{v}\right)\right)\mathit{\alpha }+\dot{\mathit{\alpha }}-\mathit{\alpha }{\mathit{L}}^T=-\mathrm{curl}(\mathit{\alpha }\times \mathit{V}+{\mathit{L}}^p),\end{array}$$

(4a)

$$\begin{array}{c}\begin{array}{c}\mathit{W}=\mathit{\chi }+\mathrm{grad}\mathit{f},{\mathit{F}}^e={\mathit{W}}^{-1},\\ \mathrm{curl}\mathit{\chi }=-\mathit{\alpha },\\ \mathrm{div}\mathit{\chi }=\mathbf{0},\end{array}\}\end{array}$$

(4b)

$$\begin{array}{c}\mathrm{div}\left(\mathrm{grad}\dot{\mathit{f}}\right)=\mathrm{div}\left(\mathit{\alpha }\times \mathit{V}+{\mathit{L}}^p-\dot{\mathit{\chi }}-\mathit{\chi }\mathit{L}\right),\end{array}$$

(4c)

$$\begin{array}{c}\mathrm{div}\left[\mathit{T}\left(\mathit{W}\right)\right]=\{\begin{array}{c}\mathbf{0},\mathrm{quasistatic}\hfill \\ \rho \dot{v},\mathrm{dynamic}.\hfill \end{array}\end{array}$$

(4d)

Here, Equation (4a) is the evolution equation for the dislocation density tensor, and (4b) is the incompatibility equation in terms of $\mathit{\chi }$ with a condition that $\mathrm{div}\mathit{\chi }=\mathbf{0}$, while f denotes the plastic position vector, and the evolution equation for f is given by (4c). ${\mathit{L}}^p$ is the additional mesoscale field which represents the averaged rate of plastic straining due to all dislocations that cannot be represented by $\mathit{\alpha }\times \mathit{V}$, where both fields in the product represent space-time running averages of the corresponding microscopic fields. This is due to the fact that the average of products is not the same as the product of averages, in general. Hence, all fields in the MFDM framework are running space-time averages of the corresponding fields of the field dislocation mechanics (FDM) theory; except for ${\mathit{L}}^p$, as there is no corresponding field in FDM theory to ${\mathit{L}}^p$.

To close the above set of equations in (4), constitutive statements for $\mathit{V},{\mathit{L}}^p,\mathit{T}$ are used, consistent with the mechanical dissipation being non-negative [11].

Initial conditions for $\mathit{\alpha },\mathit{f}$ and boundary conditions for $\mathit{\alpha },\mathit{f},\mathit{\chi }$, and v are specified, for a well-set evolution in MFDM.

2.1. Boundary Conditions

• The $\mathit{\alpha }$ evolution equation (shown in Equation (4a)) has a convective boundary condition of the form $(\mathit{\alpha }\times \mathit{V}+{\mathit{L}}^p)\times \mathit{n}=\mathit{\Phi }$, where $\mathit{\Phi }$ is a second order tensor valued function of position and time on the boundary, characterizing the flux of dislocations at the surface satisfying the constraint $\mathit{\Phi }\mathit{n}=\mathbf{0}$. Here, n is the outward unit normal field on the boundary.
  The boundary condition is specified in two ways: (a) Plastically constrained: $\mathit{\Phi }(x,t)=\mathbf{0}$ is specified at a point x on the boundary for all times, which ensures that there is no outflow of dislocations at that point of the boundary, and only parallel motion along the boundary is allowed. (b) Plastically unconstrained: A less restrictive boundary condition where ${\widehat{\mathit{L}}}^p\times \mathit{n}$ is simply evaluated at the boundary (akin to an outflow condition), along with the specification of dislocation flux $\mathit{\alpha }(\mathit{V}·\mathit{n})$ on the inflow part of the boundary. Additionally, for all calculations presented in this paper $(\mathrm{curl}\mathit{\alpha }\times \mathit{n})=\mathbf{0}$ is imposed, a particular specification of a boundary condition that arises from simple mathematical modeling of the manifestation of dislocation core energy at the mesoscale.
• For the incompatibility equation, $\mathit{\chi }\mathit{n}=\mathbf{0}$ is applied on the outer boundary of the domain, which along with the system (4b) ensures that $\mathit{\chi }$ vanishes when $\mathit{\alpha }$ is zero in the entire domain.
• The $\mathit{f}$ evolution equation requires a Neumann boundary condition, i.e., $\left(\mathrm{grad}\dot{\mathit{f}}\right)\mathit{n}=(\mathit{\alpha }\times \mathit{V}+{\mathit{L}}^p-\dot{\mathit{\chi }}-\mathit{\chi }\mathit{L})\mathit{n}$ on the outer boundary of the domain.
• The material velocity boundary conditions are applied based on the loading type, which is discussed later in Section 4.

2.2. Initial Conditions

• The initial condition $\mathit{\alpha }(\mathit{x},0)=\mathbf{0}$ is assumed for all simulations here.
• In general, the initial condition for $\mathit{f}$ is obtained by solving for $\mathit{\chi }$ from the incompatibility equation and solving for $\mathit{f}$ from the equilibrium equation, for prescribed $\mathit{\alpha }$ on the given initial configuration. We refer to this scheme as the elastic theory of continuously distributed dislocations (ECDD). For the initial conditions on $\mathit{\alpha }$ considered above, this step is trivial, with $\mathit{f}=\mathit{X}$, where X is the position field on the initial configuration.
• The model admits an arbitrary specification of $\dot{\mathit{f}}$ at a point to uniquely evolve f using (4c) in time, and this rate is prescribed to vanish.

Interface conditions for plastic distortion rate: The interface jump condition for the $\mathit{\alpha }$ solve is $\left[[\mathit{\alpha }\times \mathit{V}+{\widehat{\mathit{L}}}^p+l^2\widehat{\gamma }\mathrm{curl}\mathit{\alpha }]\right]\times \mathit{n}=\mathbf{0}$, and this gets imposed in the finite element formulation. Here, $\left[\right[(·)\left]\right]$ denotes the jump in $(·)$ across the interface. This interface condition allows for the discontinuity of certain components of plastic strain rate across the interface [19].

Traction continuity and natural interface conditions consistent with the PDE for $\mathit{\chi },\dot{\mathit{f}}$ are also enforced by our finite element formulation.

2.3. Constitutive Relations

Constitutive relations in MFDM are required for the stress T, the plastic distortion rate ${\mathit{L}}^p$, and the dislocation velocity V. The details of the thermodynamically consistent constitutive formulations can be found in Section 3.1 of [11].

Table 1. Constitutive relations for Cauchy stress and core energy density.

Saint-Venant-Kirchhoff Material	${\displaystyle \varphi \left(\mathit{W}\right)=\frac{1}{2{\rho }^{*}}{\mathit{E}}^{\mathit{e}}:\mathbb{C}:{\mathit{E}}^{\mathit{e}},\mathit{T}={\mathit{F}}^{\mathit{e}}[\mathbb{C}:{\mathit{E}}^{\mathit{e}}]{{\mathit{F}}^{\mathit{e}}}^{\mathit{T}}}$
Core energy density	${\displaystyle {\rm Y}\left(\mathit{\alpha }\right)=\frac{1}{2{\rho }^{*}}ϵ\mathit{\alpha }:\mathit{\alpha }}$

presents the constitutive relation for Cauchy stress and mesoscopic core energy density for the material. Here, ${\rho }^{*}$ denotes the mass density of the pure, un-stretched elastic lattice, and $\mathbb{C}$ is the fourth-order elasticity tensor, while ${\mathit{E}}^e$ is the elastic Lagrangian strain tensor. ${\rm Y}\left(\mathit{\alpha }\right)$ is a heuristic representation of the averaged microscopic core energy at the mesoscale, and $ϵ$ is a constant that has physical dimensions of $stress\times {length}^2$.

Table 2. Constitutive relations for plastic strain rate ${\mathit{L}}^p$.

Crystal plasticity

${\displaystyle {\widehat{\mathit{L}}}^{\mathit{p}}=\mathit{W}{\left(\sum _{\mathit{k}}^{{\mathit{n}}_{\mathit{s}\mathit{l}}}{\widehat{\gamma }}^{\mathit{k}}{\mathit{m}}^{\mathit{k}}\otimes {\mathit{n}}^{\mathit{k}}\right)}_{\mathit{sym}};{\widehat{\gamma }}^{\mathit{k}}=\mathit{sgn}\left({\tau }^{\mathit{k}}\right){\widehat{\gamma }}_0^{\mathit{k}}{\left(\frac{|{\tau }^{\mathit{k}}|}{\mathit{g}}\right)}^{\frac{1}{\mathit{m}}}}$

${\displaystyle {\mathit{L}}^p={\widehat{\mathit{L}}}^p+\frac{l^2}{n_{sl}}\sum _k^{n_{sl}}\left|{\widehat{\gamma }}^k\right|\mathrm{curl}\mathit{\alpha }}$

shows the constitutive relations for plastic distortion rate. Classical crystal plasticity is assumed here and the plastic distortion ${\mathit{L}}^p$ is given by the sum of slipping distortions on prescribed slip systems [11]. In Table 2, ${(·)}_{sym}$ denotes the symmetric part of $(·)$, ${\widehat{\gamma }}_0^k$ is a reference strain rate, ${\widehat{\gamma }}^k$ is the slipping rate due to statistically stored dislocations on the $k^{th}$ slip system, $n_{sl}$ is the total number of slip system, $sgn\left({\tau }^k\right)$ denotes the sign of the scalar ${\tau }^k$, g is the material strength, m is the material rate sensitivity. ${\mathit{m}}^k$ and ${\mathit{n}}^k$ represents the slip direction and slip plane normal for the $k^{th}$ slip system, in the current configuration. The resolved shear stress on the $k^{th}$ slip system is given as ${\tau }^k={\mathit{m}}^k·\mathit{T}{\mathit{n}}^k$. The slip direction and slip plane normal in the current configuration are related to the unstretched slip direction (${\mathit{m}}_0^k$) and the slip plane normal (${\mathit{n}}_0^k$) by ${\mathit{m}}^k={\mathit{F}}^e{\mathit{m}}_0^k$, and ${\mathit{n}}^k={{\mathit{F}}^e}^{-T}{\mathit{n}}_0^k$. The classical plastic distortion rate ${\widehat{\mathit{L}}}^p$ from crystal plasticity, is augmented with an additive term in $\mathrm{curl}\mathit{\alpha }$, motivated from the mechanical dissipation expression in [11] [Equation (8)], with a mobility coefficient given by $l^2$ multiplied by the average slipping rate of the total slip systems [11]. Here, l is a material length related to the gross modeling of mesoscale effects of dislocation core energy and is defined as $l^2=ϵ/g_0$, with $g_0$ being the initial strength (initial yield stress in shear).

Table 3. Constitutive relations for dislocation velocity V.

${\displaystyle {\mathit{T}}^{{}^{\prime }}=\mathit{T}-\frac{\mathit{t}\mathit{r}\left(\mathit{T}\right)}{3}\mathit{I}}$;	${\displaystyle \mathit{a}=\frac{1}{3}\mathit{X}\left(\mathit{t}\mathit{r}\left(\mathit{T}\right){\mathit{F}}^{\mathit{e}}\mathit{\alpha }\right);}$	${\displaystyle \mathit{c}=\mathit{X}\left({\mathit{T}}^{{}^{\prime }}{\mathit{F}}^{\mathit{e}}\mathit{\alpha }\right)}$
${\displaystyle \mathit{d}=\mathit{c}-\left(\mathit{c}-\frac{\mathit{a}}{\left|\mathit{a}\right|}\right)\frac{\mathit{a}}{\left|\mathit{a}\right|};}$	${\displaystyle \mathit{V}=\zeta \frac{\mathit{d}}{\left|\mathit{d}\right|};}$	${\displaystyle \zeta =\frac{{\mu }^2{\eta }^{2}b}{g^2{n}_{sl}}\sum _k^{n_{sl}}\left|{\widehat{\gamma }}^k\right|}$

shows the constitutive relations for dislocation velocity. Here, d denotes the direction of dislocation velocity, and its magnitude is given by $\zeta$, ${\mathit{T}}^{\prime }$ denotes the deviatoric stress tensor, X is the third-order alternating unit tensor, $\mu$ is the shear modulus, g is the material strength, $\eta$ is a non-dimensional material constant in the empirical Taylor relationship for macroscopic strength vs. dislocation density, and b is the Burgers vector magnitude of a full dislocation in the crystalline material. For the motivation behind the constitutive statement for V, the interested reader is referred to [20,21].

Table 4. Constitutive relations for material strength g.

${\displaystyle \dot{\mathit{g}}=\mathit{h}(\mathit{\alpha },\mathit{g})\left(|{\mathit{F}}^{\mathit{e}}\mathit{\alpha }\times \mathit{V}|+\sum _{\mathit{k}}^{{\mathit{n}}_{\mathit{s}\mathit{l}}}\left|{\widehat{\gamma }}^{\mathit{k}}\right|\right);}$

${\displaystyle \mathit{h}(\mathit{\alpha },\mathit{g})=\frac{{\mathit{\mu }}^2{\mathit{\eta }}^2\mathit{b}}{2(\mathit{g}-{\mathit{g}}_0)}{\mathit{k}}_0\left|\mathit{\alpha }\right|+{\mathit{\Theta }}_0\left(\frac{{\mathit{g}}_{\mathit{s}}-\mathit{g}}{{\mathit{g}}_{\mathit{s}}-{\mathit{g}}_0}\right)}$

shows the evolution equation for material strength (g). Here, $g_s$ denotes the saturation strength, h is the total hardening rate, ${\Theta }_0$ is the Stage II hardening rate, $k_0$ (non-dimensional) characterizes the hardening rate due to geometrically necessary dislocations (GNDs) [22].

All parameters in our model, except $k_0$ and l, are part of the constitutive structure of well-accepted models of classical crystal plasticity theory. The parameter $k_0$ was introduced in [22]. The length scale l simply controls the refinement of the GND microstructure and does not play a physically significant role in our results.

3. Review of Results from MFDM

The MFDM framework is based on the space-time averaging of microscopic dislocation mechanics, and the obtained governing PDEs shown in (4) are closed by the phenomenological specification of ${\mathit{L}}^p$, V, T and g. This approach has been quite successful in addressing some challenging problems in modern plasticity theory related to the computation of patterning [11], dislocation internal stress [15], size effects in micropillar confined thin metal films [23], polygonization [16], and slip transmission at grain boundaries [24,25] among others [18,19,20,26,27,28,29,30,31,32,33,34,35,36,37,38,39], including long-standing and recent fundamental challenges in the prediction of large-deformation, dislocation mediated elastic and elastic-plastic response [11,15,16,23]. The benefits of using such a framework are as follows:

• It is difficult to deduce a physical connection between the plastic strain/distortion in classical plasticity theory to the mechanics of dislocation, beyond modeling in 1-D. The MFDM framework brings out an explicit connection between the plastic strain rate and the motion and geometry of an evolving microscopic array of dislocations. Obviously, this has many benefits, even for a phenomenological specification of the macroscopic plastic strain rate. Moreover, the MFDM framework has allowed for a first unification between phenomenological ${\mathit{J}}_2$ and crystal plasticity theories and quantitative dislocation mechanics.
• With a single extra material fitting parameter beyond a classical plasticity model (and two in the finite deformation setting), the MFDM framework has enabled a significant variety of phenomena to be modeled, in qualitative and quantitative accord with experimental results.

4. Results and Discussion

The initial configuration of NMLs, as shown in Figure 2a, is such that there are 12 alternating layers of Cu and Nb. Each layer has a thickness of $l=80$ nm with the width of the composite being $W=960$ nm. The aspect ratio of the NML composite is 1:2, with an undeformed height of the sample being $L=1920$ nm. The nominal compression loading rate is $|\widehat{e}|=0.001s^{-1}$, and at any time t, the nominal compressive strain is $\left|e\right|=|\widehat{e}|t$. The corresponding true compressive strain is denoted as $|e_T|$. The average true compressive stress is denoted by $\sigma$, while the average nominal compressive stress is denoted by ${\sigma }_N$. Plastically unconstrained boundary conditions are applied to all four (external) boundaries of the domain, for all results shown in this work. Quasi-static force equilibrium is solved for all simulations reported here.

The boundary conditions for the material velocity for the simulations of compression of NMLs are as follows:

• $v_2=0$ at the bottom boundary of the domain.
• $v_2=-\left|\widehat{e}\right|L$ at the top boundary of the domain, where L is the length of the laminate in the (undeformed) initial configuration at $t=0$.
• The applied traction in the horizontal direction is zero on both the top and the bottom boundary of the domain.
• $v_1=0$ at a single node of the bottom boundary. This along with the above conditions suffices to constrain rigid motion.

As discussed earlier, we assume 3 slip systems for both Cu and Nb. One of the reference slip directions is aligned along the layer parallel direction, while the other two are defined to initially accommodate symmetric double slip with respect to the loading axis. With reference to Figure 2a, ${\theta }_1$ and ${\theta }_2$ are both assumed to be ${30}^{\circ }$ for Cu, while for Nb, both are ${45}^{\circ }$, unless mentioned otherwise.

Table 5. Material parameters for Cu and Nb.

Parameter	${\widehat{\gamma }}_0$	m	$\eta$	b	${\mathit{g}}_0$	${\mathit{g}}_s$	${\mathit{\Theta }}_0$	${\mathit{k}}_0$	l	E	$\mathit{\nu }$
	($s^{-1}$)			(Å)	(MPa)	(MPa)	(MPa)		($\mathsf{\mu }$m)	(GPa)
Cu	0.001	0.03	$\frac{1}{3}$	2.556	210	1890	273	20	$\sqrt{3}\times 0.1$	144.58	0.324
Nb	0.001	0.03	$\frac{1}{3}$	2.86	262.5	2362.5	198	20	$\sqrt{3}\times 0.1$	110.25	0.392

shows the material parameters used for the numerical simulations of NMLs under compression in this study. The $g_0$ values for Cu and Nb are taken from [5], and the saturation strength is assumed to be $g_s=9g_0$. We use these material properties for all of our simulations unless otherwise stated, although these values of initial and saturation yield strength for single crystal materials seem somewhat out of the norm. Accordingly, we also perform simulations with more commonly used values without a qualitative change in results.

4.1. Kink Band Formation

Figure 2b,c shows the true and engineering stress-strain curve obtained from our simulation and the experimentally obtained stress-strain curve from [8], and it is clear that the numerically obtained stress-strain curves do not match with the experimental ones.

The yield strength obtained from our simulations is between 0.45–0.5 GPa, while in the experimental results from [8], it is close to 1.0 GPa. The strain-hardening response is more prominent in our numerical result as compared to the experimental observation, where after yielding it stays at a more-or-less constant level of stress. However, in the current work or our past works like [16,23], we have not employed any fitting of material parameters to obtain the numerical results that are in qualitative or quantitative agreement with experimental observations. This is discussed in a broader context in Section 4.6, where we compare strain gradient plasticity (SGP) based models and MFDM for modeling problems in mesoscale plasticity.

Figure 3 shows the field plot of geometrically necessary dislocation (GND) density (${\rho }_g=\left|\mathit{\alpha }\right|/b\left(m^{-2}\right)$) and the norm of logarithmic strain tensor $\left|ln\mathit{V}\right|$. Here, V is obtained from the polar decomposition of the deformation gradient tensor ($\mathit{F}=\mathit{V}\mathit{R}$), and F is defined with respect to the initial configuration at $t=0$. GNDs are generated in the domain due to the inhomogeneity in the material response after yielding from the material inhomogeneity of the laminate structure in NMLs, and the imposition of jump condition on plastic strain rate at the interfaces.

It is evident from Figure 3 that the kink-band forms in accordance with the experimental observations in [4,8]. It is evident that the deformation field localizes within the band, and the kink band is inclined (not perpendicular) to the layering in NML as shown in Figure 1. Also, as per Figure 3b that the principal stretches of V are maximum within the kink band formed, due to intense deformation in the band region.

The inclination of the kink band shown in Figure 3 with respect to the layer normal (in the initial configuration), i.e., $\beta$ (shown in Figure 1a) is estimated roughly to be ${32}^{\circ }$, and, curiously, it is close to $\beta ={35}^{\circ }$ reported in the experimental work in [8]. The angle $\beta$ of the inclination of the kink band is calculated based on the actual deformed coordinates of the kink band boundary (note that the aspect ratio of the plot box shown in Figure 3 is close to 1:1).

Figure 4 shows the comparison of the stress-strain response with different values of $g_0$ (initial yield), but with the same ratio between the saturation yield strength and initial yield i.e., $g_s=9g_0$. As expected, the initial yield strength for the NML is less, with reduced values of $g_0$, as compared to the experimentally observed yield limit [8], but the ultimate strength, in this case, is closer to the experimental value. Importantly, we observe the formation of kink bands with similar shapes in our numerical simulations for smaller values of $g_0$ and $g_s$, as well.

The orientation of the kink band depends on the Schmid factors of the 3 slip systems considered, which is the largest for the 2 oblique slip systems. Based on the deformed shape of the NML as shown in Figure 5, the slip direction that is inclined $\left({\theta }_1\right)$ w.r.t. to the layer parallel slip direction becomes the most active slip system based on the Schmid factors in the deformed lattice. However, it is not so clear what breaks the initial symmetry between the double-symmetric slip systems to make the average resolved shear stress on the ${\theta }_1$ oriented slip plane larger than the ${\theta }_2$ oriented plane. The absolute value of the averaged Schmid factor (over the whole domain) for the 3 slip systems in the initial configuration is 0 (layer-parallel slip vector) and 0.4665 (${\theta }_1$ and ${\theta }_2$ oriented slip vectors), while the absolute values of the averaged Schmid factors over the domain at $|e_T|=3.21\%$ are $0.0021$ (layer parallel slip vector), $0.4739$ (${\theta }_1$ oriented slip vector), and $0.4731$ (${\theta }_2$ oriented slip vector). At $|e_T|=23\%$, the absolute values of the averaged Schmid factors for the 3 slip systems are $0.1407$ (layer-parallel slip vector), $0.42246$ (${\theta }_1$ oriented slip vector), and $0.3357$ (${\theta }_2$ oriented slip vector).

The rotation of the deformed slip vectors w.r.t. their respective initial orientations is shown in Figure 5 at $|e_T|=23\%$, and it can be deduced that near the kink band region, the rotation of slip vectors is relatively large. The average rotation for the 3 slip vectors is $9.{99}^{\circ }$ (layer parallel slip vector), $10.{27}^{\circ }$ (${\theta }_1$ oriented slip vector), and $9.{64}^{\circ }$ (${\theta }_2$ oriented slip vector). The negative curl of the inverse elastic distortion tensor is the GND density, and the rotation of lattice vectors is a direct function of the (inverse) elastic distortion tensor. As can be observed from Figure 5, there are large gradients in the rotation field of slip vectors near the kinked region, which results in relatively large magnitudes of GND density in the kink band (refer to Figure 3a).

4.2. Effect of Different Slip Systems, Ordering of Metallic Layers, and Slip Systems Orientations

We report various numerical studies with respect to changing the ordering of metallic layers, the effect of removal of slip vectors, and the slip vector orientations on the kink band formation.

• Flipping the order of Cu and Nb layers in the nano-laminate structure resulted in a change in the orientation of the band as shown in Figure 6a, due to a geometrical asymmetry in the initial configuration.
• Removal of ${\theta }_1$ inclined slip system, and just considering layer parallel slip system and ${\theta }_2$ inclined slip system, resulted in the orientation of the band as shown in Figure 6b.
• Removal of the ${\theta }_2$ oriented slip plane did not change the orientation of the kink band formed, and it is similar to as shown in Figure 3.
• Removal of the layer parallel slip system and just considering the double symmetric slip system with the same orientations, as before, (for which the stress-strain response is shown in Figure 7) did not change the orientation of the kink band formed, as shown in Figure 8.
• We consider different orientations of slip vectors (${\theta }_1={\theta }_2$) from ${20}^{\circ }$ to ${45}^{\circ }$ for Cu, and keep the slip vector orientation for Nb to be the same i.e., ${\theta }_1={\theta }_2={45}^{\circ }$. Kink bands form for slip vector orientations from ${20}^{\circ }$ to ${30}^{\circ }$, but for ${35}^{\circ }$ to ${45}^{\circ }$ kink bands do not form, and NMLs undergo more uniform compression.

The stress-strain response with just 2 slip systems is relatively stronger as compared to the 3-slip system case, as shown in Figure 7. As concerns, the ∼ flat experimental stress-strain response beyond yielding, the availability of more slip systems in realistic modeling of the NML would lead to more plastic straining, but it would be accompanied by more gain in strength as well. So, in general, it is difficult to say with certainty that the availability of more slip systems will result in less strain-hardening, but that is what we observe in our limited simulations reported here.

We further compare the field plots of GND density (${\rho }_g$) with/without the layer parallel slip system as shown in Figure 8 at $|e_T|=23\%$ strain. It is evident that the shape of the kink band and the GND density are similar in the two cases. This is a departure from the results reported in the experimental work by Zhang et al. [8], where it was reported that the layer-parallel slip system is responsible for the evolution of the shape of the kink band. We leave a more realistic investigation of the effects of layer parallel slip systems for future work, where we will be modeling the 3-D geometry of NML and considering all the slip systems in Cu (FCC) and Nb (BCC).

Another interesting observation is that the domain-averaged rotation, with respect to the initial orientation, decreases with an increase in the initial orientation angle of the slip vectors from the loading axis, in the range of initial orientations where kink bands form. This shows that in order to accommodate the same level of plastic deformation, slip vectors with smaller angular deviation from the loading axis rotate more on average during plastic deformation.

4.3. Effect of Layer Thickness Length Scale on the Formation of Kink Band

We investigate the effect of layer thickness on the formation of kink bands in NMLs. It was experimentally observed in the work by Nizolek et al. [4] that kink bands do not form and homogeneous deformation is observed for layer thickness greater than $250$ nm. We run our simulation with the same parameters as stated in Table 5, and take the layer thickness $l=500$ nm, and the NML is composed of a total of 12 layers. Hence, the width of the NML composite is $W=6\mathsf{\mu }$m, and keeping the same aspect ratio of 1:2, the length of the NML composite is $L=12\mathsf{\mu }$m. It can be deduced from the deformed shape of NML (shown in Figure 9b) that kink bands do not form with larger layer thickness, and it is getting compressed homogeneously without the formation of any kink bands.

4.4. Conventional Plasticity within Our Model

In this section, we investigate the response of NML in three cases. For case (i), we take the usual parameters ($k_0=20,l=0.01732,\mathit{V}\ne \mathbf{0}$), while for case (ii), the parameters are ($k_0=0,l=0,\mathit{V}\ne \mathbf{0}$), and for case (iii), the parameters taken are ($k_0=0,l=0,\mathit{V}=\mathbf{0}$). Physically, case (ii) corresponds to just Voce’s law-based hardening (without the effects of GND hardening) and no plastic strain rate due to mesoscopic core energy. For case (iii), in addition to case (ii), there is no plastic straining due to the motion of GNDs, and this case is closest to conventional plasticity within our unified crystal plasticity and dislocation mechanics framework. As expected, with no hardening from GNDs accounted for in case (ii) and case (iii), we obtain a weaker strain-hardening response, as compared to case (i). It was observed in numerical simulations that kink bands do not form for case (ii) and case (iii), as shown in Figure 10b, which shows the GNDs density field plot for case (ii) at $|e_T|=12\%$ true strain.

4.5. NML Subjected to Layer-Perpendicular Compression

In this section, the layer direction of NML is aligned perpendicular to the compression loading direction, and the stress-strain response for the layer perpendicular compression of NML is compared with the layer parallel compression as shown in Figure 11. As observed in the experiments [3], kink bands do not form for layer perpendicular compression and the deformation of the NML is shear-driven parallel to the layer direction. This results in material flow under shearing within the individual layers as shown in Figure 12. For the layer perpendicular micropillar compression of NMLs in [3], it was reported that diagonally inclined shear bands are formed, from the left or right boundary to the right/left boundary. We leave such studies for future work where we model the slender micropillar shape of the NML using a plane strain idealization with the undulating interfaces. Up to the rotation of the layers towards the loading axis, similar behavior is observed in the work [23].

4.6. Comparison of Strain Gradient Plasticity (SGP) Models and MFDM

• The numerical results that we have obtained are with continuous plastic flow across the interfaces, imposing continuity of certain components of plastic strain rate across the interface. Moreover, we see the accumulation of GND density across the interface (refer to Figure 3a), without having to do constrained plastic flow across the interface, which models the impenetrable interfaces in MFDM framework, as discussed in [25]. This is unlike the recent work by Zecevic et al. [5], where ‘micro-hard’ boundary conditions are implemented on the interfaces for their numerical study. From the current state-of-the-art of SGP models [5,40] there appears to be a high degree of indeterminacy in the nature of boundary/interface conditions to be imposed, with a significant impact on model predictions.
  Using ‘micro-hard’ boundary conditions within SGP models would make it difficult to reproduce drastically different scalings in the micropillar compression experiments for two different configurations (${90}^{\circ }$ and ${45}^{\circ }$ oriented metal thin film with respect to compression loading axis), as experimentally observed in [41]. In the work of Kuroda et al. [40], the plastic constraints on the metal-ceramic interface are relaxed beyond a certain level of plastic strain gradient on the interface/boundary, in order to obtain the same scalings as experimentally observed. However, our model is able to reproduce similar scalings as compared to experimental ones, without any ad-hoc modifications to the boundary conditions, as shown in [23]. Moreover, our model also predicts the formation of kink bands in the compression of NMLs, as shown in this current (first) simplified study, again with no special fitting of material parameters or changes to the structure of the theory (which, of course, includes the nature of boundary and interface conditions). The jump conditions of our model have also been shown to be successful in studies of texture evolution and recrystallization, among others [27,42,43].
• Our model is not able to reproduce size dependence on the initial yield strength, as observed in the micropillar compression experiments of single crystal pure Ni [44]; to our knowledge, SGP-based theories have not modeled these experiments either. Most SGP theories are able to predict a significant size effect at initial yield in the presence of boundary constraints or in the presence of inhomogeneous deformation, without being able to fundamentally distinguish boundary constraints arising from the kinematics of dislocation slip [16,23].
• SGP studies of micropillar confined thin film plasticity and kink banding [5,40,45], usually employ a $\sim 0$ value of work hardening to reproduce experimental behavior, raising the question of accurate representation of macroscopic behavior by the models. Moreover, the SGP studies of [5,46] employ several variants of the core energy function and show a sensitive and significant dependence on stress-strain response to such a choice. In our work, no such choices are required, and the work hardening rates we use in our strength evolution ensure, at least in the confined thin film problem and for modeling macroscopic response, that a physically appropriate mechanical response is obtained. For the present kink-band problem, the obtained stress-strain response (without any fitting) is unsatisfactory compared to the experiment (but no worse than the SGP result of [5]), and this is an issue that requires further work.

5. Conclusions

This work reports a preliminary study of the modeling of kink-banding in NML within the MFDM framework. We report reasonable success, which further strengthens the ability of the MFDM framework to model important and challenging problems in the field of modern plasticity. Our simulations of a minimalistic model recover important experimental results in the compression of NMLs such as the following:

• dependence of kink-band formation on the layer thickness.
• kink bands do not form when the layer direction is aligned perpendicular to the compression loading direction, and there is shear-driven deformation within the layers of NML.

At this time, we are unable to provide a detailed analysis of the location of the kink bands in our numerical simulations. Rough speculation in this regard is as follows: the bands cannot cross the top and bottom boundaries because of the applied velocity boundary conditions. The corners of the sample result in stress concentrations due to geometric and changes in velocity boundary conditions on intersecting external boundaries. These stress concentrations may be the source of the location of the bands with the orientation governed by the most active slip systems. Of course, there are subtle details of the effects of the GND evolution that definitely affect the localization of the band and the layer length scale effects. Such a detailed understanding of the localization of the kink band is a fundamental challenge for future work as is the modeling of the full 3-D geometry of NML samples and the crystallography of actual slip systems present in Cu and Nb.

The prediction of experimentally observed kink bands in NMLs [4,8] along with the micropillar confined thin film plasticity [41,47] form a discerning, critical set of key experimental tests of mesoscale plasticity theories, to be accomplished without ad-hoc model modifications on a case-by-case basis. The MFDM model has been shown to be in qualitative accord with these experimental tests, building confidence in its adequacy as a model for the intricate phenomena involved in mesoscale plasticity.