# Standard gradient plasticity: formulation and numerical integration

Duy-Hao Luu, Habibou Maitournam, Son Nguyen

► To cite this version:

Duy-Hao Luu, Habibou Maitournam, Son Nguyen. Standard gradient plasticity: formulation and numerical integration. 11e colloque national en calcul des structures, CSMA, May 2013, Giens, France. hal-01717037

HAL Id: hal-01717037

https://hal.science/hal-01717037

Submitted on 25 Feb 2018

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

Public Domain
# Standard gradient plasticity: formulation and numerical integration

D.H. LUU<fnref n="1" />, M.H. MAITOURNAM<fnref n="2" />, Q.S. NGUYEN<fnref n="3" />

**Résumé** — The paper is devoted to the numerical implementation of a strain gradient plasticity model proposed by Q.S Nguyen. The model combining isotropic and kinematic hardening uses gradient of accumulated plastic strain to reproduce micro-mechanical features related to size effect. Its formulation is presented within the general framework of the *classical continuum thermodynamics*. Its numerical integration meets difficulties by the presence of strain gradients in the local theory formulation, which leads to a Laplacian equation solved at global level and to boundary value problems governed by partial differential equations of higher order with non-standard boundary conditions. A methodology of a diffusion-like problem is used to overcome this difficulty. This work is numerically solved via UMATs constructed in common available FE codes (e.g. Cast3M).

**Mots clés** — Strain Gradient Plasticity, Material length scales, Size Effect, Numerical integration

## 1 Introduction

In recent years there has been an increasing interest in modeling of plastic deformations at small scales. Examples are found in microelectronic components, micro-electro-mechanical systems (MEMS) and thin film applications. A number of experiments have demonstrated that length scale effect, or also called size effect, is of importance for sufficiently small test specimens of various geometries and loading situations (Venkatraman and Bravman 1992 [12]; N.A. Fleck et al. 1994 [4]; Nix and Gao 1998 [9]; Stolken and Evans 1998 [11]). The association with size effect is usually gradient effect.

Classical plasticity theory predicts no size neither gradient effects in micro-materials. New models are needed to capture relevant effects in micromechanics. To our knowledge, all these models are of gradient or nonlocal class, proposed since some recent decades. The existing models used in micro-mechanics can be enumerated as, after the synthesis of G.Z. Voyiadjis [13] : discrete dislocation dynamics simulations, molecular dynamics simulations, crystal plasticity theories, and gradient plasticity theories such as the one proposed here, henceforward briefly called *Gradient Model* (see works of N.A. Fleck & J.W. Hutchinson 2001 [3], H.B. Muhlhaus, G.Z. Voyiadjis, M.E. Gurtin, C. Polizzotto, R. de Borst, P. Gudmundson, S. Andrieux, E. Lorentz, etc., and especially those by Q.S. Nguyen 2005, 2010, 2011, 2012 which will be used as a basis for the development in the present work).

Thanks to their simplicity in numerically implementation, the gradient plasticity theory, pioneered by Aifantis (1984, 1987) [1, 2], is much developed and applied in micro-mechanics. In the framework of this work, a *strain gradient plasticity (SGP) model* proposed by Q.S. Nguyen [7, 5], will be addressed, in both *mathematical modeling and numerical implementation*.

The gradient formulation may be considered a higher-order extension of the local plasticity theory. This SGP theory still lies within the general framework of the classical continuum thermodynamics. However, the introduction of strain gradients into the local theory formulations leads to a Laplacian equation needing to be solved at global level and to boundary value problems governed by partial differential equations of higher order with non-standard boundary conditions. A methodology of a diffusion-like problem is used to overcome this difficulty. Indeed, this work numerically solve the *Laplacian equation at global level* via UMATs built in available Finite Element codes (Cast3M).

For the Gradient Model of mechanical behavior and time-independent processes, two models are considered : Gradient-Independent Dissipation Potential (Model 1), and Gradient-Dependent Dissipation Potential (Model 2). The last one, the most difficult due to its indetermination of dissipation forces for

<footnote n="1">¹ PhD student, Laboratoire de Mcanique des Solides (LMS), Ecole Polytechnique, Palaiseau, luuduyhao@gmail.com</footnote>

<footnote n="2">² CNRS UMR 7649, Laboratoire de Mcanique des Solides (LMS), Ecole Polytechnique, Palaiseau, habibou@lms.polytechnique.fr</footnote>

<footnote n="3">³ CNRS UMR 7649, Laboratoire de Mcanique des Solides (LMS), Ecole Polytechnique, Palaiseau, son@lms.polytechnique.fr</footnote>
null associated fluxes, is overcome with Q.S. Nguyen [7, 5] proposal of Energy Regularization Method.

The Gradient Model of accumulated plasticity with isotropic and kinematic hardening is examined in both mathematical and numerical aspects. The problem of thin wire torsion is treated as an illustration.

In numerical aspect, an implicit algorithm corresponding to the deformation theory is used. The question of convergence was also considered. A Ramberg-Osgood uniaxial stress-strain curve causes non-convergence in some problems. A conversion to another smoother exponential curve is recommended. Using the later, the convergence seems to be always achieved and even rather quickly. Application to typical problems show a good agreement between numerical results and references.

# 2 Formulation of the gradient plasticity model : application to Gradient Model of accumulated plasticity with isotropic-kinematic hardening

## 2.1 State and complementary laws

In small transformation, a Gradient Model of accumulated plasticity $\gamma$ with isotropic-kinematic hardening consists of the internal variable $\phi = (\varepsilon^p, \gamma)$, which represents the plastic strain and the equivalent plastic strain, such that $\varepsilon_{kk}^p = 0$ (plastic incompressibility) and $\gamma \ge 0$ (cumulated plastic strain), and of the following energy and dissipation potentials (started with the case of linear hardening)

$$ \left\{ \begin{aligned} W &= W(\varepsilon, \varepsilon^p, \gamma, \nabla\gamma) = W_e + W_c + W_i + W_g \\ &= \frac{1}{2}(\varepsilon - \varepsilon^p) : \mathbb{L} : (\varepsilon - \varepsilon^p) + \frac{1}{2}\varepsilon^p : \mathbb{H} : \varepsilon^p + \frac{1}{2}J\gamma^2 + \frac{1}{2}G\nabla\gamma \cdot \nabla\gamma \end{aligned} \right. \quad (1a) $$

$$ D = D(\dot{\varepsilon}^p, \dot{\gamma}, \nabla\dot{\gamma}, \varepsilon^p, \gamma) \iff \text{gradient-dependent (model II)} \quad (1b) $$

and under the regularized form :

$$ \left\{ \begin{aligned} W^* &= W^*(\varepsilon, \varepsilon^p, \gamma, \beta) = W(\varepsilon, \varepsilon^p, \gamma, \nabla\gamma) + \frac{1}{2}r(\beta - \nabla\gamma)^2 \\ &= \frac{1}{2}(\varepsilon - \varepsilon^p) : \mathbb{L} : (\varepsilon - \varepsilon^p) + \frac{1}{2}\varepsilon^p : \mathbb{H} : \varepsilon^p + \frac{1}{2}J\gamma^2 + \frac{1}{2}G\nabla\gamma \cdot \nabla\gamma + \frac{1}{2}r(\beta - \nabla\gamma)^2 \end{aligned} \right. \quad (2a) $$

$$ D^* = D^*(\dot{\varepsilon}^p, \dot{\gamma}, \dot{\beta}, \varepsilon^p, \gamma) \iff \text{gradient-independent (model I)} \quad (2b) $$

State and complementary laws offer the non-dissipative forces in this case :

$$ \left\{ \begin{aligned} \mathbf{X}_e^p &= W_{,\varepsilon^p}^* = -\sigma + \mathbb{H} : \varepsilon^p & \mathbf{Y}_e^p &= W_{,\nabla\varepsilon^p}^* = 0 \\ \mathbf{X}_e^\gamma &= W_{,\gamma}^* = J\gamma & \mathbf{Y}_e^\gamma &= W_{,\nabla\gamma}^* = G\nabla\gamma - r(\beta - \nabla\gamma) \\ \mathbf{X}_e^\beta &= W_{,\beta}^* = r(\beta - \nabla\gamma) & \mathbf{Y}_e^\beta &= W_{,\nabla\beta}^* = 0 \end{aligned} \right. \quad (3) $$

and the dissipative forces :

$$ \left\{ \begin{aligned} \mathbf{X}_d^p &= \partial D_{\varepsilon^p}^* \neq 0 & \mathbf{Y}_d^p &= \partial D_{\nabla\varepsilon^p}^* = 0 \\ \mathbf{X}_d^\gamma &= \partial D_\gamma^* \neq 0 & \mathbf{Y}_d^\gamma &= \partial D_{\nabla\gamma}^* = 0 \quad (\text{actually substituted with } \mathbf{X}_d^\beta) \\ \mathbf{X}_d^\beta &= \partial D_\beta^* \neq 0 & \mathbf{Y}_d^\beta &= \partial D_{\nabla\beta}^* = 0 \end{aligned} \right. \quad (4) $$

Applying generalized local equilibrium equation for forces associated with $\varepsilon_p$, $\gamma$ and $\beta$ yields

$$ (\mathbf{X}_e^p - \nabla \cdot \mathbf{Y}_e^p) + (\mathbf{X}_d^p - \nabla \cdot \mathbf{Y}_d^p) = 0 \implies \mathbf{X}_e^p + \mathbf{X}_d^p = 0 \quad (5a) $$

$$ (\mathbf{X}_e^\gamma - \nabla \cdot \mathbf{Y}_e^\gamma) + \mathbf{X}_d^\gamma = 0 \quad (\text{i.e., it comes back model I}) \quad (5b) $$

$$ (\mathbf{X}_e^\beta - \nabla \cdot \mathbf{Y}_e^\beta) + (\mathbf{X}_d^\beta - \nabla \cdot \mathbf{Y}_d^\beta) = 0 \implies \mathbf{X}_e^\beta + \mathbf{X}_d^\beta = 0 \quad (5c) $$

and then the dissipative forces are expressed as :

$$ \mathbf{X}_d^p = -\mathbf{X}_e^p = \sigma - \mathbb{H} : \varepsilon^p \quad (6a) $$

$$ \mathbf{X}_d^\gamma = -\mathbf{X}_e^\gamma + \nabla \cdot \mathbf{Y}_e^\gamma = -J\gamma + (G+r)\nabla^2\gamma - r\nabla \cdot \beta \quad (6b) $$

$$ \mathbf{X}_d^\beta = -\mathbf{X}_e^\beta = r(\nabla\gamma - \beta) \quad (6c) $$

Mises-like plastic criterion of the following form is considered :

$$ f = f(\mathbf{X}_d^p, \mathbf{X}_d^\gamma, \mathbf{X}_d^\beta) = \|\mathbf{X}_d^p\| + \mathbf{X}_d^\gamma + \frac{1}{\ell} \|\mathbf{X}_d^\beta\| - k \le 0 \quad (7) $$
## 2.2 Evolution laws

$$
\left\{
\begin{aligned}
\dot{\varepsilon}^p &= \lambda \frac{\partial f(\mathbf{X}_d^p, \mathbf{X}_d^\gamma, \mathbf{X}_d^\beta)}{\partial \mathbf{X}_d^p} = \lambda \frac{\mathbf{X}_d^p}{\|\mathbf{X}_d^p\|} \\
\dot{\gamma} &= \lambda \frac{\partial f(\mathbf{X}_d^p, \mathbf{X}_d^\gamma, \mathbf{X}_d^\beta)}{\partial \mathbf{X}_d^\gamma} = \lambda \\
\dot{\beta} &= \lambda \frac{\partial f(\mathbf{X}_d^p, \mathbf{X}_d^\gamma, \mathbf{X}_d^\beta)}{\partial \mathbf{X}_d^\beta} = \frac{1}{\ell} \lambda \frac{\mathbf{X}_d^\beta}{|\mathbf{X}_d^\beta|}
\end{aligned}
\right.
\quad \text{where } \lambda \ge 0; f \le 0; \lambda f = 0 \qquad (8)
$$

The additional boundary condition in this case is :

$$ Y^{\gamma} \cdot n = 0 \iff (Y_{e}^{\gamma} + Y_{d}^{\gamma}) \cdot n = 0 \quad (\text{using (3) and (4)}) \implies ((G+r)\nabla\gamma - r\beta) \cdot n = 0 \quad (9) $$

Finally, for a solid V subjected to a classical loading path and a controlled displacement on a portion of its boundary, in small and quasi-static transformation from a given initial state $u_o, f_o$, under the constitutive insulation condition following Polizzotto, the governing equations are :

$$ \bullet \forall x \in V $$

$$ \mathbf{X}_d^p = \sigma - \mathbb{H} : \varepsilon^p $$

$$ \mathbf{X}_d^\gamma = -J\gamma + (G+r)\nabla^2\gamma - r\nabla\cdot\beta $$

$$ \mathbf{X}_d^\beta = r(\nabla\gamma - \beta), \quad (\text{used to replace the role of } \mathbf{Y}_d^\gamma \text{ in loading functions as well as flow laws of model-II, thus leading to model-I problem}) $$

$$ \begin{aligned} f &= f(\mathbf{X}_d^p, \mathbf{X}_d^\gamma, \mathbf{X}_d^\beta) = \|\mathbf{X}_d^p\| + \mathbf{X}_d^\gamma + \frac{1}{\ell} \|\mathbf{X}_d^\beta\| - k \le 0 \\ \dot{\varepsilon}^p &= \lambda \frac{\partial f(\mathbf{X}_d^p, \mathbf{X}_d^\gamma, \mathbf{X}_d^\beta)}{\partial \mathbf{X}_d^p} = \lambda \frac{\mathbf{X}_d^p}{\|\mathbf{X}_d^p\|} \end{aligned} \qquad (10) $$

$$ \dot{\gamma} = \lambda \frac{\partial f(\mathbf{X}_d^p, \mathbf{X}_d^\gamma, \mathbf{X}_d^\beta)}{\partial \mathbf{X}_d^\gamma} = \lambda \quad \text{where } \lambda \ge 0; f \le 0; \lambda f = 0 $$

$$ \dot{\beta} = \lambda \frac{\partial f(\mathbf{X}_d^p, \mathbf{X}_d^\gamma, \mathbf{X}_d^\beta)}{\partial \mathbf{X}_d^\beta} = \frac{1}{\ell} \lambda \frac{\mathbf{X}_d^\beta}{|\mathbf{X}_d^\beta|} $$

$$ \nabla \cdot \boldsymbol{\sigma} + f_{vu} = 0 $$

$$ \bullet \forall x \in \partial V $$

$$ \boldsymbol{\sigma} \cdot \boldsymbol{n} = f_{su} \quad \text{on } \partial V_{fu}, \quad u = u_d \quad \text{on } \partial V_u, \quad \partial V_{fu} \cap \partial V_u = \partial V $$

$$ ((G + r)\nabla\gamma - r\beta) \cdot n = 0 \quad \text{on } \partial V $$

# 3 Numerical implementation

This section is devoted to the numerical implementation of the proposed Gradient Model, in which the *deformation plasticity method (deformation theory)* is used. First, a linear hardening is examined and then the generalization for nonlinear hardening is treated.

## 3.1 Governing equations

1. von Mises-like plastic criteria (7) :

$$ f = f(\mathbf{X}_d^p, \mathbf{X}_d^\gamma, \mathbf{X}_d^\beta) = \|\mathbf{X}_d^p\| + \mathbf{X}_d^\gamma + \frac{1}{\ell} \|\mathbf{X}_d^\beta\| - k \le 0 \qquad (11a) $$

$$ \text{with : } \|\mathbf{X}_d^p\| = \sqrt{\frac{3}{2}(\boldsymbol{\xi} : \boldsymbol{\xi})} = \sqrt{\frac{3}{2}} \|\boldsymbol{\xi}\| $$

$$ \boldsymbol{\xi} = \text{dev}(\mathbf{X}_d^p) = \text{dev}(\boldsymbol{\sigma} - \mathbb{H}\boldsymbol{\varepsilon}^p) = \boldsymbol{s} - \mathbb{H}\boldsymbol{\varepsilon}^p \quad (\mathbb{H} : \text{4-order tensor}) $$

$$ \mathbf{X}_d^\gamma = -J\gamma + (G+r)\nabla^2\gamma - r\nabla\cdot\beta \qquad (11b) $$

$$ \mathbf{X}_d^\beta = r(\nabla\gamma - \beta) \qquad (11c) $$

$$ \text{thus } f = \sqrt{\frac{3}{2}} \|\boldsymbol{\xi}\| - J\gamma + (G+r)\nabla^2\gamma - r\nabla\cdot\boldsymbol{\beta} + \frac{1}{\ell} |r(\nabla\gamma - \boldsymbol{\beta})| - k \le 0 \qquad (11d) $$
2. Evolution laws (8), with the boundary conditions as in (10):

$$
\left\{
\begin{array}{ll}
\bullet \dot{\epsilon}^p = \lambda \frac{\partial f(\mathbf{X}_d^p, \mathbf{X}_d^\gamma, \mathbf{X}_d^\beta)}{\partial \mathbf{X}_d^\beta} = \lambda \sqrt{\frac{3}{2}} \frac{\xi}{\|\xi\|} & \\
\bullet \dot{\gamma} = \lambda \frac{\partial f(\mathbf{X}_d^p, \mathbf{X}_d^\gamma, \mathbf{X}_d^\beta)}{\partial \mathbf{X}_d^\gamma} = \lambda & \\
\bullet \dot{\beta} = \lambda \frac{\partial f(\mathbf{X}_d^p, \mathbf{X}_d^\gamma, \mathbf{X}_d^\beta)}{\partial \mathbf{X}_d^\beta} = \frac{1}{\ell} \lambda \frac{\mathbf{X}_d^\beta}{|\mathbf{X}_d^\beta|} & \\
\bullet \text{the boundary conditions as in (10)} & \text{where } \lambda \ge 0; f \le 0; \lambda f = 0 \quad (12)
\end{array}
\right.
$$

It is well known that the associated model of deformation plasticity (Hencky) is obtained from
the previous incremental equations by replacing $\dot{\epsilon}^p$, $\dot{\gamma}$ and $\dot{\beta}$ by $\epsilon^p$, $\gamma$ and $\beta$. Thus, the following
equations hold instead of (12):

$$
\bullet \quad \epsilon^p = \lambda \sqrt{\frac{3}{2}} \frac{\xi}{\|\xi\|} \qquad (13a)
$$

$$
\bullet \gamma = \lambda \qquad \text{where } \lambda \ge 0; f \le 0; \lambda f = 0 \quad (13b)
$$

$$
\bullet \beta = \frac{1}{\ell} \lambda \frac{\mathbf{X}_d^{\beta}}{|\mathbf{X}_d^{\beta}|} \qquad (13c)
$$

$$
\bullet \text{ the boundary conditions as in (10)} \tag{13d}
$$

The deformation theory in which a 1-step scheme is used (implicit Hencky algorithm) will be dealt
with in the next item.

**3.2 Algorithm (implicit scheme)**

The case of deformation plasticity is thus obtained from the incremental description by a 1-step increment from the initial state. However, the problem is solved by iterations in which the total forces are updated as $F^{tot} = F^{ext} + F^p(\epsilon^p)$ after the value of $\epsilon^p$ of the previous iteration; $F^{ext}$ are the external forces and $F^p(\epsilon^p)$ the plastic forces associated to $\epsilon^p$.

For each iteration, all quantities $u, \varepsilon, \gamma, \varepsilon^p, \sigma, \text{etc}$ are re-calculated from the initial state 0, using the updated total forces $F^{tot}$.

All the following quantities are understood corresponding to the last increment with regard to the
considered problem, therefore accompanied without any mechanical iteration index.

For each mechanical iteration, one has to :

1. Calculate $u, \varepsilon$ as elastic solution

2. Compute $\gamma$ and $\beta$ by solving a Laplacian equation of unknown $\gamma$ after an explicit schema with internal iterations between $\gamma$ and $\beta$:

- Calculate $\beta$:

From (13c), $\beta$ is collinear with $\mathbf{X}_d^\beta$; combining with (11c) and noting the positivity of $r$, $\beta$ is
also collinear with $\nabla\gamma$. Finally one gets:

$$
\beta = \frac{\gamma \nabla \gamma}{\ell |\nabla \gamma|} = \beta(\gamma) \qquad (14)
$$

- Calculate $\gamma$:

$$
\xi = s - \mathbb{H}\epsilon^p = 2\mu e - (\mathbb{H} + 2\mu)\epsilon^p \quad (15)
$$

where $s = \text{dev}(\boldsymbol{\sigma})$; $\boldsymbol{\sigma} = \lambda_L \text{tr}(\boldsymbol{\varepsilon} - \boldsymbol{\varepsilon}^p)\mathbf{1} + 2\mu(\boldsymbol{\varepsilon} - \boldsymbol{\varepsilon}^p)$; and $e = \text{dev}(\boldsymbol{\varepsilon})$; $\lambda_L$ and $\mu$ are Lamé coefficients.

$\mathbb{H}\boldsymbol{\varepsilon}^p = \mathbb{H}\boldsymbol{\varepsilon}^p$ in the case of a spatially isotropic material (H : hardening modulus).

From (13a), $\xi$ and $\epsilon^p$ are collinear. Combining with (15), $\xi$ and $e$ are collinear too. Therefore
the last allows to write :

$$
\|\xi\| = 2\mu \|e\| - (H + 2\mu) \|\epsilon^p\| = 2\mu \|e\| - (H + 2\mu) \lambda \sqrt{\frac{3}{2}} \quad (16)
$$
Thus (11d) follows, with the mind that $\gamma = \lambda$ :

$$ -\ell^2 f_1(1+f_2)\nabla^2\gamma + \left( \frac{J+1.5H+r/\ell^2}{2\mu} + \frac{3}{2} \right)\gamma = \sqrt{\frac{3}{2}} \|e\| - \frac{k+r(\nabla\cdot\beta+|\nabla\gamma|/\ell)}{2\mu} \quad (17) $$

Here $r=f_2G$ with $G=2\mu\ell^2 f_1$; Coefficients $f_1$ and $f_2$ are defined and used instead of $r$ and $G$ for more flexibility. For example, a particular case, the non-gradient model ($r=0$) with linear isotropic hardening, i.e. $J(\gamma) = J\gamma$, gives :

$$ \gamma = \frac{2\mu}{J+3\mu} \left\langle \sqrt{\frac{3}{2}} \|e\| - \frac{k}{2\mu} \right\rangle_+ \quad (18) $$

$\langle \cdot \rangle_+$ denotes the positive part.

Consequently, (17) in a more familiar form, is :

$$ \left\{ \begin{array}{l} -K\nabla^2\gamma + c\gamma = f_{imp} \\ K = \ell^2 f_1(1+f_2) \quad \text{and} \quad c = \left(\frac{J+1.5H}{2\mu} + \frac{3}{2}\right) \\ f_{imp} = \left\langle \sqrt{\frac{3}{2}} \|e\| - \frac{k}{2\mu} \right\rangle_+ - \frac{r(\nabla \cdot \beta + |\nabla \gamma|/\ell)}{2\mu} \\ \text{Boundary conditions as in (10)} \end{array} \right. \quad (19) $$

With a given deformation field, $f_{imp} = f(\gamma)$ is a function of $\gamma$ as its includes $\beta$ which is a nonlinear function with respect to $\gamma$ due to the expression (14).

The second of the boundary conditions in (10) in this case, is :

$$ [(G+r)\nabla\gamma - r\beta] \cdot n = 0 \implies K\nabla\gamma \cdot n = \frac{r\beta}{2\mu} \cdot n = \Phi_d; \quad \Phi_d \text{ defined for use later in (25)} \quad (20) $$

In the more general case, a Gradient Model with nonlinear isotropic and kinematic hardening, the corresponding energy potential parts $W_i$ and $W_c$ in (1a) are no longer so simple, they are only known under derivative forms ($W_i' = J(\gamma)$ and $W_c' = H(\varepsilon^p)$). Here $J(\gamma)$ and $H(\varepsilon^p)$ are nonlinear functions with respect to arguments $\gamma$ and $\varepsilon^p$, respectively.

For example, a function of isotropic hardening $J(\gamma)$ can be deduced from a given Ramberg-Osgood uniaxial tensile stress-strain curve which form, denoting $E$ as Young modulus, is :

$$ \frac{\varepsilon}{\varepsilon_0} = \frac{\sigma}{\sigma_0} + \left(\frac{\sigma}{\sigma_0}\right)^n, \quad n = 5 \sim 10, \quad \varepsilon_0, \sigma_0 \text{ are material parameters, } \varepsilon_0 = \frac{\sigma_0}{E} \quad (21) $$

Or, a popular form of $J(\gamma)$ is an exponential function with two parameters $J_\infty$ and $b$ :

$$ J(\gamma) = J_\infty(1 - e^{-b\gamma}) \quad (22) $$

In this work, the second form of $J(\gamma)$ is used because it offers an easy and rapid numerical convergence. The first one causes divergence in some cases (see details in the next section). Then (17) becomes :

$$ -\ell^2 f_1(1+f_2)\nabla^2\gamma + \left(\frac{r}{2\mu\ell^2} + \frac{3}{2}\right)\gamma = \sqrt{\frac{3}{2}} \|e\| - \frac{k+J(\gamma)+1.5H(\varepsilon^p)+r(\nabla\cdot\beta+|\nabla\gamma|/\ell)}{2\mu} $$

And (19) becomes :

$$ \left\{ \begin{array}{l} -K\nabla^2\gamma + c\gamma = f_{imp} \\ K = \ell^2 f_1(1+f_2); \quad c = \frac{3}{2} \\ f_{imp} = \left\langle \sqrt{\frac{3}{2}} \|e\| - \frac{k}{2\mu} \right\rangle_+ - \frac{J(\gamma)+1.5H(\varepsilon^p)+r(\nabla\cdot\beta+|\nabla\gamma|/\ell)}{2\mu} \\ \text{Boundary conditions still as in (10) which second one is (20) too} \end{array} \right. \quad (23) $$
Note that $\xi$ and $e$ are collinear and that $e$ is known, (13a) allows to express :

$$ \varepsilon^p = \lambda \sqrt{\frac{3}{2}} \frac{\xi}{\|\xi\|} = \gamma \sqrt{\frac{3}{2}} \frac{e}{\|e\|} = \varepsilon^p(\gamma) \quad (24) $$

Then $H(\varepsilon^p) = H(\varepsilon^p(\gamma))$, a nonlinear function with respect to $\gamma$.

Thus $f_{imp}$ is a nonlinear function of $\gamma$, includes the nonlinear terms $\nabla \cdot \beta(\gamma)$, $J(\gamma)$, $H(\varepsilon^p(\gamma))$, $|\nabla \gamma|$, not only $\nabla \beta(\gamma)$ as before. Once $f_{imp}$ is known (via the value of $\gamma$ know from the previous iteration), a familiar diffusion-like equation is found which can be classically solved in global way. Eq. 23 is thus solved by a Newton method, using *internal iterations*.

Finally, (19) or (23) are the constitutive equations which give $\gamma$ and $\varepsilon^p$ in terms of $\varepsilon$ (or $e$ rather) as the solution of a diffusion problem in which the right member $f_{imp}$ is assumed known.

Hence, an *explicit scheme* is needed where $f_{imp} = f(\gamma)$ of the current iteration is known via the value $\gamma$ from the previous internal iteration; and then the Laplacian equation (19) or (23) are solved as in diffusion problem.

The principal steps for the internal iterations (index "i" is assigned to the $i^{th}$ iteration) are:

* Initializing $\gamma_1 = 0$ for the 1st internal iteration.

* Calculating for the $i^{th}$ internal iteration the value of $f_{imp}^i = f(\gamma_{i-1})$.

* Calculating $\gamma^i$ by solving the laplacian equation (as thermal problem) with $f_{imp}^i$ known.

* Repeating the above steps (steps 2nd, 3rd) until the convergence of $\gamma$.

Finally after calculating $\gamma$ and then the other quantities such as $\beta$ (with the help of (14)), $\varepsilon^p, \sigma, \dots$, and $F^{tot} = F_d^{ext} + F^p(\varepsilon^p)$ is updated for the next global iteration.

3. Restart a new global iteration until the convergence.

## 3.3 Numerical implementation

By finite elements, the following Galerkin representation is considered:

$$ u(x) = u^m U_m(x); \gamma(x) = \gamma^n N_n(x) $$

The nodal values $u^m$ and $\gamma^i$ lead to column matrices $[u]$ and $[\gamma]$ while $\varepsilon, \varepsilon^p$ and $\sigma$ are fields defined numerically at Gauss points. For a given field of plastic strain $\varepsilon^p$, the equilibrium equations under a given loading and plastic strain lead to a global matrix equation for the displacement $u$.

$$ [K]\{U\} = F_d^{ext} + F^p(\varepsilon^p) $$

where $[K]$ denotes the elastic stiffness matrix.

For a given displacement $[u]$, the associated fields of $\varepsilon^p$ and $\sigma$ and the equivalent plastic $[\gamma]$ are determined by a diffusion-like problem since (19) or (23) depending on the case of linear or nonlinear hardening, respectively.

As in diffusion problem, matrix equations are:

$$ \left\{ \begin{array}{l} [[K] + [C]] \{\gamma\} = \{F\} \\ [K] = \int_V \{\nabla N_i\}^T K \{\nabla N_j\} dV \quad \text{and} \quad [C] = \int_V \{N_i\}^T \rho c \{N_j\} dV \\ \{F\} = \int_V \{N_j\} f_{imp} dV + \int_{\partial V_q} \{N_j\} \Phi_d d\Gamma \end{array} \right. \quad (25) $$

where K, c and $f_{imp}$ can be directly found at (19) or (23), $\Phi_d$ at (20).

This equation is classical and can be found in computer codes, cf. CAST3M for example. The presence of nonlinear term $H(\varepsilon^p(\gamma), J(\gamma), \beta(\gamma), |\nabla \gamma|)$ in the expression of $f_{imp}$ and $\Phi_d$ must be taken into
account by purely numerical internal iterations. The positive condition $γ ≥ 0$ leads finally to the following algorithm, in the general case of nonlinear hardening :

$$ 
\left\{ 
\begin{array}{ll}
[\mathbb{K}] + [\mathbb{C}]\{\gamma^*\} = \{\mathbb{F}\} & \\
\gamma = 0 & \text{if } \|e\| < \sqrt{\frac{2}{3}}\frac{k}{2\mu} \text{ or } \gamma^* < 0, \text{ after (23)} \\
\gamma = \gamma^* & \text{otherwise}
\end{array}
\right. 
$$

## 3.4 Numerical convergence

It is noteworthy to note the following features of the numerical implementation.

* In mechanical aspect : in the present work, the deformation plasticity method with the implicit algorithm is used.

* In purely numerical aspect : an explicit scheme, as presented in the page 6, is required to numerically solve the Laplacian equation in which the variable $γ$ also appears in the right side.

* Numerical convergence question : The difficulty in convergence when dealing with the nonlinear hardening problem is that, once a functional of uniaxial tensile stress-strain such as Ramberg-Osgood is used, it may sometimes cause non-convergence due to its very sharp slope around the zone close to the origin. In order to overcome this inconvenience, a smoother functional as an exponential type as presented in (22), is strongly recommended (see Fig.1).

<figure><img src="image_1.png" /><figcaption>Fig. 1 – Conversion of a "uniaxial tensile stress-strain curve" of Ramberg-Osgood to an exponential curve (cf. N.A. Fleck and J.W. Hutchinson 2001 [3])</figcaption></figure>

# 4 Applications : thin wire torsion

Some most typical examples in applying SGP models for mechanical problems at micro scale are : (1) shearing of a thin film layer sandwiched between two substrates, (2) thin wire torsion, and (3) expansion of a spherical micro void. These problems permit to draw conclusions about the effectiveness of Gradient Model as well as the role the material length parameter. Here, only thin wire torsion problem is treated.

A solid cylindrical wire of radius R is considered. When the material is isotropic and the wire is twisted monotonically, the total shear strain is a known function of position according to $ε_s = αr$ where $α$ is the twist per unit length. It is assumed that there is no constraint of the plastic flow at the surface of the wire. Simulation results are presented in (Fig.2).

# 5 Conclusion

In this work, a Standard Gradient Plasticity Model in the context of Generalized Standard Models, is used to capture some relevant effects which are not included yet in the classical mechanics. Both mathematical modeling and numerical implementation aspects are addressed. The attention is focussed on time-independent processes such as incremental plasticity and brittle damage. In particular, the theoretical difficulty concerning the case of gradient-dependent dissipation is underlined.
<figure><img src="image_2.png" /><figcaption>Fig. 2 – Thin wire torsion simulation : torque versus twist for a solid wire of radius R (cf. N.A. Fleck and J.W. Hutchinson 2001 [3])</figcaption></figure>

The gradient formulation is established with the introduction of a *material characteristic length scale* *l*, and may be considered a higher-order extension of the local plasticity theory. In addition to the inherent boundary conditions of classical plasticity, a secondary boundary condition is introduced. The presence of strain gradient in the local formulation leads to a Laplacian equation solved at global level and to non-standard boundary value problem with partial differential equations of higher order. A computational method based on diffusion-like problem is used to overcome these difficulties in usual finite element. A Gradient Model of accumulated plasticity with isotropic and kinematic hardening, is detailed in the computational aspect. The thin wire torsion problem is treated as an illustration.

The formulation simplicity and distinct computational advantages of the class of models make it
attractive for applications in micro-mechanical problems.

Références

[1] Aifantis E.C. On the microstructural origin of certain inelastic models. J. Eng. Mat. Tech., 106 :326-330, 1984

[2] Aifantis E.C. The physics of plastic deformations. Int. J. Plasticity 3, 211-247, 1987

[3] Fleck N.A. and Hutchinson J.W. A reformulation of strain gradient plasticity. Journal of the Mechanics and Physics of Solids, 49(10) :2245-2271, 2001.

[4] Fleck N.A., Muller G.M., Ashby M.F., and Hutchinson J.W. Strain gradient plasticity : Theory and experiment. Acta Metallurgica et Materialia, 42(2) :475-487, 1994.

[5] Nguyen Q.S. Standard gradient models and gradient plasticity. Unpublished manuscript : A revisit of Q.S. Nguyen 2005 and 2011.

[6] Nguyen Q.S. Stability and Nonlinear Solid Mechanics. Wiley, Chichester, 2000.

[7] Nguyen Q.S. Variational principles in the theory of gradient plasticity. Comptes Rendus Mécanique, 339(12) :743-750, 2011.

[8] Nguyen Q.S. and Andrieux S. The non-local generalized standard approach : a consistent gradient theory. Comptes Rendus Mécanique, 333(2) :139-145, 2005.

[9] Nix W.D. and Gao H. Indentation size effects in crystalline materials : A law for strain gradient plasticity. Journal of the Mechanics and Physics of Solids, 46(3) :411-425, 1998.

[10] Polizzotto C. Unified thermodynamic framework for nonlocal/gradient continuum theories. European Journal of Mechanics - A/Solids, 22(5) :651-668, 2003.

[11] Stolken J.S. and Evans A.G. A microbend test method for measuring the plasticity length scale. Acta Materialia, 46(14) :5109-5115, 1998.

[12] Venkatraman R. and Bravman J.C. Separation of film thickness and grain boundary strengthening effects in al thin films on Si. Journal of Materials Resistance, 7 :2040-2048, 1992.

[13] Voyiadjis G.Z., Pekmezi G., and Deliktas B. Nonlocal gradient-dependent modeling of plasticity with anisotropic hardening. International Journal of Plasticity, 26(9) :1335-1356, 2010.