A gradient plasticity model for the simulation of shear localization

Abstract

The numerical simulation of shear localization under high strain rates can be modeled by a system of four partial differential equations including conservation of momentum, conservation of energy, elastic, and inelastic constitutive relations. This article introduces the gradient terms of the equivalent plastic strain to the inelastic equation based on the implicit gradient theory of plasticity to preserve the ellipticity for the shear band modeling. The model is constructed by the mixed finite element formulation with B-bar to reduce shear locking effects considering displacement, stress, equivalent plastic strain, and temperature as the solution field and thereafter solving the entire nonlinear governing system simultaneously. The performance of the gradient plasticity model is verified by two benchmark shear band problems, and the obtained numerical results are tested with the high-rate experimental results.

Introduction

Shear band is a dynamic failure, common to be observed in metallic materials as a highly localized zone of intense plastic deformation when subjected to the high strain rates.^1,2 As the fact that shear band is considered as a material instability propagating at the time scale of microseconds and its width is on the order of tens of microns, ³ modeling shear band problems can be a hard core in these aspects. Molinari and Clifton ⁴ proposed an analytical model based on the thermal perturbation concerning the competition between strain rate hardening and thermal softening. During the dynamic loading, different mechanisms such as strain hardening, strain rate hardening, and thermal softening are competing which results in three distinct stages in the process of shear band formation.^5,6 For numerical modeling, plastic flow models derived from experimental results dependent on strain rate, strain rate hardening, and temperature are recommended to describe the formation of shear localization.^7,8

Mesh sensitivity is another problem for numerical modeling of shear localization which leads to different results with respect to the mesh size as the absence of convergence. Since the shear band formation is extremely fast, on the order of microseconds, ³ while heat conduction is a relatively slow process, the adiabatic condition is widely assumed in many models^7,9,10 with the terminology of “adiabatic shear bands.” It is accepted that there is not enough time for the heat to be conducted away ¹¹ when applied with high strain rates, whereas this assumption has shown to lead to mesh-dependent results.

From the physical standpoint and considering the effect of thermal diffusion, the problem is revealed by the inclusion of the intrinsic length scale created by the competition between heat conduction and heating.^12–14 Moreover, as explained in the works by McAuliffe and Waisman,^15,16 the diffusive term and the heat source are of the same order of magnitude. Strain gradient theories are other regularization techniques to alleviate the mesh sensitivity for strain softening materials.^17,18 Higher order gradients of strain were first introduced into the plastic constitutive equations by Aifantis ¹⁹ and Coleman and Hodgdon ²⁰ for rigid plastic materials accompanied with a gradient parameter called length scale. A phenomenological strain gradient plasticity model is proposed ²¹ and supported by the size effect in the torsion experiments of thin copper wires with micrometer diameters.^22,23 The length scale acts as a material property connecting the microstructure of the metallic material such as the grain size, dislocation density, and microstresses to the continuum^10,24,25 and shows the ability of the microstructure to transfer information with neighboring points in a certain distance. Shi et al. ²⁶ used the theory of mechanism-based strain gradient (MSG) plasticity to study plastic flow localization problems in ductile materials and determined the thickness analytically based on a bifurcation analysis. Similarly, Chen et al. ²⁷ studied the plastic flow localization in ductile materials using the couple stress (CS) strain gradient theory and the strain gradient hardening (softening) law (C-W) and showed that the strain rate was non-uniform within the band and reached a maximum at the center of the shear band. Voyiadjis et al. ²⁸ proposed a higher order thermodynamically strain gradient plasticity model which accounts for the thermomechanical behavior of the microstructure and examined the numerical results with experimental data.

In this work, four partial differential equations (PDEs) describing conservation of momentum, conservation of energy, elastic constitutive relation, and inelastic constitutive relation are used to model shear band problems. We first demonstrate the theory of gradient plasticity on the theoretical basis of nonlocal theories and then based on the implicit gradient theory of plasticity, and we introduce the gradient term of equivalent plastic strain (EQPS) to the PDE’s system in shear band modeling. Furthermore, we detail the weak form of the governing system and solve the entire nonlinear PDEs simultaneously using the mixed finite element ²⁹ taking displacement, stress, EQPS, and temperature as the solution field. Finally, we investigate the performance of the gradient plasticity model on a benchmark shear band problem using the mixed finite element formulation with B-bar to reduce the shear locking effects.

The article is organized as follows. Section “Problem statement and mechanical model” presents the gradient plasticity theory and the mechanical model. In this section, higher order gradients of the EQPS are considered in the inelastic constitutive relation which preserve the ellipticity for the simulation of shear band problems. ³⁰ The weak form of the governing equations is discussed as well in this section. In section “Numerical results,” a two-dimensional (2D) benchmark shear band problem where a shear band is formed at $45\circ$ and a pure shearing problem are selected to verify the performance of the gradient plasticity model. Moreover, the numerical results for $45\circ$ shear band are compared with the results obtained from the experiments conducted by Marchand and Duffy. ⁵ Final concluding remarks and perspectives for future work are given in section “Conclusion.”

Problem statement and mechanical model

Gradient plasticity theory

Shear band width is usually on the order of tens of micrometers which is a macroscopic phenomenon,^3,31 but its evolution lies in the material microstructure. When the deformation is homogeneous in space, its gradient terms are small and negligible. But when it comes to the highly localized plastic deformation or the observed plasticity phenomenon which endures a size effect, the gradient terms play a dominant role. ²²

The gradient plasticity theory introduces the higher order plasticity terms into the conventional governing differential equations together with the intrinsic length scales to characterize the microstructure behavior.

In the nonlocal regime, the nonlocal plastic strain ${\overline{\gamma }}_p^{\mathrm{nl}}$ around a certain material point $x_0$ is computed as a weighted average of the local plastic strain

$$\begin{array}{c}{\overline{\gamma }}_p^{\mathrm{nl}}(x_0)=\frac{1}{\mathrm{\Psi }(\xi )}\underset{\mathrm{\Omega }}{\int }\psi (\xi ){\overline{\gamma }}_p(x_0+\xi )d\mathrm{\Omega }\hfill \\ \mathrm{with}\mathrm{\Psi }(\xi )=\underset{\mathrm{\Omega }}{\int }\psi (\xi )d\mathrm{\Omega }\hfill \end{array}$$

(1)

where ${\overline{\gamma }}_p$ is the local plastic strain, $\xi$ is the infinitesimal vector denoted by $x-x_0$ , and $\psi (\xi )$ is the weight function, which is determined by $\xi$ and is assumed to be homogeneous and isotropic. It is usually defined as the Gauss distribution, $\psi (\xi )=\exp (-{\xi }^2/2l_c^2)$ .

With the assumption of isotropy and neglecting the fourth and higher order terms, the local plastic strain is expanded into a Taylor series

$${\overline{\gamma }}_p(x_0+\xi )\stackrel{.}{=}{\overline{\gamma }}_p(x_0)+\nabla {\overline{\gamma }}_p(x_0)\xi +\frac{{\nabla }^2{\overline{\gamma }}_p(x_0)}{2!}{\xi }^2$$

(2)

Plugging equation (2) into equation (1), the following relation is obtained

$${\overline{\gamma }}_p^{\mathrm{nl}}={\overline{\gamma }}_p+\frac{l_c^2}{2}{\nabla }^2{\overline{\gamma }}_p$$

(3)

where ${\nabla }^2$ denotes the Laplacian operator, and $l_c$ is the intrinsic length scale or gradient coefficient which is required for dimensional consistency.

Equation (3) inevitably leads to a fourth-order derivative of the displacement fields and a C¹-continuity requirement of displacement. ³² In order to avoid this problem, equation (3) is differentiated twice and then substituted into equation (3) again; also, neglecting fourth-order terms yields

$${\overline{\gamma }}_p={\overline{\gamma }}_p^{\mathrm{nl}}-\frac{l_c^2}{2}{\nabla }^2{\overline{\gamma }}_p^{\mathrm{nl}}$$

(4)

The aforementioned “implicit gradient theory of plasticity” formulation results in $C^0$ displacement and EQPS shape functions which are adopted in the inelastic constitutive relation in this gradient plasticity model.

Problem statement

The shear band problem can be described by a system of four PDEs including two balance equations, momentum and energy, elastic, and inelastic constitutive equations. The momentum equation which ignores the body force is written as follows

$$\rho \stackrel{·}{\underset{¯}{\nu }}=\nabla ·\underset{¯}{\underset{¯}{\sigma }}$$

(5)

where $\underset{¯}{\nu }$ is the velocity field, $\underset{¯}{\underset{¯}{\sigma }}$ is the stress tensor, and $\rho$ is the material density. Super imposed dot $(·)\mathrm{denotes}$ a derivative with respect to time $\mathrm{and}\nabla ·$ is the divergence operator.

The elastic constitutive relation in rate form can be given as follows

$$\underset{¯}{\underset{¯}{\stackrel{·}{\sigma }}}={\underset{¯}{\underset{¯}{\underset{¯}{C}}}}^{\mathrm{elas}}:(\underset{¯}{\underset{¯}{D}}-{\underset{¯}{\underset{¯}{D}}}^T-{\underset{¯}{\underset{¯}{D}}}^p)$$

(6)

where ${\underset{¯}{\underset{¯}{\underset{¯}{C}}}}^{\mathrm{elas}}$ is the elasticity modulus tensor. The rate of deformation $\underset{¯}{\underset{¯}{D}}$ has been additively decomposed into three components: the rate of elastic deformation ${\underset{¯}{\underset{¯}{D}}}^e$ , the rate of thermal deformation ${\underset{¯}{\underset{¯}{D}}}^T$ , and the rate of plastic deformation ${\underset{¯}{\underset{¯}{D}}}^p$ . Under small deformation assumption, the rate of deformation tensor is as follows

$$\underset{¯}{\underset{¯}{D}}=\underset{¯}{\underset{¯}{\stackrel{·}{\epsilon }}}=\frac{1}{2}[\nabla \underset{¯}{\nu }+{(\nabla \underset{¯}{\nu })}^T]={\nabla }^s\underset{¯}{\nu }$$

(7)

$\mathrm{where}\nabla$ and ${\nabla }^s$ are the gradient and the symmetric part of the gradient operators, respectively.

The rate of thermal deformation is defined as follows

$${\underset{¯}{\underset{¯}{D}}}^T=\alpha \stackrel{·}{T}\underset{¯}{\underset{¯}{I}}$$

(8)

where $\alpha$ is the thermal expansion coefficient, $T$ is the temperature, and $\underset{¯}{\underset{¯}{I}}$ is the second-order identity tensor.

The rate of plastic deformation is defined using a $J2$ plasticity model as follows

$${\underset{¯}{\underset{¯}{D}}}^p=\frac{3}{2}\frac{{\stackrel{·}{\overline{\gamma }}}_p}{\overline{\sigma }}\underset{¯}{\underset{¯}{S}}$$

(9)

where ${\overline{\gamma }}_p$ is the local EQPS and $\overline{\sigma }$ is the equivalent stress.

The local EQPS rate is expressed as follows

$${\stackrel{·}{\overline{\gamma }}}_p=\sqrt{\frac{2}{3}{\underset{¯}{\underset{¯}{D}}}^p:{\underset{¯}{\underset{¯}{D}}}^p}=\frac{g(\overline{\sigma },T,{\overline{\gamma }}_p)}{\overline{\sigma }}\sqrt{\frac{3}{2}\underset{¯}{\underset{¯}{S}}:\underset{¯}{\underset{¯}{S}}}=g(\overline{\sigma },T,{\overline{\gamma }}_p)$$

(10)

Based on the aforementioned gradient theory (see equation (4)), the inelastic constitutive relation is derived as follows

$${\stackrel{·}{\overline{\gamma }}}_p={\stackrel{·}{\overline{\gamma }}}_p^{\mathrm{nl}}-\frac{l_c^2}{2}{\nabla }^2{\stackrel{·}{\overline{\gamma }}}_p^{\mathrm{nl}}=g(\overline{\sigma },T,{\overline{\gamma }}_p)$$

(11)

where $g(\overline{\sigma },T,{\overline{\gamma }}_p)$ is a flow law and the deviatoric stress $\underset{¯}{\underset{¯}{S}}$ is defined as follows

$$\underset{¯}{\underset{¯}{S}}=\underset{¯}{\underset{¯}{\sigma }}-\frac{\underset{¯}{\underset{¯}{I}}:\underset{¯}{\underset{¯}{\sigma }}}{3}\underset{¯}{\underset{¯}{I}}$$

(12)

The effective stress is given as follows

$$\overline{\sigma }=\sqrt{\frac{3}{2}\underset{¯}{\underset{¯}{S}}:\underset{¯}{\underset{¯}{S}}}$$

(13)

The energy balance equation is written as follows

$$\rho c_p\stackrel{·}{T}=\chi \underset{¯}{\underset{¯}{\sigma }}:{\underset{¯}{\underset{¯}{D}}}^p$$

(14)

where $c_p$ and $\chi$ are the specific heat and Taylor–Quinney coefficient, respectively. The term, $\chi \underset{¯}{\underset{¯}{\sigma }}:{\underset{¯}{\underset{¯}{D}}}^p$ , in the right-hand side (RHS) of equation (16) assumes that part of the plastic work produces heat and the rest is stored in the metal’s lattice^?. $\chi$ represents the ratio of the rate of heat produced by plasticity and the rate of the plastic work.

Considering the following relation

$$\begin{array}{c}\multicolumn{1}{c}{\underset{¯}{\underset{¯}{\sigma }}:{\underset{¯}{\underset{¯}{D}}}^p=[\underset{¯}{\underset{¯}{S}}+\frac{\mathrm{tr}(\underset{¯}{\underset{¯}{\sigma }})}{3}\underset{¯}{\underset{¯}{I}}]:\frac{3}{2}\frac{g(\overline{\sigma },T,{\overline{\gamma }}_p)}{\overline{\sigma }}\underset{¯}{\underset{¯}{S}}}\\ =\frac{3}{2}\frac{g(\overline{\sigma },T,{\overline{\gamma }}_p)}{\overline{\sigma }}\underset{¯}{\underset{¯}{S}}:\underset{¯}{\underset{¯}{S}}=\overline{\sigma }g(\overline{\sigma },T,{\overline{\gamma }}_p)\hfill \end{array}$$

(15)

the energy equation can be simplified as follows

$$\rho c_p\stackrel{·}{T}=\chi \overline{\sigma }g(\overline{\sigma },T,{\overline{\gamma }}_p)$$

(16)

Finally, the four governing equations which describe four unknown fields of velocity, EQPS, stress, and temperature, over the problem domain $\mathrm{\Omega }$ are summarized as follows

$$\{\begin{array}{c}\begin{array}{c}\multicolumn{1}{c}{\rho \underset{¯}{\stackrel{·}{\nu }}=\nabla \underset{¯}{\underset{¯}{\sigma }}}\\ \multicolumn{1}{c}{{\stackrel{·}{\overline{\gamma }}}_p^{\mathrm{nl}}=g(\overline{\sigma },T,{\overline{\gamma }}_p)+\frac{l_c^2}{2}{\nabla }^2{\stackrel{·}{\overline{\gamma }}}_p^{\mathrm{nl}}}\\ \underset{¯}{\underset{¯}{\stackrel{·}{\sigma }}}={\underset{¯}{\underset{¯}{\underset{¯}{C}}}}^{\mathrm{elas}}:({\nabla }^s\underset{¯}{\nu }-\frac{3}{2}\frac{g(\overline{\sigma },T,{\overline{\gamma }}_p)}{\overline{\sigma }}\underset{¯}{\underset{¯}{S}}-\alpha \stackrel{·}{T}\underset{¯}{\underset{¯}{I}})\hfill \\ \multicolumn{1}{c}{\rho c_p\stackrel{·}{T}=\chi \overline{\sigma }g(\overline{\sigma },T,{\overline{\gamma }}_p)}\end{array}\hfill \end{array}$$

(17)

The boundary conditions are as follows

$$\{\begin{array}{cc}\underset{¯}{\nu }=\overline{v}\hfill & \mathrm{on}{\Gamma }^v\hfill \\ \multicolumn{1}{c}{T=\overline{T}}& \mathrm{on}{\Gamma }^T\hfill \\ \multicolumn{1}{c}{\underset{¯}{n}·\underset{¯}{\underset{¯}{\sigma }}=\overline{t}}& \mathrm{on}{\Gamma }^t\hfill \\ \multicolumn{1}{c}{\underset{¯}{n}·\nabla {\overline{\gamma }}_p^{\mathrm{nl}}=0}& \mathrm{on}{\Gamma }^{{\overline{\gamma }}_p}\hfill \end{array}$$

(18)

where $\underset{¯}{n}$ denotes the external normal unit vector and the initial conditions

$$\{\begin{array}{cc}v(\underset{¯}{x},t=0)\hfill & =v_0\hfill \\ \multicolumn{1}{c}{T(\underset{¯}{x},t=0)}& =T_0\hfill \\ \multicolumn{1}{c}{\sigma (\underset{¯}{x}t=0)}& ={\sigma }_0\hfill \\ \multicolumn{1}{c}{{\overline{\gamma }}_p(\underset{¯}{x},t=0)}& ={\overline{\gamma }}_{p,0}\hfill \end{array}$$

(19)

Note that when the gradient term of plastic strain is ignored, equation (17) comes to the adiabatic assumption of shear band simulation.

Discretization of the gradient model

The weak form of the system in equation (17) is obtained by multiplying the momentum equation, inelastic constitutive relation, energy equation, and elastic constitutive equation by the corresponding weight functions: $w_v$ , $w_{{\overline{\gamma }}_p}$ , $w_{\sigma }$ , and $w_T$ and then integrating over the problem domain $\mathrm{\Omega }$ . Hence, the weak form of the residual $\mathbf{R}=[{\mathbf{R}}_u,{\mathbf{R}}_{{\overline{\gamma }}_p},{\mathbf{R}}_{\sigma },{\mathbf{R}}_T{]}^T$ is as follows

$${\mathbf{R}}_u=\underset{\mathrm{\Omega }}{\int }(w_v\rho \underset{¯}{\stackrel{·}{\nu }}+\nabla w_v\underset{¯}{\underset{¯}{\sigma }})d\mathrm{\Omega }-\underset{{\Gamma }^t}{\int }{w}_v\overline{t}d{\Gamma }^t$$

(20)

$${\mathbf{R}}_{{\overline{\gamma }}_p}=\underset{\mathrm{\Omega }}{\int }{w}_{{\overline{\gamma }}_p}[{\stackrel{·}{\overline{\gamma }}}_p-g(\overline{\sigma },T,{\overline{\gamma }}_p)]+\frac{l_c^2}{2}\nabla w_{{\overline{\gamma }}_p}\nabla {\stackrel{·}{\overline{\gamma }}}_{p}d\mathrm{\Omega }$$

(21)

$$\begin{array}{c}{\mathbf{R}}_{\sigma }=\underset{\mathrm{\Omega }}{\int }{w}_{\sigma }[\underset{¯}{\underset{¯}{\stackrel{·}{\sigma }}}-{\underset{¯}{\underset{¯}{\underset{¯}{C}}}}^{\mathrm{elas}}:({\nabla }^s\underset{¯}{\nu }-\frac{3}{2}\frac{g(\overline{\sigma },T,{\overline{\gamma }}_p)}{\overline{\sigma }}\underset{¯}{\underset{¯}{S}}-\alpha \stackrel{·}{T}\underset{¯}{\underset{¯}{I}})]d\mathrm{\Omega }\hfill \end{array}$$

(22)

$${\mathbf{R}}_T=\underset{\mathrm{\Omega }}{\int }{w}_T(\rho c_p\stackrel{·}{T}-\chi \overline{\sigma }{\stackrel{·}{\overline{\gamma }}}_p)d\mathrm{\Omega }$$

(23)

In this model, the weak form contains spatial derivatives of velocity and EQPS; therefore, these two fields must be discretized with bilinear $C^0$ continuous functions, that is, $(\underset{¯}{\nu },{\overline{\gamma }}_p)\in C^0(\mathrm{\Omega })$ , and the stress and temperature fields are approximated by piecewise $C^{-1}$ functions, $(\underset{¯}{\underset{¯}{\sigma }},T)\in C^{-1}(\mathrm{\Omega })$ .

This model is discretized in time using the backward Euler method and the Galerkin approximation is assumed in each element. Hence, the linearization of the nonlinear system is required to solve the nonlinear problem at every time step. The Jacobian matrix with respect to the four variable fields is calculated analytically and defined by the Gâteaux derivative $\delta {\mathbf{x}}_{n+1}$ in the direction of each solution increment ³⁴

$$\begin{array}{c}\multicolumn{1}{c}{\mathbf{J}({\mathbf{x}}_{n+1})\delta {\mathbf{x}}_{n+1}={\frac{d\mathbf{R}({\mathbf{x}}_{n+1}+\epsilon \delta {\mathbf{x}}_{n+1})}{d\epsilon }|}_{\epsilon =0}}\\ =(\frac{\mathbf{M}}{\Delta t}+\mathbf{K})\delta {\mathbf{x}}_{n+1}\hfill \end{array}$$

(24)

where $\epsilon$ is an infinitesimal scalar and $\mathbf{K}$ is the consistent Jacobian matrix of ${\mathbf{F}}^{\mathrm{int}}-{\mathbf{F}}^{\mathrm{ext}}$ , and by taking the calculation of the first and second horizontal block as an example, the following equations are obtained

$$\begin{array}{c}\multicolumn{1}{c}{{\mathbf{J}}_{\mathrm{uu}}\delta \underset{¯}{u}={\frac{d{\mathbf{R}}_u(\underset{¯}{u}+\epsilon \delta \underset{¯}{u})}{d\epsilon }|}_{\epsilon =0}}\\ =\frac{d}{d\epsilon }{[\underset{\mathrm{\Omega }}{\int }({\omega }_u\rho \frac{{\partial }^2}{\partial t^2}(\underset{¯}{u}+\epsilon \delta \underset{¯}{u})+\nabla {\omega }_u·\underset{¯}{\underset{¯}{\sigma }})d\mathrm{\Omega }-\underset{{\Gamma }_u}{\int }{\omega }_u\overline{t}d{\Gamma }_t]|}_{\epsilon =0}\hfill \\ =\underset{\mathrm{\Omega }}{\int }{\omega }_u\rho \frac{{\partial }^2\delta \underset{¯}{u}}{\partial t^2}d\mathrm{\Omega }=\underset{\mathrm{\Omega }}{\int }{\omega }_u\rho \delta \underset{¯}{\stackrel{··}{u}}d\mathrm{\Omega }\hfill \end{array}$$

(25)

$${\mathbf{J}}_{u{\overline{\gamma }}_p}\delta {\overline{\gamma }}_p={\frac{d{\mathbf{R}}_u({\overline{\gamma }}_p+\epsilon \delta {\overline{\gamma }}_p)}{d\epsilon }|}_{\epsilon =0}=0$$

(26)

$$\begin{array}{c}\multicolumn{1}{c}{{\mathbf{J}}_{u\sigma }\delta \sigma ={\frac{d{\mathbf{R}}_u(\underset{¯}{\underset{¯}{\sigma }}+\epsilon \delta \underset{¯}{\underset{¯}{\sigma }})}{d\epsilon }|}_{\epsilon =0}}\\ ={\frac{d}{d\epsilon }[\underset{\mathrm{\Omega }}{\int }({\omega }_u\rho \frac{{\partial }^2\underset{¯}{u}}{\partial t^2}+\nabla {\omega }_u·(\underset{¯}{\underset{¯}{\sigma }}+\epsilon \delta \underset{¯}{\underset{¯}{\sigma }}))d\mathrm{\Omega }-\underset{{\Gamma }_u}{\int }{\omega }_u\overline{t}d{\Gamma }_t]|}_{\epsilon =0}\hfill \\ =\underset{\mathrm{\Omega }}{\int }\nabla {\omega }_u·\delta \underset{¯}{\underset{¯}{\sigma }}d\mathrm{\Omega }\hfill \end{array}$$

(27)

$${\mathbf{J}}_{\mathrm{uT}}\delta T={\frac{d{R}_u(T+\epsilon \delta T)}{d\epsilon }|}_{\epsilon =0}=0$$

(28)

$${\mathbf{J}}_{{\overline{\gamma }}_{p}u}\delta \underset{¯}{u}={\frac{d{\mathbf{R}}_{{\overline{\gamma }}_p}(\underset{¯}{u}+\epsilon \delta \underset{¯}{u})}{d\epsilon }|}_{\epsilon =0}=0$$

(29)

$$\begin{array}{c}\multicolumn{1}{c}{{\mathbf{J}}_{{\overline{\gamma }}_p{\overline{\gamma }}_p}\delta {\overline{\gamma }}_p={\frac{d{\mathbf{R}}_{{\overline{\gamma }}_p}({\overline{\gamma }}_p+\epsilon \delta {\overline{\gamma }}_p)}{d\epsilon }|}_{\epsilon =0}}\\ =\frac{d}{d\epsilon }{[\underset{\mathrm{\Omega }}{\int }{\omega }_{{\overline{\gamma }}_p}({\stackrel{·}{\overline{\gamma }}}_p+\epsilon \delta {\stackrel{·}{\overline{\gamma }}}_p)-{\omega }_{{\overline{\gamma }}_p}g(\overline{\sigma },T,{\overline{\gamma }}_p+\epsilon \delta {\overline{\gamma }}_p)+\frac{l_c^2}{2}\nabla {\omega }_{{\overline{\gamma }}_p}\nabla ({\stackrel{·}{\overline{\gamma }}}_p+\epsilon \delta {\stackrel{·}{\overline{\gamma }}}_p)]d\mathrm{\Omega }|}_{\epsilon =0}\hfill \\ =\underset{\mathrm{\Omega }}{\int }[{\omega }_{{\overline{\gamma }}_p}\delta {\stackrel{·}{\overline{\gamma }}}_p-{\omega }_{{\overline{\gamma }}_p}{\frac{d}{d\epsilon }g(\overline{\sigma },T,{\overline{\gamma }}_p+\epsilon \delta {\overline{\gamma }}_p)|}_{\epsilon =0}+\frac{l_c^2}{2}\nabla {\omega }_{{\overline{\gamma }}_p}\nabla \delta {\stackrel{·}{\overline{\gamma }}}_p]d\mathrm{\Omega }\hfill \end{array}$$

(30)

$$\begin{array}{c}\multicolumn{1}{c}{{\mathbf{J}}_{{\overline{\gamma }}_p\sigma }\delta \sigma ={\frac{d{\mathbf{R}}_{{\overline{\gamma }}_p}(\underset{¯}{\underset{¯}{\sigma }}+\epsilon \delta \underset{¯}{\underset{¯}{\sigma }})}{d\epsilon }|}_{\epsilon =0}}\\ ={\frac{d}{d\epsilon }\underset{\mathrm{\Omega }}{\int }-{\omega }_{{\overline{\gamma }}_p}g(\overline{\sigma }(\underset{¯}{\underset{¯}{\sigma }}+\epsilon \delta \underset{¯}{\underset{¯}{\sigma }}),T,{\overline{\gamma }}_p)d\mathrm{\Omega }|}_{\epsilon =0}\hfill \\ =\underset{\mathrm{\Omega }}{\int }-{\omega }_{{\overline{\gamma }}_p}\frac{d}{d\epsilon }{g(\overline{\sigma }(\underset{¯}{\underset{¯}{\sigma }}+\epsilon \delta \underset{¯}{\underset{¯}{\sigma }}),T,{\overline{\gamma }}_p)|}_{\epsilon =0}\hfill \end{array}$$

(31)

$$\begin{array}{c}\begin{array}{c}\multicolumn{1}{c}{{\mathbf{J}}_{{\overline{\gamma }}_{p}T}\delta T={\frac{d{\mathbf{R}}_{{\overline{\gamma }}_p}(T+\epsilon \delta T)}{d\epsilon }|}_{\epsilon =0}}\\ ={\frac{d}{d\epsilon }\underset{\mathrm{\Omega }}{\int }-{\omega }_{{\overline{\gamma }}_p}g(\overline{\sigma },T+\epsilon \delta T,{\overline{\gamma }}_p)d\mathrm{\Omega }|}_{\epsilon =0}\hfill \\ =\underset{\mathrm{\Omega }}{\int }-{\omega }_{{\overline{\gamma }}_p}\frac{d}{d\epsilon }{g(\overline{\sigma },T+\epsilon \delta T,{\overline{\gamma }}_p)|}_{\epsilon =0}\hfill \end{array}\hfill \end{array}$$

(32)

Here, the derivatives of the flow law with respect to ${\overline{\gamma }}_p$ , $\underset{¯}{\underset{¯}{\sigma }}$ , and $T$ , ${\mathrm{dg}(\overline{\sigma },T,{\overline{\gamma }}_p+\epsilon \delta {\overline{\gamma }}_p)/d\epsilon |}_{\epsilon =0}$ , ${\mathrm{dg}(\overline{\sigma }(\underset{¯}{\underset{¯}{\sigma }}+\epsilon \delta \underset{¯}{\underset{¯}{\sigma }}),T,{\overline{\gamma }}_p)/d\epsilon |}_{\epsilon =0}$ , and ${\mathrm{dg}(\overline{\sigma },T+\epsilon \delta T,{\overline{\gamma }}_p)/d\epsilon |}_{\epsilon =0}$ depend on the type of flow law, hence leading to the following linear system

$$\begin{array}{c}\mathbf{J}\delta \mathbf{x}=-\mathbf{R}\iff \hfill \\ \multicolumn{1}{c}{\left[\begin{array}{cccc}{\mathbf{M}}_u\hfill & 0\hfill & \Delta t{\mathbf{K}}_u\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& {\mathbf{M}}_{{\overline{\gamma }}_p}+\Delta t{\mathbf{G}}_{{\overline{\gamma }}_p}\hfill & \Delta t{\mathbf{G}}_{{\overline{\gamma }}_p\sigma }\hfill & \Delta t{\mathbf{G}}_{{\overline{\gamma }}_{p}T}\hfill \\ \multicolumn{1}{c}{\Delta t{\mathbf{K}}_{\sigma }}& \Delta t{\mathbf{G}}_{\sigma {\overline{\gamma }}_p}\hfill & {\mathbf{M}}_{\sigma }+\Delta t{\mathbf{G}}_{\sigma }\hfill & \Delta t{\mathbf{G}}_{\sigma T}+{\mathbf{M}}_{\sigma T}\hfill \\ \multicolumn{1}{c}{0}& {\mathbf{M}}_{T{\overline{\gamma }}_p}\hfill & \Delta t{\mathbf{G}}_{T\sigma }\hfill & {\mathbf{M}}_T+\Delta t{\mathbf{G}}_T\hfill \end{array}\right]}\\ \multicolumn{1}{c}{\times \left[\begin{array}{c}\delta \underset{¯}{u}\hfill \\ \multicolumn{1}{c}{\delta {\overline{\gamma }}_p}\\ \multicolumn{1}{c}{\delta \underset{¯}{\underset{¯}{\sigma }}}\\ \multicolumn{1}{c}{\delta T}\end{array}\right]=-\left[\begin{array}{c}{\mathbf{R}}_u\hfill \\ \multicolumn{1}{c}{{\mathbf{R}}_{{\overline{\gamma }}_p}}\\ \multicolumn{1}{c}{{\mathbf{R}}_{\sigma }}\\ \multicolumn{1}{c}{{\mathbf{R}}_T}\end{array}\right]}\end{array}$$

(33)

where $\mathbf{x}$ is the vector of unknowns defined as $\mathbf{x}=(u,{\overline{\gamma }}_p,\sigma ,T)$ . $\mathbf{G}$ represents the stiffness matrices produced from linearization of the inelastic constitutive law. Note that the Jacobian matrix is not symmetric.

Numerical results

In this section, we verify the role of the second-order gradient term of the EQPS in the inelastic constitutive equation which justified the PDEs for shear band problems being solved. The introduced intrinsic length scale in the model produces a material nonlocality acting as a localization limiter and hence regularizing the problem. We study the behavior of this gradient plasticity model on two shear band benchmark problems including a pure shear problem and a $45\circ$ problem where a shear band is formed at $45\circ$ due to high strain rate tension of a metallic plate and then compare the obtained numerical results with the experimental results from Marchand and Duffy. ⁵

The PDE model is implemented using the finite element analysis program (FEAP). ³⁵ The linearized systems of equations are solved with PETSc’s ³⁶ LU factorization algorithm using a quotient minimum degree ordering.^37,38 The postprocessing and graphics are generated using ParaView ³⁹ and Matplotlib. ⁴⁰

$45\circ$ shear band under uniform tension

In this section, a square plate is considered with uniform tension impact at its top and bottom edges. An imperfection is assumed at the center of the plate in order to generate X-shaped shear bands. Considering the symmetry to both the x-axis and the y-axis, only the top-right quarter of the plate is studied in the analysis. The geometry of the square plate with the edge dimension of $2H=100\mu \mathrm{m}$ , together with the profile of the applied velocity, $v_{\mathrm{BC}}$ , is given in Figure 1. The velocity loading increases linearly from 0 to 20 m/s in 1 μs and remains constant afterward which corresponds to a nominal strain rate of $2\times {10}^5{\mathrm{s}}^{-1}$ for this plate. Zero flux for temperature field is applied on all edges, and plane strain conditions are assumed.

Figure 1.

The geometry of the square plate under tension and the applied velocity profile. The imperfection is represented as a red dashed circle.

A hyperbolic secant type of imperfection is considered during the shear band simulation. Hence, the values of the yield stress and yield strain are reduced by ${\beta }_{2D}$ as follows

$${\sigma }_{\mathrm{yield}}={\sigma }_0{\beta }_{2D}(x,y)$$

(34)

$${\gamma }_{\mathrm{yield}}={\gamma }_0{\beta }_{2D}(x,y)$$

(35)

$${\beta }_{2D}(x,y)=1-\beta {\left[\mathrm{sech}\left(\frac{\sqrt{{(x-x_0)}^2+{(y-y_0)}^2}}{r_0}\right)\right]}^2$$

(36)

Here, $\beta =4\%$ is the reduction in the material parameters in percentage at the center of the imperfection; $x_0=0$ , $y_0=0$ represent the center of the imperfection with $r_0=1.5\times {10}^{-6}\mathrm{m}$ the radius in meters.

The material used for the analysis is HY-100 steel with a Johnson–Cook flow law ⁴¹ shown as follows

$$\overline{\sigma }=[A_s+B_s{({\gamma }_p)}^N][1-{\left(\frac{T-T_0}{T_m-T_0}\right)}^m][1+c\ln \left(\frac{{\stackrel{·}{\gamma }}_p}{{\stackrel{·}{\gamma }}_0}\right)]$$

(37)

where $\overline{\sigma }$ is the material flow stress. The material parameters associated with the Johnson–Cook parameters ⁴² for the material are detailed in Tables 1 and 2.

Table 1.

Material properties for HY-100 steel. ³⁸

Physical quantity name	Symbol	Value	Unit
Young’s modulus	$E$	207	GPa
Poisson’s ratio	$v$	0.3	−
Mass density	$\rho$	7746	kg m−3
Specific heat	$c_p$	502	m2 s−2 K−1
Thermal conductivity	$\kappa$	34	kg m s−3 K
Coefficient of thermal expansion	$\alpha$	$14e^{-6}$	K−1
Taylor–Quinney coefficient	$\chi$	0.9	−

Table 2.

Johnson–Cook parameters for HY-100 steel. ³⁸

Reference strain rate	${\stackrel{·}{\gamma }}_0$	1.0	s−1
Rate sensitivity parameter	$c$	0.011	−
Yield shear stress	$A_s$	$758e^6$	kg m−1 s−2
Shear stress hardening parameter	$B_s$	$402e^6$	kg m−1 s−2
Strain hardening exponent	$N$	0.26	−
Reference temperature	$T_0$	293	K
Reference melt temperature	$T_m$	1793	K
Thermal softening exponent	$m$	1.13	−

The top-right quarter plate is discretized into $50\times 50$ , $80\times 80$ , and $100\times 100$ structured quad elements. The total duration of the simulation is 3 μs. The natural material length scale $l=\mu b/S_0=0.55\mu \mathrm{m}$ , proposed by Gurtin, ⁴³ is used in our simulation where $\mu$ is the shear modulus, $b$ is the magnitude of Burgers vector, and $S_0$ is the slip resistance. The average stress and strain curve calculated within the shear band zone for the three different meshes are illustrated in Figure 2. The experimental curve from Marchand and Duffy ⁵ is also included for comparison. It is clear to see that the curves for different meshes are in close agreement indicating a mesh-independent analysis and a satisfactory level of convergence.

Figure 2.

The average stress–strain curves along the diagonal line of the quarter plate for three different meshes compared with the experimental results. ⁵

The EQPS and the temperature distribution along the diagonal line of the quarter plate (orange dashed line in Figure 1) for three different meshes at the final time are detailed in Figure 3.

Figure 3.

The equivalent plastic strain and temperature distribution along the diagonal line of the quarter plate plotted at the final time for mesh $50\times 50$ , $80\times 80$ , and $100\times 100$ .

These plots indicate that the EQPS and temperature within the shear band are much higher than that outside the band zone, and the width of the plasticity zone is not broadened with respect to the mesh size. Moreover, with the shear band simulation based on the irreducible mixed finite element formulation, the maximum local temperature variation is almost over 800 K. According to the experimental results conducted by Marchand and Duffy, ⁵ the surface temperature changes are at least 428°C (802.4°K) measured by a high-speed infrared radiometer. Hence, this gradient plasticity model captures the trend that is observed in the experiments subjected to high strain rates.

Plane strain shearing

The plane strain shearing problem is modeled in this section. The geometry of the plate is shown in Figure 4 with $H=1\times {10}^{-5}\mathrm{m}$ . The material parameters are the same as those used for the previous $45\circ$ shear band problem. No imperfection is introduced in this example, and the band is triggered by the stress concentration created by the boundary conditions. A horizontal velocity loading is applied on the upper left edge of the plate with $v_0=5\mathrm{m}/\mathrm{s}$ . The mesh for the analysis can be seen in Figure 5 with 11 structured elements (the smallest element size is $0.16\mu \mathrm{m}$ ) across the shear band zone corresponding to 4100 elements in total.

Figure 4.

The geometry of the shear plate with the impact loading. No imperfection is required.

Figure 5.

The mesh of the shear plate.

Figure 6 elaborates the EQPS distribution along the vertical line at the center of the plate plotted at three different times $t=1.0\mu \mathrm{s}$ , $t=1.6\mu \mathrm{s}$ , and $t=2.0\mu \mathrm{s}$ showing the evolution of shear band when enduring the impact loading. The snapshots of the EQPS field corresponding to different times are listed in Figure 7. The figures show that the width of the localized plasticity zone remains almost the same and does not widen during the shear localization of the plate.

Figure 6.

The EQPS evolution along vertical line at the center of the plate.

Figure 7.

Plots of the equivalent plastic strain at three different times for plane strain shearing problem.

Conclusion

This article developed an implicit gradient plasticity model considering the second-order gradient term of equivalent plastic strain in the PDE system for dynamic shear localization simulation. This model solves the entire PDEs simultaneously using the mixed finite element method with B-bar to reduce shear locking effects. Herein, the EQPS and temperature are interpreted as additional degree of freedoms and they are regarded as the unknown field together with displacement and stress. The residual of the nonlinear system is discretized with the Galerkin approximation spatially and with a backward Euler method in time. The Jacobian matrix is calculated analytically as the Gâteaux derivative of the residual avoiding the process of numerical differentiation. This implicit gradient plasticity model is tested by simulating two benchmark shear band problems. One is the $45\circ$ shear band problem where the X-shaped shear band is triggered by the central imperfection. Compared to the high-rate experiments, the numerical modeling using the gradient plasticity model captures the trends observed in the experiments. Another is the plane strain shearing problem where the shear band evolution is elaborated.