The thermodynamically consistent framework accounting for the thermomechanical behavior of the microstructure is addressed using the finite-element implementation. In particular, two different classes of the strain gradient plasticity (SGP) theories are proposed: In the first theory, the dissipation potential is dependent on the gradient of the plastic strain, as a result, the nonrecoverable microstresses do not have a value of zero. In the second theory, the dissipation potential is independent of the gradient of the plastic strain, in which the nonrecoverable microstresses do not exist. Recently, Fleck et al. pointed out that the nonrecoverable microstresses always generate the stress jump phenomenon under the nonproportional loading condition. In this work, a one-dimensional finite-element solution for the proposed strain gradient plasticity model is developed for investigating the stress jump phenomenon. The proposed strain gradient plasticity model and the corresponding finite-element code are validated by comparing with the experimental data from the two sets of microscale thin film experiments. In both experimental validations, it is shown that the calculated numerical results of the proposed model are in good agreement with the experimental measurements. The stretch-passivation problems are then numerically solved for investigating the stress jump phenomenon under the nonproportional loading condition.

# Higher-Order Thermomechanical Gradient Plasticity Model With Energetic and Dissipative Components OPEN ACCESS

**George Z. Voyiadjis**

**Yooseob Song**

**Taehyo Park**

^{1}Corresponding author.

Contributed by the Materials Division of ASME for publication in the JOURNAL OF ENGINEERING MATERIALS AND TECHNOLOGY. Manuscript received May 29, 2016; final manuscript received October 6, 2016; published online February 7, 2017. Assoc. Editor: Xi Chen.

*J. Eng. Mater. Technol*139(2), 021006 (Feb 07, 2017) (12 pages) Paper No: MATS-16-1154; doi: 10.1115/1.4035293 History: Received May 29, 2016; Revised October 06, 2016

## Abstract

## Introduction

It is well known that the classical continuum plasticity theory cannot capture the size effect of the microstructure during the course of plastic deformation [1–3]. A number of theoretical, experimental, and numerical works based on the SGP theories [4–7] have been proposed to describe the size-dependent plastic behavior of the microstructure since the work of Aifantis [8].

At the microlevel approach of the material, it is inevitable to develop a constitutive model based on the thermodynamics in order to have the physical understanding of the material behavior [9]. In their work, it is pointed out that the recoverable microforces are totally defined by the specific free energy function, and the nonrecoverable microforces are determined to be orthogonal to the dissipation function. Based on this decomposition, in-depth studies for the thermodynamically consistent gradient plasticity have been investigated by Gurtin and Anand [5,10–13]. Two sources of hardening, kinematic hardening in terms of the energetic backstress from an increase in free energy due to the density of the geometrically necessary dislocations (GNDs) and dissipative hardening in terms of an increase in slip resistance, are introduced by Gurtin and Anand [12]. Gurtin and Anand [13] discussed the physical interpretation of the nonlocal terms in the flow rules of Ref. [14] based on the thermodynamic principles and concluded that the flow rule of Ref. [14] does not satisfy the fundamental thermodynamic requirement (the plastic dissipation must be positive) under certain nonproportional loading history. The same violation on the plastic dissipation is also pointed out by Gudmundson [15] that the formulation of Ref. [14] does not guarantee the non-negative dissipation. Hutchinson [16] modified the formulation of Ref. [14] to meet the thermodynamic requirement by partitioning the thermodynamic microforces into recoverable and nonrecoverable components. In addition, the specific phenomenon in the work of Refs. [12,13], and [15], where a finite stress jump occurs due to an infinite change in the strain loading, and its physical interpretation are investigated by Hutchinson et al. [16–18]. Hutchinson [16] and Fleck et al. [17,18] pointed out that the nonrecoverable microstress quantities are the main cause of this phenomenon, and it is shown in their work that the elastic stress jump disappears if the nonrecoverable microstress quantities are dropped by assuming the dissipation potential function to be independent of the gradient of the plastic strain.

There are a couple of important issues that need to be addressed in the SGP theory. One of them is the effect of free surface and interface. It has been proven through the experimental works of Refs. [7] and [19] that the effect of interface and free surface is non-negligible. It is well known that the internal interface may act as an obstacle to block the dislocations, which lead to the enhancement of plastic resistance, while the free surface may act as a source for the development of dislocations. Moreover, when it comes to the numerical implementation, it is essential to incorporate the proper modeling of the interface into the higher-order SGP formulation for the well-posedness of the governing equations. Another important issue that should be considered, especially in the development of the constitutive equations of the small-scale metallic materials, is the coupling effect of the thermal activation energy along with the interaction mechanism of dislocations in terms of the SGP theories. It is well-known that a combination of applied stress and thermal activation is required to allow the dislocations to move through the lattice by overcoming an energy barrier rather than their individual effect. In addition, Gurtin and Reddy [20] pointed out that a theory that differentiates between recoverable and nonrecoverable microstresses requires a framework that allows for thermal variation. In this sense, special attention is given to the present framework of the small-scale metallic materials in order to incorporate the effect of the thermal variation into the constitutive models.

Meanwhile, numerous studies based on the numerical implementation of the SGP theory, in particular the finite-element simulation, have been carried out to qualitatively reproduce the material behavior of the microstructure such as the size effects observed in experiments [21–27]. In the initial stages, the primary objective of the finite-element simulation of the SGP theory is to demonstrate the mesh independency of the solution during plastic deformations [28]. The finite-element simulation for the SGP theory is too complicated, and a number of material parameters should be determined without adequate information. In spite of these drawbacks, it is worthwhile to keep investigating the SGP theory using the finite-element method because of its apparent advantages, flexibility and the ease in the simulation. In this work, the finite-element solution for the SGP theory is developed and applied to the stretch-passivation problem in order to investigate the elastic stress jump phenomenon under the nonproportional loading history. For the validation of the proposed model, the comparison with two sets of experiments, which have been carried out by Xiang and Vlassak [19] and Han et al. [29], is also performed. In addition, the comparison of the material responses between two classes of the SGP theories with and without the nonrecoverable higher-order microstresses under the nonproportional loading is addressed.

The main objective of the present investigation is to propose the thermodynamically consistent framework accounting for the thermomechanical behavior of the microstructure on the fast transient time using the finite-element method and to apply to the stretch-passivation problem in order to investigate the elastic stress jump phenomenon. To the best of our knowledge, no such studies have been investigated yet. The derived model is able to be used on solving very complicated phenomena such as the dynamic localization of inelastic flow in adiabatic shear bands and the perforation of metal plates by the projectiles at various impact speeds.

The outline of this work is as follows: The enhanced gradient plasticity framework for the bulk is presented in Sec. 2 by using the temperature and its gradient as an additional state variable. In particular, the formulations for two classes of the SGP theories with the dissipation potential dependent on $\epsilon \u0307ij,kp$ and with the dissipation potential independent of $\epsilon \u0307ij,kp$ are proposed along the specific form of the free energy function. The generalized thermodynamically consistent framework for the interface is discussed in Sec. 3 in order to develop the continuum-based microscopic behavior of the interface. A one-dimensional finite-element solution scheme for the proposed SGP theory is developed in Sec. 4. The proposed model and the corresponding finite-element code are validated by comparing with two sets of microscale experiments, which have been investigated by Han et al. [29] and Xiang and Vlassak [19], in Sec. 5. The comparison of the responses on the stress–strain curve between two classes of the SGP theories under the nonproportional straining history is carried out in Sec. 6, and some important conclusions are summarized in Sec. 7.

## Thermodynamically Consistent Nonlocal Strain Gradient Plasticity Theory: Bulk

The principle of virtual power enhanced by involving the temperature and its gradient as an additional state variable is used here. In this sense, the internal virtual power $PINT$ is expressed in terms of the energy contribution in the arbitrary volume $V$ as follows:

where $\sigma ij$ is the Cauchy stress tensor, $Rij$ is the backstress conjugate to the plastic strain rate $\epsilon \u0307ijp$, $Qijk$ is the higher-order microstress conjugate to the gradient of the plastic strain rate $\epsilon \u0307ij,kp$, and $A$ and $Bi$ are the generalized stresses conjugate to the rate of the temperature $T \u0307$ and the gradient of the temperature rate $T \u0307,i$, respectively. $PINT$ is balanced with the external virtual power $PEXT$ expended by the tractions $ti$ conjugate to the macroscopic velocity $vi$, the microtractions $mij$ and $a$ conjugate to $\epsilon \u0307ijp$ and $T \u0307$, and the body force $bi$ conjugates to $vi$ as follows:

where $S$ is the external surface. By the statement of the principle of virtual power $PINT=PEXT$ and using the divergence theorem, the equations for balance of linear momentum and nonlocal microforce balance can be represented, respectively, for volume $V$, as follows:

where $\tau ij=\sigma ij\u2212\sigma kk\delta ij/3$ is the deviatoric part of the Cauchy stress tensor, and $\delta ij$ is the Kronecker delta. The equations for local surface traction conditions and nonlocal microtraction conditions for the external surface $S$ can be expressed, respectively, with the outward unit normal vector to $S$, $nk$, as follows:

The second law of thermodynamics (entropy production inequality) requires that the rate of the free energy increase is not greater that the rate at which the work is done. Based on this, the entropy production inequality for the bulk is written as follows:

where $\rho $ is the mass density, $s$ is the specific entropy, $e$ is the specific internal energy, and $qi$ is the heat flux. By substituting the time derivative of $e$ from the equation for the Helmholtz free energy $\Psi =e\u2212Ts$, one can obtain the nonlocal free energy inequality, or Clausius–Duhem inequality as it is often called, as follows:

Decomposition of the microstresses $Rij$, $Qijk$, and $A$ into recoverable and nonrecoverable components gives

By taking the time derivative of the Helmholtz free energy function, which is assumed to be a smooth function of $\epsilon ije,\epsilon ijp,\epsilon ij,kp,T$, and $T,i$, and substituting it into Eq. (10), one can obtain the recoverable part of the thermodynamic forces as follows:

It is assumed in this work that the temperature and its gradient affect both the stored and the dissipated energy in order to address the nonequilibrium heat transfer due to the smallness of the structure and short period of time.^{2}

In the present framework, the Helmholtz free energy is assumed to have the following form:

where $Eijkl$ is the elastic modulus tensor, $\alpha t$ is the coefficient of linear thermal expansion, $Tr$ is the reference temperature, $h$ and $r$$(0<r<1)$ are the hardening material constants, $\epsilon p=\epsilon ijp\epsilon ijp$ is the accumulated plastic strain, $Ty$ and $n$ are the thermal material constants that need to be calibrated by comparing with the experimental data, $G$ is the shear modulus in the case of isotropic linear elasticity, $\u2113R$ is the recoverable length scale that describes the feature of the short-range interaction of GNDs, $c\epsilon $ is the specific heat capacity, and $a$ is a material constant that characterizes the interaction between phonon and electron. In addition, the term $(1\u2212(T/Ty)n)$ in Eq. (15) accounts for the thermal activation mechanism for overcoming the local obstacles to dislocation motion.

The recoverable counter parts of the thermodynamic forces such as $\sigma ij$, $AR$, $Bi$, $RijR$, and $QijkR$ can be obtained by using Eq. (12) as follows:

In this section, the dissipation potential functions for the aforementioned two classes of the SGP theories are postulated, respectively. The first class is derived from the dissipation potential dependent on the gradient of the plastic strain, which leads to the nonzero nonrecoverable thermodynamic microstress $QijkNR\u22600$, while the other class is derived from the dissipation potential independent on the gradient of the plastic strain, which leads to $QijkNR=0$.

Coleman and Gurtin [31] pointed out that the dissipation potential function is composed of two parts: the mechanical part, which is dependent on the plastic strain and its gradient, and the thermal counterpart, which shows the purely thermal effect such as the heat conduction. In this sense, the functional form of the dissipation potential for the former class can be put forward as

where $\sigma y$ is a material constant accounting for the yield strength; $m1$ and $m2$ are the non-negative material constants for the rate sensitivity parameter, in which the limit $m\u21920$ corresponds to rate-independent material behavior; $p \u03070$ is a constant for the reference flow rate; $n1$, $n2$, and $Ty1$ and $Ty2$ are the thermal material constant that need to be calibrated by comparing with the experimental data; $\u2113NR$ is the nonrecoverable length scale that corresponds to the dissipative effects in terms of the gradient of the plastic strain rate; $b$ is the material constant accounting for the energy exchange between phonon and electron; and $h(T)$ is the thermal conductivity coefficient. The generalized nonrecoverable effective plastic strain measure $p \u0307$ is defined as a function of the plastic strain rate, the gradient of the plastic strain rate, and the nonrecoverable length scale as follows [24]:

By using the dissipation potential given in Eq. (21) along with Eq. (14) and the assumption $h(T)/T=h0=constant$, the nonrecoverable thermodynamic forces for the former class can be obtained as follows:

On the other hand, the functional form of the dissipation potential for the latter class can be postulated by setting $\u2113NR=0$ in Eq. (21) as follows:

Meanwhile, the recoverable microstress quantities $RijR$ can be further decomposed into $(RijR)IR$, which mimics dissipative behavior by describing irreversible loading process, and $(RijR)R$, which describes reversible loading process due to the energy carried by dislocations [20]. Thus, the effective dissipation for the framework presented here can be defined as follows (see Ref. [30] for details):

The evolution of the temperature field is governed by the law of conservation of energy (the first law of thermodynamics). The terms addressing heating as a result of the inelastic dissipation and thermomechanical coupling are involved for describing the evolution of the temperature field. The equation for the conservation of energy in the present work is put forward as follows [32]:

By considering the effective dissipation potential given in Eq. (32) along with the equations for the entropy production (the second law of thermodynamics) previously derived in Eqs. (9) and (10), the relationship for the evolution of the entropy, which describes the irreversible process, can be derived as follows [33]:

By using Eq. (17) for solving the rate of the entropy $s \u0307$, the evolution of the temperature field can be obtained as follows [33]:

By substituting the constitutive equations of the recoverable microstresses given in Eqs. (16)–(20) and the nonrecoverable microstresses for the theory with the dissipation potential dependent on $\epsilon \u0307ij,kp$ given in Eqs. (23)–(26) into Eq. (35) and defining three additional terms $teff=h0/\rho cv$, $tR=\u0251Tr/h0$, and $tIR=bTr/\rho cv$, the rate of temperature can be obtained as follows [33]:

where $teff$ is the effective thermal diffusivity in the phonon–electron; $tR$ and $tIR$ are the microstructural (reversible) and transient (irreversible) time scales, respectively; and $\u03d1$ is a material constant dependent on the material microstructure. From the definition of the time scales, it is readily shown that the positive time scales are able to guarantee ɑ $>0$ and $b>0$; therefore, the Helmholtz free energy is guaranteed to be concave with respect to the temperature $T$ and its gradient $T,i$. The rate of temperature for the theory with the dissipation potential independent on $\epsilon \u0307ij,kp$ can be readily obtained by setting $\u2113NR=0$ in Eq. (36) as follows:

## The Effect of the Microscopic Boundaries on the Small-Scale Materials: Interface

As it is mentioned in Sec. 1, the free surface and interface can strongly influence the mechanical properties and the behavior of the material, especially when it comes to the small-scale metallic materials. Also, the effect of the grain boundaries of the microstructure on the thermomechanical material behavior increases as the surface/volume ratio increases. Besides the aforementioned physical viewpoint, nonstandard higher-order boundary conditions are required at the grain boundary for the well-posedness of the governing equations in the higher-order SGP theory from the mathematical point of view. In this sense, it is essential to incorporate the proper modeling of the microscopic boundary conditions into the higher-order SGP formulation.

The principle of virtual power is used to determine the balance of the microforces for the interfaces as it is used for the bulk in the previous section. An interface (grain boundary) separating grains $G1$ and $G2$ is taken into account here and it is assumed that the displacement field is continuous, i.e., $uiG1=uiG2$, across the grain boundary (Fig. 1). As shown in this figure, a dislocation moving toward the grain boundary in grain $G1$ cannot pass through the grain boundary but it is trapped and accumulated at the grain boundary due to the misalignment of the grains $G1$ and $G2$ that are contiguous to each other. In this sense, the grain boundary acts as an obstacle to block the dislocation movement; therefore, the yield strength of the material increases as the number of grain boundaries increases. By considering the effect of the interfacial surface energy dependent on the plastic strain rate at each side of the plastically deforming phase, the internal power for the interface can be expressed as follows [15,35,36]:

where $SI$ is an arbitrary subsurface of the grain boundary with an infinitesimal thickness, and $MijIG1$ and $MijIG2$ are the interfacial microforces expending power over the interfacial plastic strain rates at the $SIG1$ and $SIG2$, respectively. The interfacial external virtual power $PEXTI$, which is balanced with the interfacial internal power $PINTI$, is expended by the macrotractions $\sigma ijG1(\u2212njI)$ and $\sigma ijG2(njI)$ conjugate to the macroscopic velocity $vi$, and the microtractions $QijkIG1(\u2212nkI)$ and $QijkIG2(nkI)$ conjugate to $\epsilon \u0307ijp\u2009IG1$ and $\epsilon \u0307ijp\u2009IG2$, respectively, as follows:

Equating $PINTI=PEXTI$ by considering the arbitrary variation of the plastic strain at the interface, the interfacial macro- and microforce balances can be obtained as follows [30]:

The same procedure, which was used for development of the constitutive equations for the single grain, is carried out again in this section. The interfacial Clausius–Duhem inequality can be retrieved with the interfacial microtraction $MijI$ and the interfacial Helmholtz free energy per unit surface area $\Psi I$ as follows:

By assuming $\Psi I=\Psi I(\epsilon ijpI)$ to be convex with respect to $\epsilon ijpI$, Eq. (43) can preserve the positive value of the dissipation as it is reduced to the following expression:

It is further assumed here that the interfacial microtraction $MijI$ is decomposed into the recoverable $MijI,R$ and nonrecoverable $MijI,NR$ components, such as

where the components $MijI,R$ and $MijI,NR$ represent the mechanisms for the pre- and post-slip transfer, and thus involve the plastic strain at the interface prior to the slip transfer $\epsilon ijp\u2009I(e)$ and the one after the slip transfer $\epsilon ijp\u2009I(p)$, respectively. The overall plastic strain at the interface can be obtained by the summation of both plastic strains, such as

The recoverable component of the interfacial microforces can be obtained by substituting Eq. (45) into Eq. (44) as follows:

The nonrecoverable component of the interfacial microforces can then be defined as

Thus, the thermodynamic requirement that the plastic dissipation must be non-negative is satisfied if $MijI,NR$ is determined from Eq. (49) with the interfacial dissipation potential $DI$ to be convex with respect to $\epsilon \u0307ijpI$.

It is well known that the interface plays a role as the barrier to plastic slip in the early stages of plastically deforming phase, while it acts as a source of the dislocation nucleation in the later stages. The energetic condition in the area around the interface is affected by the long-range internal stress fields associated with constrained plastic flow, which lead to the accumulated and pile-up of dislocations near the interface. Thus, the condition at the interface is determined by a surface energy that depends on the plastic strain state at the interface [22].

The interfacial Helmholtz free energy per unit surface area of the interface in the present work is put forward under the guidance of Ref. [22], such as^{3}

By substituting the interfacial Helmholtz free energy per unit surface given in Eq. (50) into Eq. (47), the interfacial recoverable microstresses $MijI,R$ can be obtained as follows:

As can be seen in Eq. (51), $MijI,R$ do not involve the plastic strain rate, which is related to the dislocation slip and the temperature since the interfacial recoverable microstresses are activated by the recoverable stored energy.

In addition to the stored energy due to dislocation pile-up, there are two main mechanisms affecting the energy dissipation during the dislocation movement in the grain boundary area. The first mechanism is related to an energy change in the grain boundary region. The macroscopic accumulated plastic strain at the interface can be connected to the microscopic deformation of the interface through the quadratic mean of the deformation gradient. Thus, the energy change after the onset of slip transmission to the adjacent grain is able to be approximately determined by a quadratic function of the deformation gradient at the microscale and hence the interfacial plastic strain at macroscale. The other mechanism introduces the energy involved in the deformation of the grain boundary. This energy is mainly due to the energy dissipation during the dislocation movement and can be taken as a linear function of the interfacial plastic strain.

The interfacial dissipation potential $DI$ in the current study is postulated by combining the above-mentioned mechanisms as follows:

where $\u2113NRI$ is the interfacial nonrecoverable length scale, $mI$ and $p \u03070I$ are the viscous related material parameters, $\sigma yI$ is a constant accounting for the interfacial yield stress at which the interface starts to deform plastically, $hI$ is an interfacial hardening parameter representing the slip transmission through the interface, $TyI$ is the scale-independent interfacial thermal parameter at the onset of yield, $nI$ is the interfacial thermal parameter, and $\epsilon pI(p)=\epsilon ijp\u2009I(p)\epsilon ijp\u2009I(p)$ and $\epsilon \u0307pI(p)=\epsilon \u0307ijpI\u2009(p)\epsilon \u0307ijpI\u2009(p)$ are defined, respectively, with the plastic strain at the interface after the slip transfer $\epsilon ijp\u2009I(p)$ and its rate $\epsilon \u0307ijp\u2009I(p)$. The rate-dependency and temperature-dependency of the interfacial dissipation energy are clearly shown in Eq. (52) through the terms $(\epsilon \u0307pI(p)/p \u03070I)mI$ and $(1\u2212TI/TyI)nI$, respectively.

The interfacial nonrecoverable microstresses $MijI,NR$ can be obtained by substituting Eq. (52) into Eq. (49) as follows:

By substituting Eqs. (51) and (53) into Eq. (45), one can obtain the interfacial microtraction $MijI$ as follows:

As can be seen in Eq. (54), a free surface, i.e., microfree boundary condition, at the grain boundary can be described by setting $\u2113RI=\u2113NRI=0$ and it is also possible to describe a surface passivation, i.e., microclamped boundary condition, by setting $\u2113RI\u2192\u221e$ and $\u2113NRI\u2192\u221e$.

## A One-Dimensional Finite-Element Implementation for the Proposed Strain Gradient Plasticity Theory

In this section, a one-dimensional version of finite-element implementation for the SGP theory proposed in Secs. 2 and 3 is presented because of the complexity of the proposed governing partial differential equations. The developed finite-element solution is then validated in Sec. 5 by comparing with two sets of microscale experiments, which have been conducted by Han et al. [29] and Xiang and Vlassak [19].

An initially uniform single grain with two grain boundaries is solved under macroscopically uniform uniaxial stress as shown in Fig. 2. The grain with the size of $L$ is assumed to be infinitely long along the *x*-direction and initially homogeneous; therefore, the solution depends only on the *x*-direction.

In a one-dimensional scheme, the macroscopic partial differential equations for balance of linear momentum equation (3) and the corresponding macroscopic boundary condition, i.e., either the prescribed spatially uniform stress $\sigma \u2020$ or the prescribed displacement $u\u2020$, are represented in the absence of body forces $bi$ as follows:

On the other hand, the microscopic partial differential equations for nonlocal force balances (Eq. (4)) and the corresponding microscopic boundary conditions at the grain boundaries of both sides are given as

The strong forms of the macroscopic and microscopic partial differential equations, i.e., Eqs. (55) and (58), can be represented in weak forms after applying the macroscopic boundary conditions, i.e., Eqs. (56) and (57), and the microscopic boundary conditions, i.e., Eqs. (59) and (60), as follows:

where the arbitrary virtual fields $u\u0303$ and $\epsilon \u0303p$ are assumed to be kinematically admissible weighting functions in the sense that $u\u0303x=0=u\u0303x=L=0\u2009and\u2009\epsilon \u0303x=0p=\epsilon \u0303x=Lp=0$.

The “User-Element” subroutine UEL in the commercial finite-element package abaqus/standard [37] is developed in the present work in order to numerically solve the weak forms of the macroscopic and microscopic force balances, Eqs. (61) and (62), respectively. In the finite-element solution of the current work, the grain with the size of $L$ is approximately discretized using the finite elements with the number of $el$, i.e., $L=\u2211e=1el\Omega e$. The displacement field $u(x)$ and the plastic strain field $\epsilon p(x)$ are discretized independently and both of the fields are taken as fundamental unknown nodal degrees-of-freedom. In this regard, the displacement field and the corresponding strain field $\epsilon (x)$ and the plastic strain field and the corresponding plastic strain gradient field $\epsilon ,xp(x)$ can be obtained by using the interpolation functions as follows:

where $\mathbb{N}u\xi $ and $\mathbb{N}\epsilon p\xi $ are the shape functions, and $Uu\xi $ and $E\epsilon p\xi $ are the nodal values of the displacements and the plastic strains at node $\xi $, respectively. The terms $nu$ and $n\epsilon p$ represent the number of nodes per single element for the displacement and the plastic strain, respectively.^{4}

Meanwhile, in the present work, the implicit Backward-Euler time integration scheme is adopted to preserve the adiabatic condition with respect to the sharp increase in temperature in the system. The short total time with a sizable number of small time steps in the course of the simulation is taken into account in order to introduce the effect of high strain rate loading.

Based on the previously-mentioned spatial and time discretization, the element-level system of equations for the nodal residuals for the displacement and the plastic strain of each finite-element $\Omega e$ are given by substituting Eqs. (63) and (64) into Eqs. (61) and (62) as follows:

where $(ru)\xi $ and $(r\epsilon p)\xi $ are the contribution of each element to the element-level residuals for the displacement and plastic strain, respectively, and are required to be defined in the UEL for every element during the whole of iteration. The last term in Eq. (66), i.e., $MI\mathbb{N}\epsilon p\xi $, is applied only for the nodes on the interface which is at $x=0$ and $x=L$ in the present work.

The global coupled system of equations, $(ru)\xi =0$ and $(r\epsilon p)\xi =0$, are solved by abaqus/standard [37] based on the Newton–Raphson iterative scheme. Occasionally, the modified Newton–Raphson method, referred to as quasi Newton–Raphson method, is employed in the case that the numerical solution suffers a divergence during the initial increment immediately after an abrupt change in loading. In quasi Newton–Raphson method, a specific correction factor, which is less than one, is multiplied by one portion of the stiffness matrix. By using this method, a divergence problem can be overcome, however, convergence is expected to be slow because of the intense computational cost.

The Taylor expansion of the residuals with regard to the current nodal values can be expressed by assuming the nodal displacement and the plastic strain in iteration $\zeta $ as $Uu\zeta $ and $E\epsilon p\zeta $, respectively, as follows:

where $\Delta Uu\xi =(Uu\zeta +1)\xi \u2212(Uu\zeta )\xi $, $\Delta E\epsilon p\xi =(E\epsilon p\zeta +1)\xi \u2212(E\epsilon p\zeta )\xi $, and $O((\Delta Uu\xi )2,(\Delta E\epsilon p\xi )2)$ is the big *O* notation to represent the terms of higher-order than the second degree. The residual is ordinarily calculated at the end of each time step, and the values of the nodal displacements and the plastic strains are updated during the iterations. The increments in nodal displacements and the plastic strains can be computed by solving the system of linear discretized equations shown in Eq. (69) with the Newton–Raphson iterative method

where $K\Omega e$ is the Jacobian (stiffness) matrix that represents the contribution of each element to the element-level residuals for the displacement and plastic strain. This matrix is required to be defined in the user-subroutine for every element during all the iterations.

From Eqs. (67) and (68) along with the discretization for the displacements (Eq. (63)) and the plastic strains (Eq. (64)) at the end of a time step and the functional forms of the recoverable and nonrecoverable higher-order stresses defined in Secs. 2.3, 2.4, and 3.2, the Jacobian matrix for the two classes of the SGP theory can be obtained, respectively, as follows:

and

## Experimental Validation of the Proposed Strain Gradient Plasticity Model

Thin films have been widely used, e.g., to provide insulating layers between conductors, diffusion barriers, hardness coating with the aim of protecting from scratch, and passivation across a wide range of technology-based industries. Moreover, thin films are highly important in the material science because thin films play a critical role in the investigation and development of the new materials. In this sense, the proposed strain gradient plasticity model and the corresponding finite-element code are validated, and the material parameters for the proposed model are also calibrated by comparing with the experimental observations from two sets of microscale thin film experiments in this section. In addition, the comparison between the proposed SGP model and Voyiadjis and Faghihi [30] is carried out to show the increase in accuracy of the proposed model. The first experiment, which was carried out by Han et al. [29], is the microtensile test on the nickel (Ni) thin films at elevated temperature. The experimental work of Xiang and Vlassak [19] on the size effect in the electroplated copper (Cu) thin films with the various microscale thickness is selected for the second validation of the proposed SGP model.

Han et al. [29] developed the microscale tensile testing system, which is composed of a high temperature furnace, a micromotor actuator, and the digital image correlation (DIC) system, for evaluating the mechanical properties of the Ni thin films at high temperatures. Dog bone shaped specimens used in their experiments were made by microelectromechanical system (MEMS) processes, and the primary dimensions of the specimen are shown in Fig. 3. The calibrated material parameters as well as the general parameters for Ni are presented in Table 1. The results of microscale tensile tests at four different temperatures, i.e., $25\u2009\xb0C$, $75\u2009\xb0C$, $145\u2009\xb0C$, and $218\u2009\xb0C$, and the corresponding numerical results from the proposed model are shown in Fig. 4. As shown in this figure, it is clear from both the experimental and numerical results that the Young's modulus is not affected by variations in temperature, while the yield and tensile strength decrease as the specimen temperature increases. In addition, Fig. 4 clearly shows that the Bauschinger effect is not affected very much by variations in the specimen temperature. Meanwhile, the calculated results of the proposed model compare better to the experimental data than that of Voyiadjis and Faghihi [30] (Fig. 4).

Xiang and Vlassak [19] investigated the size effect with a variety of film thicknesses on the plastic behavior of the freestanding electroplated Cu thin films by performing the plane-strain bulge test. In their experimental technique, the film of interest is deposited on a silicon (Si) frame, and long rectangular thin films are fabricated using standard micromachining technology. In this plane-strain bulge test, the rectangular freestanding membranes surrounded by a rigid Si frame are deformed in plane strain by applying a uniform pressure to one side of the membrane as shown in Fig. 5. The stress–strain response of either a single material membrane or a stack of multiple membranes adhered on a rigid Si frame can be obtained by using the following two simple equations:

where $P$ is the applied stress, $\Delta $ is the corresponding film deflection, $t$ is the film thickness, $2d$ is the film window width, and $\epsilon r$ is a residual strain in the film.

As can be seen in the work of Xiang and Vlassak [19], the stress–strain curves of the Cu thin films with a passivation layer on both surfaces clearly show the size effect due to the presence of a boundary layer with high dislocation density near the film–passivation layer interfaces. In this sense, the bulge test of electroplated Cu films with both surfaces passivated by $20nm$ of titanium (Ti) is considered here for the experimental validation of the proposed SGP model. In order to describe the passivation effect on the proposed SGP model, the microclamped condition, which causes the dislocations to be completely blocked at the grain boundary, is imposed at both surfaces of the Cu films. Meanwhile, the experiments are performed with the various thicknesses of the Cu films of $1.0\mu m$, $1.9\mu m$, and $4.2\mu m$. The average grain sizes in all the cases are given by $1.5\xb10.05\mu m$, $1.51\xb10.04\mu m$, and $1.5\xb10.05\mu m$, respectively, which are almost identical to each other. The displacement and pressure resolutions for this bulge tests system are $0.3\mu m$ and $0.1kPa$, respectively.

The calibrated and general material parameters for Cu are presented in Table 1, and the comparison between the experimental measurements from the bulge tests and the calculated results from the proposed SGP model is shown in Fig. 6. As it is clearly shown in this figure, the size effects according to the variation of the Cu film thicknesses are well observed in both the proposed SGP model and the experimental work of Xiang and Vlassak [19]. Moreover, the numerical results from the proposed model are in good agreement with the experimental measurements.

## Nonproportional Loading

Hutchinson [16] classified the strain gradient version of J_{2} flow theories into two classes: incremental theory developed by Fleck and Hutchinson and nonincremental theory developed by Gudmundson et al. (cf. see Refs. [12,13], and [15–18] for details). The specific phenomenon in the nonincremental theory that exhibits a finite stress jump due to infinitesimal changes in plastic strain that may occur under the nonproportional loading is noted and its physical acceptance is also discussed in the work of Refs. [17] and [18]. Hutchinson [16] concluded that discontinuous changes with infinitesimal changes in boundary loads are physically suspect. Fleck et al. [17,18] showed this phenomenon with two plane strain problems, stretch-passivation problem and stretch-bending problem, for nonproportional loading condition. In the work of Refs. [17] and [18], it is noted that the nonrecoverable higher-order microstress quantities $QijkNR$ always generate the stress jump for nonproportional loading problems.

In Sec. 2, the two classes of the higher-order strain gradient plasticity theories with and without the nonrecoverable microstresses $QijkNR$ were proposed in order to investigate the size-dependent material behavior in small-scale metallic volumes based on a thermomechanical version of the second thermodynamics law and a system of microscopic force balances. In this section, the stretch-surface passivation problem is taken into account for the nonproportional loading problem since the difference in predictions of the two theories is clearly revealed under the nonproportional loading condition, while the predictions of the two theories differ only slightly under the proportional loading condition [17].

In the one-dimensional stretch-passivation problem, the grain is deformed into the plastic regime by uniaxial tensile stretch with no constraint on plastic flow at the grain boundaries, i.e., microfree boundary conditions are imposed. At a certain point, the plastic flow is then constrained by blocking off the dislocations from passing out of the grain boundary, which leads the further plastic strain not to occur at the grain boundary, i.e., microclamped boundary conditions are imposed.

In this section, the finite-element implementation validated in Sec. 5.1 by the experimental work of Han et al. [29] and the corresponding material parameters presented in Table 1 are used again for introducing the finite stress jump phenomenon due to infinitesimal changes in the plastic strain. Further evaluation of this phenomenon is also provided here in order to compare the two proposed strain gradient plasticity models with and without the nonrecoverable higher-order microstresses under the nonproportional loading condition.

First, the microtensile test on Ni thin films at elevated temperatures solved in Sec. 5.1 is revisited to investigate the occurrence of the aforementioned stress jump phenomenon under the abrupt surface passivation. Figure 7 shows the numerical results on the stress–strain behaviors of thin films for the model with the nonrecoverable microstresses. The numerical simulations are performed with $T=25\u2009\xb0C$, $75\u2009\xb0C$, $145\u2009\xb0C$, and $218\u2009\xb0C$ in common with the cases in the experimental validation (Sec. 5.1). As shown in this figure, the magnitudes of the stress jump are less than expected in all the cases since the nonrecoverable length scale $\u2113NR$ is set $0.1$ which is much smaller than the recoverable length scale $\u2113R=1.0$. Nevertheless, the very first slopes immediately after the passivation are calculated as $E25\xb0C=58.0GPa$, $E75\xb0C=59.2GPa$, $E145\xb0C=72.6GPa$, and $E218\xb0C=105.0GPa$, respectively, and this shows the responses immediately after the passivation get gradually closer to the elastic response $E=115GPa$ as the temperature increases.

The stress jump phenomenon is more clearly observed in Fig. 8. Moreover, the comparison of the results between the SGP model with the dissipative potential dependent on $\epsilon \u0307ij,kp$ and the one with the dissipative potential independent on $\epsilon \u0307ij,kp$ is addressed in Fig. 8. The results for the latter model are computed with three different values of recoverable length scales $\u2113R=0.1,0.2$, and $0.3$, and the material parameters used in this numerical work are presented in Table 2. The behavior of the SGP model with $QijkNR\u22600$ is in stark contrast with that of the SGP model with $QijkNR=0$ after the passivation point. A significant stress jump with the slope $EP$ similar to the modulus of elasticity $E$ is shown in the SGP model with $QijkNR\u22600$. On the other hand, no elastic stress jump is observed in the SGP model with $QijkNR=0$. This result is exactly in agreement with the prediction in Refs. [17] and [18]. In the case that the dissipative potential is independent of $\epsilon \u0307ij,kp$, the contribution from the gradients of the plastic strain is entirely energetic as can be seen in Sec. 2.4. Both the increase in the yield strength in the early stages of passivation and subsequent hardening due to the effects of the plastic strain gradient are observed along with the increase of the energetic length scale as shown in Fig. 8.

## Conclusions

In this work, the thermodynamically consistent framework accounting for the thermomechanical behavior of the microstructure on the fast transient time is addressed using the finite-element method. The enhanced gradient plasticity framework for the bulk is presented by employing the temperature and its gradient as an additional state variable. In particular, the formulations for two classes of the strain gradient plasticity theories are proposed along the specific form of the free energy function. The main distinction between the two theories is the presence of the nonrecoverable higher-order microstress quantities $QijkNR$. In the first theory, the dissipation potential is dependent on $\epsilon \u0307ij,kp$, as a result, $QijkNR$ do not have a value of zero in this formulation. In the second theory, the dissipation potential is independent on $\epsilon \u0307ij,kp$, in which, $QijkNR$ do not exist. Fleck et al. [17,18] noticed that $QijkNR$ always gives rise to the stress jump phenomenon, which causes a controversial dispute in the field of strain gradient plasticity theory with respect to whether it is physically acceptable or not, under the nonproportional loading condition. In addition, the generalized thermodynamically consistent framework for the interface is provided in this work in order to develop the continuum-based microscopic behavior of the interface.

A one-dimensional finite-element solution for the proposed SGP theory is developed using the commercial finite-element package abaqus/standard [37] via the user-subroutine UEL. In the developed code, the displacement field and the plastic strain field are discretized independently from each other, in other words, the plastic strain field is employed as an additional nodal degrees-of-freedom in addition to the displacement field.

The proposed strain gradient plasticity model and the corresponding finite-element code are then validated, and the material parameters for the proposed model are also calibrated by comparing with the experimental observations from the two sets of microscale thin film experiments in this work. In addition, the comparison between the proposed SGP model and Voyiadjis and Faghihi [30] is carried out to show the increase in accuracy of the proposed model. The first experiment, which was carried out by Han et al. [29], is the microtensile test on Ni thin films at elevated temperatures, i.e., $25\u2009\xb0C$, $75\u2009\xb0C$, $145\u2009\xb0C$, and $218\u2009\xb0C$. It is clearly revealed from both the experimental and numerical results that the Young's modulus and Bauschinger effect are not affected very much by variations in temperature, while the yield and tensile strength decrease as the specimen temperature increases. Meanwhile, the calculated results of the proposed model compare better to the experimental data than that of Voyiadjis and Faghihi [30]. For the second experimental validation, the experimental work of Xiang and Vlassak [19] is employed to compare the numerical results from the proposed model to the experimental observations in terms of the size-dependent behaviors in the electroplated Cu thin films based on the plane-strain bulge test technique. The experiments and finite-element implementations are performed with the various thicknesses of the Cu films of $1.0\mu m$, $1.9\mu m$, and $4.2\mu m$. The size effect according to the variation of the Cu film thicknesses is well observed in both the proposed SGP model and the experimental work of Xiang and Vlassak [19]. Moreover, the numerical results from the proposed model are in good agreement with the experimental measurements.

Based on the validated finite-element implementation, the stretch-passivation problems are numerically solved for investigating the stress jump phenomenon under the nonproportional loading condition. The microtensile test on Ni thin films at elevated temperatures solved for the experimental validation is revisited. The magnitudes of the stress jump are less than expected since the setting value for the nonrecoverable length scale $\u2113NR$ is too small. However, it is shown that the very first responses immediately after the passivation get gradually closer to the elastic response as temperature increases. Further study is carried out for comparison of the numerical results from the SGP model with the dissipative potential dependent on $\epsilon \u0307ij,kp$ to the one with the dissipative potential independent on $\epsilon \u0307ij,kp$. It is observed that a significant stress jump occurred in the SGP model with the dissipative potential dependent on $\epsilon \u0307ij,kp$, while it did not occur in the SGP model with the dissipative potential independent on $\epsilon \u0307ij,kp$.