# A time-discontinuous elasto-plasticity formalism to simulate instantaneous plastic flow bursts

M Lamari, P Kerfriden, Oguz Umut Salman, V Yastrebov, K Ammar, S Forest

► To cite this version:

M Lamari, P Kerfriden, Oguz Umut Salman, V Yastrebov, K Ammar, et al.. A time-discontinuous elasto-plasticity formalism to simulate instantaneous plastic flow bursts. International Journal of Solids and Structures, 2024, 309, pp.113171. 10.1016/j.ijsolstr.2024.113171. hal-05390265

HAL Id: hal-05390265

https://hal.science/hal-05390265v1

Submitted on 1 Dec 2025

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.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# A time-discontinuous elasto-plasticity formalism to simulate instantaneous plastic flow bursts

M. Lamari<sup>a,*</sup>, P. Kerfriden<sup>a</sup>, O.U. Salman<sup>b</sup>, V. Yastrebov<sup>a</sup>, K. Ammar<sup>a</sup>, S. Forest<sup>a</sup>

<figure><img src="image_3.png" /></figure>

## ARTICLE INFO

Dataset link: https://zenodo.org/records/14266823, 10.5281/zenodo.14266823, https://github.com/Mathias-Lamari/Time-discontinuous-plasticity

**Keywords:**

Intermittent plasticity

Instabilities

Avalenche

Dislocations

Plastic deformation bands

## ABSTRACT

Plastic flow is conventionally treated as continuous in finite element (FE) codes, whether in isotropic, anisotropic plasticity, or crystal plasticity. This approach, derived from continuum mechanics, contradicts the intermittent nature of plasticity at the elementary scale. Understanding crystal plasticity at micro-scale opens the door to new engineering applications, such as microscale machining. In this work, a new approach is proposed to account for the intermittence of plastic deformation while remaining within the framework of continuum mechanics. We introduce a material parameter, the plastic deformation threshold, denoted as $\Delta p_{min}$, corresponding to the plastic deformation carried by the minimal plastic deformation burst within the material. The incremental model is based on the traditional predictor-corrector algorithm to calculate the elastoplastic behavior of a material subjected to any external loading. The model is presented within the framework of small deformations for von Mises plasticity. To highlight the main features of the approach, the plastic strain increment is calculated using normality rule and consistency conditions, and is accepted only if it exceeds $\Delta p_{min}$. To achieve this, a time-discontinuous generalization of the Karush-Kuhn-Tucker (KKT) conditions is proposed. The simulations show that the introduction of the plastic threshold allows for the reproduction of the spatiotemporal intermittence of plastic flow, capturing the self-organization of plastic flow in complex loading scenarios within an FE model.

## 1. Introduction

Plastic flow is usually modeled as a continuous spatio-temporal deformation mechanism in engineering finite element (FE) codes, be it in the context of isotropic, anisotropic, or crystal plasticity. This assumption underlying continuum plasticity models contradicts the intermittent nature of plasticity at elementary scales (Brown, 2012). This intermittency has been experimentally observed in various ways. For instance, during the compression of micropillars, the stress-strain curve exhibits significant serrations, which is related to the abrupt activation of a small number of slip planes (Zhang et al., 2017; Uchic et al., 2009). Similar results are observed during the tension of submicrometer single crystals (Oh et al., 2009; Kiener et al., 2008; Borasi et al., 2023). Acoustic emission and high-resolution extensometry experiments on centimeter-sized single-crystal samples also reveal the organized, non-chaotic nature of plastic flow due to dislocation avalanches and their interactions (Weiss et al., 2007). Digital image correlation (DIC) at microscopic scale also demonstrates the presence of localized deformation bands at grain scale in macro-size polycrystalline materials due to local hardening (nano-precipitates) followed by softening (Charpagne

et al., 2021; Marano et al., 2019). When the tested specimen is macroscopic, or the material contains a high density of grain boundaries, dislocation avalanches are hindered, and the flow becomes quasi-continuous, explaining why it has not been necessary to model this intermittency in traditional macroscopic industrial problems. However, with advances in nanotechnology and the increasing production of micro-metric components, developing a model that accounts for the spatial and temporal intermittency of plasticity and is usable in FE codes becomes necessary (Zepeda-Ruiz et al., 2021). Understanding crystalline plasticity at this scale opens the door to a wide range of new engineering applications, such as micro-machining and the design of novel materials like hierarchical steels (Xu et al., 2022b; Zhu et al., 2016; Xu et al., 2022a).

Let us acknowledge that at atomic scales, molecular dynamics (MD) (Fernández-Castellanos et al., 2021; Patinet et al., 2011) and discrete dislocation dynamics (DDD) (Gómez-García et al., 2006; Csikor et al., 2007) models accurately and “naturally” reproduce the correlated nature of dislocation avalanches and their power law amplitude

<footnote><sup>a</sup> Mines Paris, PSL University, Centre des Matériaux, CNRS UMR 7633, BP 87, 91003, Evry, France</footnote>

<footnote><sup>b</sup> CNRS, LSPM UPR 3407, Université Sorbonne Paris Nord, 93430, Villetaneuse, France</footnote>

<footnote>* Corresponding author.
E-mail address: mathias.lamari@minesparis.psl.eu (M. Lamari).</footnote>
distribution (Weiss et al., 2007). While these models help understand the physics behind the intermittency and organization of plastic flow, they cannot be used at industrial scales due to excessively high computation times.

Several microscopic plasticity models have been recently proposed in the literature to account for the statistical structure of intermittent acoustic emission generated by plastic flow bursts without deviating (too much) from the continuum mechanics formalism. A noteworthy example of such models is the mesoscale tensorial model (MTM) (Salman and Truskinovsky, 2011, 2012; Baggio et al., 2023b,a). Models that attempt to introduce intermittency into crystal plasticity also exist, typically relying on the introduction of stochasticity in an internal system variable (dislocation density Weiss et al., 2015, additional stress added to the yield limit Wijnen et al., 2021; Vermeij et al., 2024, critical resolved shear stress Gélébart, 2021) or within the solver itself (Yu et al., 2021b,a). These modifications succeed in predicting the intense localization of plastic deformation into finite bands, and size dependence of the yield stress in microscopic specimens. However, they rely on introducing probabilities into the crystalline plasticity model itself and thus require prior knowledge of the laws governing the correlations of plastic flows. They are also computationally expensive due to the stochastic treatment at each time step. Recently, Ruestes and Segurado (2024) used a kinetic Monte Carlo selection process to control the activation of microscopic slip mechanisms, achieving similar results. Without introducing any stochasticity, Ryś et al. (2024) have demonstrated that spontaneous localization of plastic strain can occur in crystal plasticity due to the yield-vertex effect when the latent-hardening of slip systems is higher than the self-hardening.

In this work, we propose a new approach to account for the intermittency of plastic deformation that does not rely on the prior knowledge of the statistics of plastic flow bursts. We start from a standard isotropic rate-independent J2-plasticity formulation, which is implemented in a finite element solver using the standard implicit time-stepping with radial return algorithm. We introduce a single strictly positive parameter, the plastic strain threshold denoted as $\Delta p_{\min}$, which corresponds to the smallest (cumulative) plastic strain increment that the material may undergo at a particular time under load. This plastic strain threshold is physically associated to the smallest dislocation avalanche that can develop. Solving the equation of elasto-plasticity under the constraint that the plastic increment cannot be smaller than the plastic threshold is performed with very limited modification of the standard numerical formulation of departure. Typically, our implementation in the FEChiCSX FEA library (Scroggs et al., 2022) requires modifying a single line of code to take this inequality constraint into account and depart from standard continuum finite element plasticity. Although the proposed model was derived by effectively manipulating a standard continuum plasticity finite element code to trigger the intermittence effects, we propose a set of constitutive equations for the time-discontinuous plasticity model, of which the modified FEA model mentioned previously is a rigorous implicit time-stepping implementation. Our presentation of the method will then follow a rather classical mode of exposure, starting from the continuum equations before deriving the associated fully discrete finite element model. A validation of the approach is made by implementing the model in another finite element code, namely Zset (Besson and Foerch, 1997), thus showing that the model goes beyond a particular implementation.

The simulations performed using this modified, discontinuous, engineering plasticity formulation are remarquable in several ways. Noteworthy, the plastic flow arises in the form of plastic bands, which results in serrations on the load-displacement curve of tensile specimens. These bands are akin to the one observed in Portevin-Le-Chatelier (PLC) simulations (Colas et al., 2014; Wang et al., 2011; Lamari et al., 2024; Guillermin et al., 2023), but in contrast to strain ageing effects, they result from a rate-independent model formulation. In this sense, the phenomenology of the plastic flow generated by the proposed formulation shares resemblance with the Lüders effect, that would in

our case be triggered repeatedly rather than once in the case of Lüders band propagation. Remarkably, the proposed model is discontinuous but does not exhibit features of instability encountered in softening plasticity models. There is no dependency of the band statistics (intensity and width) to mesh size (providing that the level of refinement is sufficient to describe strain gradients). In addition, it will be shown that the Newton solver associated with the finite element force balance equations converges in a well-behaved manner, exhibiting the expected quadratic local convergence.

The paper is organized as follows. In Section 2, the physical motivation of the proposed model is given, followed by a presentation of the time-discontinuous model and its FE discretization. In Section 3, a first simulation under uniform fields is presented, followed by a reference simulation with complex dogbone geometry. The methodology to analyze and quantify the complex plastic activity appearing in the simulations is then given. Section 4 is dedicated to the use of this methodology to analyze the effect of different simulations and model parameters on the model response. Finally, Section 5 is devoted to the analysis of power law distributed plastic events, spatio-temporal analysis of plastic avalanches. Finally, the results of a tension test of a plate with a hole are discussed to highlight the potential of the model for application to structural computations.

## 2. Model formulation

### 2.1. Physical foundation

Recently, Perchikov and Truskinovsky have emphasized the fundamentally quantized nature of plastic deformation in crystalline materials (Perchikov and Truskinovsky, 2024), which occurs in a time-independent manner. In most metallic crystalline materials at room temperature, the most common deformation mechanism is dislocation slip. The displacement of a dislocation line induces a local displacement of $b$, which is the norm of the Burgers vector (typically $\approx 2 \times 10^{-10}$ m in iron Lamari et al., 2024). In micro-scale samples (sample size $L \approx 10^{-6}$ m) the displacement of a dislocation from one Peierls valley to the next induces a plastic deformation of $\approx 2 \times 10^{-4}$, which is small but not infinitesimal. In situ observation of dislocation gliding shows a phenomenon that is very abrupt, occurring at the speed of sound in metals (Oh et al., 2009; Yu et al., 2015). From the point of view of the observer it can be seen as instantaneous. Therefore, in this work we will present a model that is time-discontinuous, time-independent and with a minimum plastic strain increment, called the plastic threshold and denoted $\Delta p_{\min}$ with a value of $\approx 2 \times 10^{-4}$.

In single crystals, the absence of grain boundary screening of the dislocation long-range stress makes the displacement of one dislocation prone to induce motion of other dislocations in parallel slip planes, creating an avalanche of plastic deformation (Brown, 2012). In micro-scale materials, dislocation sources are rare and hindered by surface effects (Zhang et al., 2017; Uchic et al., 2009). This makes the system prone to reach self-organized criticality, leading to larger dislocation avalanches and scale-free intermittent plasticity (Dimiduk et al., 2006). Therefore, the model we propose must be scale-free and prone to promote plastic strain bursts (i.e. plastic deformation increment $> \Delta p_{\min}$).

While the model is primarily intended for crystal plasticity to be comparable with metallic crystal experiments, for simplicity we will propose our model in the framework of isotropic J2 plasticity. This is not completely unphysical, since intermittent plasticity is also found in metallic glasses (Okuyucu et al., 2023). While dislocation gliding is not the mechanism of plastic deformation in metallic glasses, the various mechanisms (Greer et al., 2013) such as the shear transformation zone (STZ) (Chevalier et al., 2018; Yan et al., 2010) also induce plastic deformation that is not locally infinitesimal. Plastic bands are also observed. The plastic threshold concept still applies and we will take its value to be $\approx 2 \times 10^{-4}$ for simplicity.
In the following, the terms plastic burst or plastic jump are primarily used to describe significant plastic events occurring over a substantial region of the specimen. The term avalanche is used occasionally, in line with its usage in the physics community, despite the fact that the simulation does not directly involve dislocation concepts and that the plastic burst occurs in only one numerical step, as will be shown below.

## 2.2. Continuum plasticity

In this section, the equations of basic continuum elasto-plasticity are recalled for a classical rate-independent J2-plasticity including an isotropic hardening law. This will clarify where the time-discontinuous elasto-plasticity model introduced in the following section diverges from the classical continuous framework.

**Linearized continuum mechanics.** The equations of elasto-plasticity are solved over a space-time domain $\Omega \times \mathcal{T}$ where $\Omega \subset \mathbb{R}^3$ and $\mathcal{T} = [0, \mathcal{T}]$. Let $x \in \mathbb{R}^3$ denote the spatial coordinates of a point of $\Omega$ and let $t \in \mathcal{T}$ denote the time variable. The boundary $\partial\Omega$ is split into a part $\partial_u\Omega$ over which time-dependent Dirichlet conditions are applied, and a part $\partial_t\Omega$ over which homogeneous Neumann conditions are applied.

The displacement $(x, t) \mapsto u(x, t)$ and Cauchy stress tensor field $(x, t) \mapsto \sigma(x, t)$ satisfies the following equations:

$$ \operatorname{div} \sigma = 0 \quad \text{in } \Omega \qquad (1) $$

$$ u = u_d \quad \text{over } \partial_u \Omega \qquad (2) $$

$$ \boldsymbol{\sigma} \cdot \boldsymbol{n}_{\partial\Omega} = 0 \quad \text{over } \partial_t \Omega \qquad (3) $$

where $n_{\partial\Omega}$ denotes the outer normal to the domain boundary and $u_d$ is a known function defined over $\partial_u\Omega \times \mathcal{T}$ and with values in $\mathbb{R}^3$. Volume forces, including inertial forces are excluded for simplicity.

The infinitesimal strain tensor is defined as

$$ \epsilon = \frac{1}{2} (\nabla u + (\nabla u)^T) \qquad (4) $$

To close the system, a constitutive relation linking $\sigma$ to the time-history of $\epsilon$ needs to be introduced.

**Constitutive equations.** Two additional unknown fields of internal variables are introduced, namely the second-order plastic strain tensor $(x, t) \mapsto \epsilon^P(x, t)$ and the scalar cumulative plastic strain $(x, t) \mapsto p(x, t)$. The constitutive relation of the elasto-plastic material is written as follows. The strain tensor can be decomposed into elastic and plastic parts

$$ \epsilon = \epsilon^e + \epsilon^P \qquad (5) $$

The Cauchy stress is given by Hooke's law for isotropic elasticity

$$ \boldsymbol{\sigma} = \mathbb{C} : (\epsilon - \epsilon^P) = \lambda \operatorname{tr}(\epsilon^e) \mathbf{I} + 2\mu \epsilon^e \qquad (6) $$

where $\mathbb{C}$ is the fourth-order elasticity tensor, $\lambda, \mu$ are Lamé coefficients and $\operatorname{tr}(\epsilon^e)$ denotes the trace of tensor $\epsilon^e$. The yield surface is defined as

$$ f(\boldsymbol{\sigma}; p) = \sigma_{vM}(\boldsymbol{\sigma}) - R(p) - \sigma_{ys} \qquad (7) $$

where $R(p)$ is the isotropic hardening function which is assumed to be monotonically non-decreasing, $\sigma_y$ is the initial yield stress and $\sigma_{vM}$ the von Mises equivalent stress ($\sigma_{vM} = \sqrt{\frac{3}{2}s} : s$ where $s$ is the deviatoric part of the Cauchy stress tensor). The normal to the yield surface is defined as $\boldsymbol{n} = \partial f / \partial \boldsymbol{\sigma}$. The normality rule for the plastic flow is:

$$ \dot{\boldsymbol{\epsilon}}^{P} = \boldsymbol{\lambda} \boldsymbol{n} \qquad (8) $$

The system is closed by the Karush-Kuhn-Tucker (KKT) conditions:

$$ f \lambda = 0 \qquad (9) $$

$$ \lambda \ge 0 \qquad (10) $$

$$ f \le 0 \qquad (11) $$

The time-integrated multiplier coincides with the cumulative plastic strain $p$ in J2-plasticity. This specificity reads as:

$$ \lambda = \dot{p} = \sqrt{\frac{2}{3} \dot{\epsilon}^{P} : \dot{\epsilon}^{P}} \qquad (12) $$

**Initial conditions.** Finally, initial conditions need to be applied to the fields of internal variables, which take the form $p|_{t=0} = 0$ and $\epsilon^p|_{t=0} = 0$.

## 2.3. Time-discontinuous plasticity model

In the time-discontinuous plasticity model, the cumulative plastic strain variable $p$ is allowed to increase only by instantaneous events (bursts) of finite magnitude. To do this, mechanical quantities are allowed to operate beyond the yield surface, i.e. violate the condition (11), and return to the yield surface when and only when the resulting instantaneous increment of cumulative plastic strain is larger than $\Delta p_{\min} \in \mathbb{R}^{**}$. The plastic threshold $\Delta p_{\min}$ is the sole additional parameter of the model, and represents the smallest plastic strain increment that the material system can accommodate.

To achieve this mathematically, Eq. (8), (9) and (11) are modified, while the others are retained. The first equation is altered as follows

$$ \llbracket \epsilon^P \rrbracket = \llbracket p \rrbracket n^{-} \qquad (13) $$

In this equation, $\llbracket \cdot \rrbracket = \lim_{h \to 0} [\cdot(t + h) - \cdot(t - h)]$, with $\cdot$ being any quantity such as $p$ or $\epsilon^P$, indicates a jump of the considered quantity at time $t$. Notation $n^{-} = \lim_{h \to 0} n(t - h)$ denotes the normal to the yield surface at time $t$, taking the limit from the left if the normal evolves in a discontinuous manner in time, as this corresponds to the stress state before the plastic burst.<fnref n="1" /> Physically, this means that the stress state before the plastic burst occurs determines the direction of plastic yield.

At each time, the plastic strain increment must satisfy new constitutive equations. The set $P^{(t)}$ is defined as the collection of admissible plastic bursts that comply with the constitutive relations presented herein. Admissible cumulative plastic strain increments that belong to the set $P^{(t)}$ will be denoted as $\llbracket p \rrbracket^*$. The corresponding evolution of yield surface after the increment $\llbracket p \rrbracket^*$ is denoted $f^{**}$, defined as $f^{**} = f(\mathbb{C} : (\epsilon - \epsilon^P - \llbracket p \rrbracket^* n^{-}), p^{-} + \llbracket p \rrbracket^*)$, with $p^{-} = \lim_{h \to 0} p(t - h)$.

The first new equation of the time-discontinuous plasticity model that controls the value of $\llbracket p \rrbracket^*$ is:

$$ \llbracket p \rrbracket^* f^{**} = 0 \qquad (14) $$

Eq. (14) is used to express the fact that if there is a non-vanishing instantaneous plastic flow increment at time $t$, its amplitude must be such that the mechanical stress returns to the yield surface.

Finally, the plastic strain increment must obey a last equation assuring that the system is closed:

$$ \llbracket p \rrbracket^* \langle \Delta p_{\min} - \llbracket p \rrbracket^* \rangle = 0 \qquad (15) $$

where $\langle \cdot, \cdot \rangle = \max\{0, \cdot, \cdot\}$. This equation means that either the plastic flow is null, or the increment of cumulative plastic strain is larger than $\Delta p_{\min}$, which introduces flow discontinuities in the model.

Then, at time $t \in \mathcal{T}$, the jump of cumulative plastic strain is chosen as the maximum of all admissible cumulative plastic strains:

$$ \llbracket p \rrbracket = \max_{\llbracket p \rrbracket^* \in P^{(t)}} \llbracket p \rrbracket^* \qquad (16) $$

Indeed, the set of Eqs. (14) and (15) introduced previously leaves an indeterminacy. This is because $\llbracket p \rrbracket = 0$ is always a trivial solution. The

<footnote n="1">¹ Notice that in J2-plasticity, the normal tensor remains unchanged during a plastic burst, i.e.</footnote>

<footnote>$\llbracket n \rrbracket = 0$.</footnote>
max operation in Eq. (16) specifies that if a jump in plastic strain is admissible and positive, then it must happen.

The term “time-discontinuous” refers to the fact that when a plastic burst occurs, certain quantities experience a discontinuity in their temporal evolution, such as p, f, εp due to Eqs. (14) and (15). However, this is not the case for all quantities, such as σ, εe, ε or n which may exhibit either time-continuous or discontinuous behavior depending on the external conditions.

## 2.4. Algorithmic implementation

This section outlines the algorithm used to implement the time-discontinuous model. A straightforward modification to the radial return algorithm ensures that plastic strain follows the new constitutive equations from Eqs. (13) to (16), with the remainder of the procedure remaining classical. A comprehensive description of the model’s implementation in the finite element solver, the time discretization of all constitutive relations, and the integration of the global Newton solver is provided in Appendix A.

The time interval $T = [0, T]$ is divided into $N$ time steps. The equation of continuum mechanics is enforced at discrete times $\bar{T} = \{t_0, t_1, ..., t_N\}$, using the finite element method. In the following, $n$ will be an integer between 0 and $N$. As is customary, the quantity $X$ at step $n$ will be denoted $X_n$.

Suppose the system’s evolution is known up to step $n$. The output of stress function $\sigma_{n+1}(\varepsilon_{n+1}; \varepsilon_n^p, p_n)$ is evaluated via an operator-splitting procedure akin to those that are traditionally employed in implicit schemes for classical rate-independent elasto-plastic constitutive laws (Besson et al., 2009; de Souza Neto et al., 2011; Simo, 1998). This is done as follows, at time $t_{n+1}$:

* Determine the direction of the plastic flow

$$ \mathbf{n}_{n+1}^{-} = \frac{\partial f}{\partial \sigma} \Big|_{\sigma_{n+1} = C:(\varepsilon_{n+1} - \varepsilon_n^p), p_{n+1} = p_n} \quad (17) $$

* Evaluate yield function $f^*$ assuming that no plastic increment occurs.

$$ f_{n+1}^* = \sigma_{eq} (C : (\varepsilon_{n+1} - \varepsilon_n^p)) - R(p_n) - \sigma_y \quad (18) $$

* Discriminate between elastic and plastic behavior:

- if $f_{n+1}^* \le 0$, set $p_{n+1} = p_n$,

- if $f_{n+1}^* > 0$, perform the following two steps:

* assume that a plastic jump occurs within time step $]t_n, t_{n+1}]$, i.e. $(p_{n+1} - p_n) \ne 0$ (in that case, Eq. (14) specifies that the stress must return to the yield limit at the end of the time step). Look for test plastic increment $\Delta p^*$ satisfying

$$ \sigma_{eq} (C : (\varepsilon_{n+1} - \varepsilon_n^p - \Delta p^* n_{n+1}^{-})) - R(p_n + \Delta p^*) - \sigma_{ys} = 0 \quad (19) $$

using a Newton-Raphson algorithm.

* Compare $\Delta p^*$ to $\Delta p_{min}$. If $\Delta p^* < \Delta p_{min}$, set $p_{n+1} = p_n$, otherwise set $p_{n+1} = p_n + \Delta p^*$

* Compute $\varepsilon_{n+1}^p = \varepsilon_n^p + (p_{n+1} - p_n)n_{n+1}^{-}$ and return $\sigma_{n+1} = C : (\varepsilon_{n+1} - \varepsilon_{n+1}^p)$

A pseudo-code of the radial return algorithm as implemented in the FENiCSX code is given in Algorithm 1. Compared to a classical radial return algorithm used in continuum plasticity, only lines 14 to 17 have been added. Removing these lines will revert to the classical plasticity approach. The full script used in our FENiCSX simulation, including

geometry and example runs, can be found in this GitHub repository. Due to limitations related to automatic differentiation in FENiCSX, which does not support loops, the maximum number of iterations in the Newton loops, $N_{Max}$, has been set to 1. However, comparisons with Zset simulations show that this does not significantly affect the simulation results, even when the work-hardening behavior is not purely linear.

### Algorithm 1 Time-discontinuous radial return (von Mises)

1: Initialize: $\epsilon, \epsilon^p, p, C, \Delta p_{min}, \mu, R, \delta, N_{max}$

2: Calculate: $\sigma^* = C : (\epsilon - \epsilon^p)$

3: Calculate: $f = \sigma_{vM}(\sigma^*) - \sigma_{ys} - R(p)$

4: Calculate: $n^- = \partial f / \partial \sigma^*$

5: Set: $i = 0$

6: Set: $\Delta p^* = 0$

7: while $|f| \ge \delta$ and $i \le N_{max}$ **do**

8: Calculate: $f' = -3\mu - \frac{\partial R}{\partial p}(p + \Delta p^*)$

9: Calculate: $\Delta p^* = -f/f'$

10: Calculate: $\sigma^* = C : (\epsilon - (\epsilon^p + \Delta p \cdot n^-))$

11: Calculate: $f = \sigma_{vM}(\sigma^*) - \sigma_{ys} - R(p + \Delta p^*)$

12: Calculate: $i = i + 1$

13: end while

14: if $\Delta p^* \ge \Delta p_{min}$ **then**

15: $\Delta p = \Delta p^*$

16: else

17: $\Delta p = 0$

18: end if

19: if $f \ge 0$ **then**

20: $\Delta \epsilon^p = \Delta p \cdot n^-$

21: else

22: $\Delta \epsilon^p = 0 \cdot n^-$

23: end if

24: Update: $\epsilon^p = \epsilon^p + \Delta \epsilon^p$

25: Update: $p = p + \sqrt{\frac{2}{3}} \Delta \epsilon^p : \Delta \epsilon^p$

26: Return: $\epsilon^p, p$

## 2.5. Equivalence with a two-yield surface model

This section demonstrates the equivalence of the time-discontinuous plasticity model to a two-yield surface model. Note that all new interpretations presented here are not additional properties of the model but are instead derived from previous equations.

It is assumed that, at time $t$, a plastic burst occurs (i.e. $\|s\| > 0$) in a region of the system. The tensors $e$ and $e^e$ are the deviatoric parts of the total and elastic strain tensors ($\epsilon^p$ is always deviatoric in von Mises plasticity). Hooke’s law provides the deviatoric stress component at the end of the plastic burst:

$$ s + \|s\| = s + 2\mu\|e^e\| = s + 2\mu(\|e\| - \|\epsilon^p\|) = s^{\text{trial}} - 2\mu\|\epsilon^p\| \quad (20) $$

In Eq. (20), the trial stress $s^{trial} = s + 2\mu\|e\|$ is the stress reached if all strain increment is elastic. The von Mises stress corresponding to $s^{trial}$ is denoted $\sigma_{vM}^{trial}$. From Eq. (20), the tensor $s^{trial}$ is the linear combination of $s + \|s\|$ and $\|\epsilon^p\|$. From Eq. (13), $\|\epsilon^p\|$ is colinear to $s + \|s\|$, since the normal to the yield surface $n$ is the same before and after plastic relaxation. Then $s^{trial}$ is colinear to $s + \|s\|$. Therefore, the following normalized tensors are equal:

$$ \sqrt{\frac{3}{2} \frac{s + \|s\|}{\sigma_{vM}^{trial}}} = \sqrt{\frac{3}{2} \frac{s^{trial}}{\sigma_{vM}^{trial}}} \quad \text{where} \quad \sigma_{vM}^{trial} = \sqrt{\frac{3}{2} (s + \|s\|) : (s + \|s\|)} \quad (21) $$

We deduce from the last Eqs. (20) and (21), and the normality rule of plastic flow Eq. (13) that:

$$ \sigma_{vM}^{trial} = \sigma_{vM}^{trial} - 3\mu \|p\| \quad (22) $$
<figure><img src="image_4.png" /><figcaption>Fig. 1. Representation of the lower and upper yield surfaces in the principal stress basis. Yielding occurs when a material point exceeds the upper yield surface, and the plastic strain increment is determined relative to the lower yield surface.</figcaption></figure>

After the plastic jump, the yield surface function is null:

$$ \sigma_{\text{vM}}^{+} = \sigma_{\text{ys}} + R(p + \llbracket p \rrbracket) \quad (23) $$

Combining Eqs. (22) and (23), the relationship between $\sigma_{\text{vM}}^{\text{trial}}$ and $\llbracket p \rrbracket$ becomes straightforward:

$$ \sigma_{\text{vM}}^{\text{trial}} = \sigma_{\text{ys}} + R(p + \llbracket p \rrbracket) + 3\mu \llbracket p \rrbracket \quad (24) $$

By substituting the condition $\llbracket p \rrbracket \ge \Delta p_{\min}$ in Eq. (24), the following condition on the von Mises stress is derived to predict the occurrence of plastic burst:

$$ \sigma_{\text{vM}}^{\text{trial}} \ge \sigma_{\text{ys}} + R(p + \Delta p_{\min}) + 3\mu \Delta p_{\min} \quad (25) $$

The latter condition is equivalent to defining an upper yield surface that describes the onset of yielding. Linear isotropic hardening is assumed for the remainder of the article to simplify the interpretation of simulations, i.e. $R(p) = H p$, $H$ being the hardening modulus. Note that the model can be applied to any form of non-decreasing hardening. Eq. (26) is then simplified:

$$ \sigma_{\text{vM}}^{\text{trial}} \ge \sigma_{\text{ys}} + H p + (3\mu + H) \Delta p_{\min} \quad (26) $$

It means that as soon as $\sigma_{\text{vM}}^{\text{trial}}$ reaches this critical level, a plastic burst takes place. In Fig. 1, a representation of the lower and upper yield surfaces in the principal stress space is given. Yielding occurs when the stress in a material point of the system exceeds the upper yield surface, and the plastic strain increment is determined relative to the lower yield surface. The difference between the two yield surfaces is constant and equal to $(3\mu + H) \Delta p_{\min}$. This constant gap explains visually why in the time-discontinuous model, the plastic strain increment cannot be infinitesimal for a given $\Delta p_{\min}$.

## 3. Results

To explore qualitatively and quantitatively the properties of our time-discontinuous model, we first simulate a simple 1D geometry under uniaxial tension. Then, an analytical solution is found for homogeneous deformation under any triaxiality condition and compared with simulations. Finally, more complex geometries are tested under uniaxial tension. The model has been implemented independently in the FE niCSX (Scroggs et al., 2022) and Zset (Besson and Foerch, 1997) solvers and both solvers deliver the same results.

### 3.1. Simulations with uniform fields

A one-dimensional geometry consisting of N hexahedral elements aligned along the $O_x$ axis is used in this section, as depicted in Fig. 2(a). On faces orthogonal to the $O_x$ axis, mixed conditions are imposed. At each numerical step, the displacement increment component along $O_x$ is fixed (0 for the left face at the origin, $\Delta u$ for the right face), while the other components are left free, except at points A and O, to suitably fix the rigid body motion (details in Fig. 2(a)). The other surfaces are free. The initial yield stress $\sigma_{\text{ys}}$ is set to 100 MPa, the linear hardening rate $H$ to 10 GPa, the Young modulus $E$ to 200 GPa, and the Poisson ratio $\nu$ to 0.3. The imposed displacement step $\Delta u$ corresponds to a strain increment $\Delta \epsilon_{xx}$ in the tensile direction of $3 \times 10^{-6}$ at each loading step. The plastic threshold $\Delta p_{\min}$ is set to $2 \times 10^{-4}$. For the simulations to accurately reflect the time-discontinuous model under displacement boundary conditions, it is necessary to have, for all component i,j, that $\Delta \epsilon_{ij} \ll \Delta p_{\min}$.

The first noteworthy feature of the simulation is the spatial homogeneity of both stress and strain fields across the entire geometry under tensile loading. Without introducing defects in the mesh or the elastic/plastic properties of the material, there is no localization of stress, regardless of the number of elements. Under tensile condition, the stress tensor has only one nonzero component along the $O_x$ axis. This test is equivalent to testing a single Gauss point, but this setup was used in order to demonstrate the absence of localization in this particular case.

The tensile curve (stress $\sigma_{xx}$ vs strain $\epsilon_{xx}$) is shown in Fig. 2(b). The evolution of the lower (red) and upper (purple) yield surfaces as defined by Eqs. (23) and (26) are also incorporated for clarity, even though it must be emphasized that they are not used explicitly in the algorithm. The specimen first deforms elastically until the von Mises stress reaches the upper yield surface. Then, the entire specimen undergoes plastic deformation, causing the stress to decrease through elastic relaxation until it reaches the lower yield surface. The plastic strain burst occurs within one numerical step. Then, the phenomenon repeats but with a hardened material due to linear isotropic hardening. Thus, periodic serrations are observed on the stress-strain curve; each stress increase corresponds to purely elastic loading of the material with a slope equal to the Young modulus; each drop corresponds to plastic deformation of the entire specimen.

Fig. 2(c) gives the evolution of equivalent plastic strain during the simulation. As expected, plastic strain increases by steps, each plastic burst occurring in one numerical increment and corresponding to a stress drop on the tensile curve. Only two values of plastic strain increment $\llbracket p \rrbracket$ are possible at each numerical step, as demonstrated in the next Section 3.2, for any triaxiality conditions. If no plastic deformation occurs, $\llbracket p \rrbracket$ is 0. If a plastic burst occurs, $\llbracket p \rrbracket$ is given by:

$$ \llbracket p \rrbracket = \frac{3\mu + H}{E + H} \Delta p_{\min} \quad (27) $$

The last Eq. (27) is valid under uniaxial tension. When the Poisson ratio $\nu = 0.5$, then $\llbracket p \rrbracket = \Delta p_{\min}$. The drop of von Mises stress at each plastic event $\llbracket \sigma_{\text{vM}} \rrbracket$ is given by:

$$ \llbracket \sigma_{\text{vM}} \rrbracket = -E \llbracket p \rrbracket \quad (28) $$

The Eqs. (27) and (28) are only valid for uniaxial straining.

### 3.2. Analytical homogeneous solution under any triaxiality

In this section, a single Gauss point is considered. The material parameters are chosen to be the same as in Section 3.1. We control the triaxiality factor $T$ defined as the ratio between the hydrostatic stress $\sigma_h$ and the equivalent stress $\sigma_{\text{vM}}$:

$$ T = \frac{\sigma_h}{\sigma_{\text{vM}}} \quad \text{with} \quad \sigma_h = \text{tr}(\sigma)/3 \quad (29) $$
<figure><img src="image_5.png" /><figcaption>Fig. 2. (a) Geometry used for the homogeneous simulations. The number N of aligned hexahedral elements can be set to any integer. Prescribed displacements are imposed on the left and right faces of the sample. Displacement components not indicated are left free. All other faces are free. (b) Evolution of stress $\sigma_{xx}$ as a function of strain $\epsilon_{xx}$ for the homogeneous simulation. The lower yield surface (red) and the upper yield surface (purple) have been indicated. (c) Corresponding evolution of the cumulative plastic strain as a function of applied strain.</figcaption></figure>

For all tested T values, which ranged from 0 (pure shear) to 2/3 (equibiaxial), the results in terms of equivalent stress and equivalent plastic strain are qualitatively identical to what is presented in Fig. 2. Periodic stress drops linked with plastic bursts are found, and both $[\mathbf{p}]$ and $[\sigma_{vM}]$ admit only one given value during plastic bursts, as in the tensile case. For pure shear, Eq. (27) is not satisfied and we have instead $[\mathbf{p}] = \Delta p_{\min}$, independently of the values of E and $\nu$.

To understand this difference and calculate the theoretical value of the stress drop and plastic strain increase during a plastic burst, we consider a time t at which a plastic burst occurs. Independently of the triaxiality factor, the two-yield surfaces concept is respected. The value of the von Mises stress just before the jump corresponds to the value of the upper yield surface prior to the plastic strain increment, i.e.:

$$ \sigma_{vM} = \sigma_{ys} + H p + (3\mu + H)\Delta p_{\min} \quad (30) $$

The value of the equivalent stress after plastic deformation corresponds to the value of the lower yield surface accounting for the plastic strain increment:

$$ \sigma_{vM} + [\sigma_{vM}] = \sigma_{ys} + H(p + [\mathbf{p}]) \quad (31) $$

Therefore:

$$ [\sigma_{vM}] = H[\mathbf{p}] - (3\mu + H)\Delta p_{\min} \quad (32) $$

In the last Eq. (32), $[\sigma_{vM}]$ and $[\mathbf{p}]$ are both unknown, and a second equation is necessary to obtain their values. This second equation is obtained by using the prescribed total deformation. As previously stated in Section 3.1, each component of the strain tensor at each numerical step must be negligible with regards to the plastic threshold $\Delta p_{\min}$. Therefore, even during the plastic burst, the total strain tensor of the system does not change:

$$ [\epsilon] = [\epsilon^c] + [\epsilon^p] = 0 \quad (33) $$

In continuum plasticity, if the strain rate is negligible, then the plastic strain rate and elastic strain rate are both individually negligible. One important difference in our time-discontinuous approach is that while $[\epsilon]$ is zero during the instantaneous plastic burst, $[\epsilon^c]$ and $[\epsilon^p]$ are not infinitesimal, and we have:

$$ [\epsilon^c] = -[\epsilon^p] \quad (34) $$

To obtain a scalar relationship between $[\mathbf{p}]$ and $[\sigma_{vM}]$, the tensor contraction of both sides of Eq. (34) with $[\sigma]$ is taken. The right-hand side is simplified by using successively the normality rule (Eq. (13)), the deviatoric nature of the normal tensor $\mathbf{n}$ and the invariance of $\mathbf{n}$ during the plastic burst:

$$ [\sigma] : [\epsilon^c] = [\mathbf{p}] ([\sigma] : \mathbf{n}) = [\mathbf{p}] ([\mathbf{s}] : \mathbf{n}) = [\mathbf{p}] [\sigma_{vM}] \quad (35) $$

The tensor contraction $[\sigma] : [\epsilon^c]$ is simplified by decomposing the stress and elastic strain tensors into their hydrostatic part $\sigma_h$ and $\epsilon_h^c$ and their deviatoric part $\mathbf{s}$ and $\mathbf{e}^c$, and using the Hooke law ($K$ is the bulk modulus):

$$ [\sigma] : [\epsilon^c] = [\sigma_h] : [\epsilon_h^c] + [\mathbf{s}] : [\epsilon^c] = \frac{1}{3K} [\sigma_h]^2 + \frac{1}{2\mu} [\mathbf{s}]^2 \quad (36) $$

The expression (36) is furthermore simplified using $[\sigma_h]^2 = 3[\sigma_h]^2 = 3T^2[\sigma_{vM}]^2$ and $[\mathbf{s}]^2 = 2/3[\sigma_{vM}]^2$:

$$ [\sigma] : [\epsilon^c] = \left( \frac{T^2}{K} + \frac{1}{3\mu} \right) [\sigma_{vM}]^2 \quad (37) $$

Combining Eqs. (35) and (37), a second relationship between $[\mathbf{p}]$ and $[\sigma_{vM}]$ is thus obtained:

$$ [\sigma_{vM}] = -M_T[\mathbf{p}] \quad \text{with} \quad M_T = \frac{1}{\frac{T^2}{K} + \frac{1}{3\mu}} = \frac{3\mu E}{(27\mu - 9E)T^2 + E} \quad (38) $$

The equivalent elastic modulus $M_T$ introduced in Eq. (37) is the harmonic mean of $K/T^2$ and $3\mu$ and represents the coefficient of proportionality between the drop in von Mises stress and the increase in equivalent plastic strain during a plastic burst. It can only be used when the system has no increase in its total strain. In uniaxial tension ($T=1/3$), $M_T$ becomes $E$, as stated in Eq. (28). In pure shear ($T=0$), $M_T$ becomes $3\mu$. When $T$ tends to infinity (pure hydrostatic state), $M_T$ becomes equivalent to $K/T^2$. Combining Eqs. (32) and (38), the relationship between $[\mathbf{p}]$ and $\Delta p_{\min}$ is found:

$$ [\mathbf{p}] = \frac{3\mu + H}{M_T + H} \Delta p_{\min} \quad (39) $$

Finally, in uniaxial tension, Eq. (39) reduces to Eq. (27), and in pure shear, $[\mathbf{p}]$ is equal to $\Delta p_{\min}$, as found in the simulations. The model response is thus entirely understood under homogeneous deformation, for any triaxiality.

### 3.3. Tensile test on the dogbone geometry

By applying the model to more complex geometries, the macroscopic response becomes very different from the homogeneous case. The example of a 3D dogbone flat specimen is described in this section, whose precise geometry and dimensions are depicted with several 3D meshes tested in Fig. 3. Simulations are performed in 3D and there are 12 elements in the thickness for the finest mesh. The actual mesh used in this section consists of linear tetrahedra and is depicted in Fig. 3(d). The impact of using the other meshes will be discussed in the following section. Uniaxial displacement is imposed on the left and right faces,
<figure><img src="image_6.png" /><figcaption>**Fig. 3.** (a) Dogbone geometry used in the reference simulation. Dimensions used for the reference simulation are: L = 20, gauge length = 14, l = 6, e = 0.25, r = 2. The element size is denoted *h*. (b) Meshes used for testing the model (tetrahedral, hexahedral and inclined hexahedral). Representation of the FE meshing in the zoomed area of (a) for (c) a tetrahedral meshing with *h*/L = 1/100, (d) a tetrahedral meshing with *h*/L = 1/200, (e) a tetrahedral meshing with *h*/L = 1/250, (f) a tetrahedral meshing with *h*/L = 1/400, (g) a hexahedral meshing with *h*/L = 1/100, (h) an inclined hexahedral meshing with *h*/L = 1/250. The mesh used in the tetrahedral cases discussed after have a *h*/L ratio of 1/400.</figcaption></figure>

which are considered clamped (no transverse displacements), and the other faces are free of forces. The elastic and plastic properties of the material, as well as the boundary conditions and the value of the plastic threshold $\Delta p_{\min}$, are the same as those presented in Section 3.1.

The tensile curve is drawn in Fig. 4(a), obtained by volume aver-
aging of $\epsilon_{xx}$ and $\sigma_{xx}$ over the gauge length, accompanied by maps of
cumulative plastic deformation (c) and the corresponding increment (d)
at several mean strain levels. Unlike the 1D case, deformation is not
homogeneous along the specimen. Plastic strain bands of finite width
appear. The width of these bands, defined as the length of elements
along the $O_x$ axis that yield at a given increment, is much larger than
the width of an element $h$ (see Fig. 4(b)). The mean plastic strain
increment in the deformation band is an order of magnitude greater
than the plastic threshold $\Delta p_{\min}$.

On the stress–strain curve in Fig. 4(a), it can be observed that, after
an elastic regime, serrations are present in the curve. Each serration
corresponds to the appearance of a macroscopic plastic strain band on
the specimen. From $\epsilon_1$ = 0.06% to $\epsilon_4$ = 0.17% strain, a stress plateau is
observed, similar to materials exhibiting static strain ageing inducing
a Lüders plateau. Beyond 0.17% strain, the serrations continue. The
plateau effect tends to disappear as the strain increases, and the average
stress evolution tends to align with the linear hardening rate $H$. This
second behavior is analogous to a type B or C PLC effect (Yilmaz,

2011; Lamari et al., 2024), notably observed in materials exhibiting
dynamic strain aging. Such a behavior is also similar to a well known
stick–slip behavior (Perfilyev et al., 2013), as observed in slip or
velocity-weakening friction laws (de Geus and Wyart, 2022).

The features of the stress–strain curve can be understood by ana-
lyzing the dynamics of appearance of plastic bands. An animation of
the evolution of cumulative plastic strain map during the simulation
is available here (Kerfriden, 2024). More videos are available in the
following zenodo repository (10.5281/zenodo.14266823). Snapshots of
cumulative plastic strain and plastic strain increment maps are respec-
tively given in Fig. 4(c) and (d). The first band nucleates at one fillet of
the specimen due to a local von Mises stress concentration. Just before
the first plastic event, a difference of 20 MPa exists between these zones
and the central part of the gauge length. These elements are thus the
first to exceed the upper yield surface and to deform plastically, as
shown at $\epsilon_1$ in Fig. 4(c) and (d). At higher deformation $\epsilon_2$ = 0.07%, they
induce the instantaneous formation of the first traversing band. This
first band induces the yield point phenomenon visible on the tensile
curve. A cut of this first band alongside the $O_x$ axis is shown in Fig. 4(b).
For each band, we define the mean plastic strain, denoted $\bar{\Delta}p$, and the
width of the plastic band, denoted $w$, corresponding to the length of
elements along the $O_x$ axis that yield. More precisely, if $L_w$ is the
<figure><img src="image_7.png" /><figcaption>Fig. 4. (a) Evolution of axial stress as a function of axial strain, obtained by integration over the gauge length of the specimen. (b) Evolution of the plastic strain increment alongside the O_x axis when a plastic burst occurs at a given numerical step. w is the band width, L is the specimen length and Δp is the mean plastic strain in the band. (c) Series of cumulative plastic strain map and (d) series of plastic strain increment over the dogbone specimen for 6 given strain levels indicated in (a).</figcaption></figure>

starting position of the band on the specimen length, l the width of the specimen and e its thickness Δp is calculated by:

$$ \overline{\Delta p} = \frac{1}{w} \int_{L_w}^{L_w+w} \llbracket p \rrbracket \left( x, \frac{l}{2}, \frac{e}{2} \right) dx \quad (40) $$

It is noteworthy that Δp is always higher than Δp_min (here 5 times higher for the first band). It is recalled that only plastic strain increment higher than Δp_min are accepted in the proposed model. At the edges of the first band, there is a local stress concentration, causing the second band to nucleate next to it. This correlated nucleation of bands persists, as illustrated at ε_3 = 0.11%. The phenomenon continues until the entire

specimen plastically deforms for the first time, marking the end of the initial plateau on the tensile curve. This step corresponds to ε_4 = 0.17%. Then, in contrast to the Lüders phenomenon, bands still appear along the gauge length until the end of the simulations, as illustrated for ε_5 = 0.28% and ε_6 = 0.52%. The localization of consecutive bands becomes more random, and they do not occur at constant stress in the presence of hardening.

The convergence of the model is good compared to the continuum model. The reference simulation has been done with time-discontinuous plasticity and with continuum plasticity (i.e., with Δp_min = 0). The computation time for the former is five times that of the latter. Additionally,
<figure><img src="image_8.png" /><figcaption>Fig. 5. Map of plastic strain increment in a region of the specimen for five given plastic events. The instant of their occurrence is indicated on the zoomed part of the tensile curve from Fig. 4. The map height corresponds to the full width of the gauge length.</figcaption></figure>

the Newton solver used in the overall problem resolution converges quadratically once the set of Gauss points that will yield and those that will remain in the elastic regime are determined.

Simulations can be run under force-controlled conditions by replacing displacement with traction boundary conditions. In this case, plastic deformation occurs in bursts at constant force, producing successive plateaus in the tensile curve, similar to force-controlled tests on alloys with the PLC effect (Guillermín et al., 2023) and as observed in micropillar compression (Uchic et al., 2009). Each plateau on the stress curve corresponds to plastic deformation bursts occurring within a single numerical step across the entire gauge length.

### 3.4. Statistical analysis methodology

A close-up view of the tensile curve from Fig. 4 is depicted in Fig. 5, alongside maps of plastic strain increment across the specimen corresponding to the indicated time. Plastic events can occur in two forms. They can either form bands extending over the whole cross-section, as discussed in the previous section, or one/few elements can deform individually. The first kind are called plastic bands, while the second will be referred to as small plastic events hereafter. Their respective impact on the tensile curve is different. While the plastic bands create a drop in stress, the small plastic events only leads to a change in slope, often barely perceptible to the naked eye.

To capture the effect of each plastic events on the stress curves, we calculate at each increment $n$ of the simulation the quantity $\Delta\sigma = -(\sigma_{xx}^{n+1} - \sigma_{xx}^n) + E\Delta\epsilon_{xx}$, with $\Delta\epsilon_{xx} = \epsilon_{xx}^{n+1} - \epsilon_{xx}^n$. Stress and strain values along $O_x$ are obtained through volumic integration over the gauge length. During the simulation, $\Delta\epsilon_{xx}$ (obtained by integration over the volume of the gauge length) remains almost constant at $3 \cdot 10^{-6}$. The quantity $\Delta\sigma$ is zero only if there is no plastic event in the gauge length. Otherwise, it is positive. Thus, $\Delta\sigma$ is an adequate measure of the effect of plastic deformation on the average stress.

The distribution of $\Delta\sigma$ obtained in the simulation detailed in Section 3.3 is presented in Fig. 6 in black in both linear (a) and logarithmic (b) scales. The inset of (b) shows the complementary cumulative density function (CCDF) in logarithmic scale. The distribution of $\Delta\sigma$ consists of two distinct parts: above a stress drop of $\Delta\sigma \approx 2.5$ MPa, it is bell-shaped, as observed on the linear scale; and below 2.5 MPa, it has a power law shape, visible on the logarithmic scale by its linear appearance. On the logarithmic scale, the shape of the distribution is reminiscent of supercriticality (Zhang et al., 2020; Richeton et al., 2005), where the distribution of events follows a power law, except for extreme events that are overrepresented due to a favoring mechanism (Dragon Kings phenomenon, Sornette and Ouillon, 2012). However, the distribution of big events following a well-defined bell curve pleads

for a difference in the physical origin of the two populations of plastic events. Therefore, in this paper we will characterize both parts of the distributions independently. The small events will be analyzed using a truncated power law. The big events, having a well-defined bell shape, will be analyzed using a Gaussian distribution. The overall distribution of $\Delta\sigma$ is divided into two parts using a cut-off value, denoted $\Delta\sigma_{cut}$, determined empirically from the simulation (here 2.5 MPa). The probability density function (PDF) of the stress drop $P(\Delta\sigma)$ is given by:

$$ \text{if } \Delta\sigma > \Delta\sigma_{\text{cut}}: \quad P(\Delta\sigma) = H_d \cdot \left( \frac{1}{\sigma_d \sqrt{2\pi}} \right) \cdot \exp \left( - \frac{(\Delta\sigma - \mu_d^\sigma)^2}{2\sigma_d^2} \right) \quad (41) $$

$$ \text{if } \Delta\sigma < \Delta\sigma_{\text{cut}}: \quad P(\Delta\sigma) = C_d \cdot \Delta\sigma^{-\alpha} e^{-\lambda\Delta\sigma} \quad (42) $$

In Eq. (42), the parameters $\mu_d^\sigma$ and $\sigma_d$ are the mean and standard deviation of the Gaussian distribution, the parameter $\alpha$ is the exponent of the power law, $\lambda$ is a parameter that accounts for the deviation from pure power law due to the non-infinite system size (Zhang et al., 2020), and $H_d$ and $C_d$ are two constants quantifying the respective weight of the two contributions to the PDF. To obtain those parameters, the distribution of small events ($< \Delta\sigma_{cut}$) is first calculated. The exponent $\alpha$ and parameter $\lambda$ are then extracted using the maximum likelihood method (Zhang et al., 2020; Clauset et al., 2009), which was implemented using the “powerlaw” library in python (Alstott et al., 2014). The latter library also allows to evaluate the likelihood of two given distributions to represent a distribution. Truncated power law was better suited than power law, exponential, stretched exponential and log-normal distribution to model the distribution. To obtain $\mu_d^\sigma$ and $\sigma_d$, the distribution of big events ($> \Delta\sigma_{cut}$) is calculated. The method of least squares was then used, with a value of $H_d = 1$ being imposed. Finally, the values of $H_d$ and $C_d$ are chosen so that the integral of the two modeled distributions is equal to the integral of the corresponding part of the distributions of all events.

A problem with the choice of the Gaussian function to model the PDF of $\overline{\Delta p}$ and $w$ is the nonzero probability of negative events, which is unphysical. However, this probability is negligible. Furthermore, simple functions with strictly positive support, such as log-normal distributions, are unsuitable to model the observed distribution because they are long-tailed and non-symmetric. Gaussian distribution was therefore used.

For simplicity, we will focus our analysis on the $\mu_d^\sigma$ and $\alpha$ parameters to characterize the evolution of each portion of $P(\Delta\sigma)$ as we change the model parameters. Respectively, $\mu_d^\sigma$ represents the mean stress drop associated with plastic bands, and $\alpha$ characterizes the nature of the plastic flow associated with small plastic events. Thus, they are the only parameters of $P(\Delta\sigma)$ that will be analyzed afterwards.
<figure><img src="image_9.png" /><figcaption>Fig. 6. Statistics of stress drops (black) determined from the reference simulation in (a) linear scale and (b) logarithmic scale. In red is indicated the modeled distribution using Eq. (42) (the dashed line is the truncated power law, while the solid line is the Gaussian distribution). The inset of (b) shows the complementary cumulative density function (CCDF) and the power law. There are 5922 events in the data set.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Fig. 7. Statistics of (a) mean plastic strain and (b) band widths (black bar chart) determined from the reference simulation. In solid line is indicated the modeled distribution using a Gaussian distribution. There are 187 events in the data set.</figcaption></figure>

Finally, the plastic activity over the $O_x$ axis has also been analyzed. The distributions of mean plastic strain $\Delta p$ (as defined in Eq. (40)) and the width of the plastic bands $w$ (see Fig. 4(b) for a graphical representation) have been determined during the simulations. With this methodology, all plastic bands are characterized, since they all cross the $O_x$ axis, but only a fraction of the small plastic events are captured, i.e. the ones situated on the $O_x$ axis. To further remove from the statistics the effect of small plastic events, a threshold is done on the value of $\bar{\Delta p}$, keeping only those above the threshold. Typically, a plastic event is considered a band when $\bar{\Delta p} > 3\Delta p_{min}$. Hence, only the plastic activity of plastic bands are analyzed afterwards. The distributions of $\bar{\Delta p}$ and $w$ obtained from the reference simulation have a Gaussian shape, as shown in Fig. 7 for the reference simulation. This is consistent with the distribution of $\Delta\sigma$ which had a Gaussian part due to plastic bands. Therefore, distributions of $\bar{\Delta p}$ and $w$ are also modeled with Gaussian distribution, as described in Eq. (42). The mean values of $\bar{\Delta p}$ and $w$ obtained by this method are respectively denoted $\mu_d^p$ and $\mu_d^w$.

## 4. Assessment of the model

In this section, we focus on the behavior of plastic bands. The effects of smaller-scale plastic events on elastic properties are analyzed in the subsequent section. If not explicitly stated otherwise, the parameters used in the simulations will be the same as those used in the preceding sections and are provided in Table 1. Complete datasets discussed after can be found at 10.5281/zenodo.14266823.

<table><caption>Table 1
Reference material and simulation parameters.</caption><thead><tr><td>Parameter</td><td>Value</td></tr></thead><tbody><tr><td>Initial yield stress σ<sub>ys</sub></td><td>100 MPa</td></tr><tr><td>Linear hardening rate H</td><td>10 GPa</td></tr><tr><td>Young's modulus E</td><td>200 GPa</td></tr><tr><td>Poisson's ratio ν</td><td>0.3</td></tr><tr><td>Plastic threshold Δp<sub>min</sub></td><td>2 × 10<sup>-4</sup></td></tr><tr><td>Prescribed deformation increment Δε<sub>ex</sub> = Δu/L</td><td>3 × 10<sup>-6</sup></td></tr><tr><td>Ratio mesh size over specimen length h/L</td><td>0.5%</td></tr></tbody></table>

### 4.1. Mesh sensitivity

The effect of meshing, with some examples displayed in Fig. 3, is studied in this section. First, the effect of mesh size is tested on linear tetrahedral elements. Then, the effect of mesh type (tetrahedron vs hexahedron), and the angles of hexahedra with respect to the tensile direction are analyzed.

Various ratios of mesh size $h$ over specimen length $L$ are tested. Specimen gauge length and heads have the same element size. The same conditions are used as described in Section 3.3. The tensile curves are given in Fig. 8(a). The effective yield strength for $h/L = 2.5\%$ is 134 MPa. It is superior to the initial yield strength $\sigma_{ys}$ because in time-discontinuous plasticity the yield surface does not define the stress at which plasticity begins (see Fig. 1). The effective yield strength decreases to 132 MPa for $h/L = 1.5\%$, and finally converges at 127 MPa for meshes finer than 0.5%. The average flow stress during the first plateau is 117 MPa and is mesh independent. The serrations are
<figure><img src="image_11.png" /><figcaption>**Fig. 8.** Effect of varying the ratio of mesh size over specimen length *h/L* on (a) the tensile curves, (b) the statistics of stress drops, (c) the statistics of mean plastic strain and (d) the statistics of band widths. In (b), (c) and (d), the dots are representative of the probability density function determined from the simulation and the solid line is the modeled distribution. The color code is indicated in (a). In (a), the tensile curves are shifted in the inset for visibility purposes. In (b), the distributions are indicated in linear scale in the inset.</figcaption></figure>

more varied in intensities for coarser meshes, while they become more periodic for meshes finer than 0.5%.

The PDF of stress drops $\Delta\sigma$ extracted from the simulations are given in Fig. 8(b) in logarithmic and linear scales. The mean value of stress drops $\mu_d^\sigma$ linked with plastic bands are indicated in Fig. 8(d). The modeled distribution are also shown. Looking at the statistics of stress drops reveals that the distribution of $\Delta\sigma$ converges for a mesh size below 1/100 of the specimen size, with the difference in the average stress drop $\mu_d^\sigma$ being about 1 MPa, which is smaller than the standard deviation of the distribution (≈2.6 MPa for purple distribution).

In Fig. 8(c) the distributions of mean plastic strain $\Delta p$ of bands obtained in the five simulations are indicated, alongside the average values of these distributions. All values are normalized by $\Delta p_{min}$. The band strain tends to be smaller for larger mesh elements. The average value $\mu_d^p$ converges as the mesh size decreases, approaching 5 times $\Delta p_{min}$. In Fig. 8(d) the distributions of band width $w$ obtained in the five simulations are indicated, alongside the average values of these distributions. All values are normalized by the length of the specimen length $L$. The band width is larger for coarser meshes. The average value $\mu_d^w$ converges as the mesh size decreases, approaching 4% of the specimen length $L$ (the standard deviation being 1% of $L$ in comparison). Overall, this analysis demonstrates that all relevant distributions (not just the average values) converge for mesh sizes below 0.5% of the total specimen length.

The results of the simulations with hexahedral elements are given in Fig. 9 with the same format as in Fig. 8. The ratio $h/L$ is set to 1/400 for hexahedral meshes, with $h$ defined in Fig. 3. Four orientations $\theta$ of hexahedra have been tested, with angles with respect to the tensile axis ranging from 45° to 90° ($\theta$ indicated in Fig. 3). The angle 54.74° has been selected since it is the orientation observed in isotropic elastoplastic materials of plastic strain localization bands (Mazière and Forest, 2015). The tensile curves of the 4 simulations (Fig. 9(a)) are in good qualitative agreement. Looking at the serrations, they appear more regular when the mesh is oriented at 90° or 60° than for smaller angles.

The stress drop distributions (Fig. 9(b)) reveals the same division, with 90° and 60° Gaussian distribution centered on 9.6 MPa and 54.74° and 45° distributions centered on 5.7 MPa. The same division can be observed in the distributions of band width in Fig. 9(d), with overall smaller bands for 54.74° and 45° orientations. However, the mean plastic strain distributions reveals that while 45° mesh phenomenology is very similar to 60° and 90° meshes (Gaussian distribution centered around $\mu_d^p = 5\Delta p_{min}$), the 54.74° distribution is significantly different. In this case, $\Delta p$ distribution (green dots in Fig. 9(c)) contains two Gaussian-like parts, one centered on $4.5\Delta p_{min}$ and the other on $7.7\Delta p_{min}$. The mean of the 54.74° distribution (calculated on simulation data, not on the Gaussian model) is $6.0\Delta p_{min}$, which is significantly higher than other distributions (standard deviation ≈ 0.6$\Delta p_{min}$).

The effect of mesh orientation on overall specimen response can be understood similarly to the effect of mesh orientation on the simulations of the Lüders phenomenon studied by Mazière and Forest (2015), within the framework of a classical continuum plasticity model and using a softening material. In this work, it is shown that meshes oriented at 54.74° allowed for discontinuities in the strain tensor (i.e. infinitesimal localization bands), which were not observed for meshes oriented at 90° (where band fronts were more diffuse). This is a well-known feature of Lagrange linear and quadratic interpolation functions which allow for strain discontinuities along edges but not inside the elements. Spurious effects arise when elements are crossed by discontinuity lines such as strain bands. This could be improved using Galerkin discontinuous types of elements. In our time-discontinuous model, similar features are observed, although bands with a width of one mesh element never appear while they were observed in Mazière and Forest (2015) for Lüders modeling. A mesh that permits perfect plastic localization alters the statistics of plastic strain, stress drops and band widths. Meshes with a 45° orientation also seem to facilitate the accommodation of more localized bands, but to a lesser extent. Fortunately, when the mesh is more randomly oriented with respect to the ideal localization band orientation, the statistics of the bands seem to converge, and the convergence values (mean and standard deviation)
<figure><img src="image_12.png" /><figcaption>**Fig. 9.** Effect of varying the orientation θ of hexahedral meshes on (a) the tensile curves, (b) the statistics of stress drops, (c) the statistics of mean plastic strain and (d) the statistics of band widths. In (b), (c) and (d), the dots are representative of the probability density function determined from the simulation and the solid line is the modeled distribution. The color code is indicated in (a). In (a), the tensile curves are shifted in the inset for visibility purposes. In (b), the distributions are indicated in linear scale in the inset.</figcaption></figure>

match those found for tetrahedral meshes. Therefore, for the model to yield a mesh-independent response, meshes must be sufficiently refined and not oriented in a way that favors the formation of localized bands.

## 4.2. Effect of the strain increment $\Delta\epsilon_{xx}$

In the simulation, axial displacement is imposed at the end of the specimens. This displacement can be linked with a increment of axial strain $\Delta\epsilon_{xx}$ in the gauge length, which is constant during the whole simulation. In this section, we examine the convergence of the model with respect to $\Delta\epsilon_{xx}$. The mesh consists of tetrahedra with a $h/L$ ratio of 5%.

In Fig. 10(a), the tensile curves are shown for five simulations. When $\Delta\epsilon_{xx}$ is too high, no pronounced serration are observed on the tensile curve. As it decreases, clear serrations appear (yellow and green). As $\Delta\epsilon_{xx}$ decreases further, serrations tend to become more regular.

The statistics of stress drops, mean plastic strains and band widths are given for the five simulations respectively in Fig. 10(b), (c) and (d). For the largest $\Delta\epsilon_{xx}$, the distribution of stress drops is not displayed as, in the absence of serration on tensile curve, they have little meaning. The distributions of band mean plastic strains and sizes for $\Delta\epsilon_{xx} = 0.5\Delta p_{min}$ reveal that bands appear in this case, despite the absence of serration on the tensile curves, and are very diverse in size and intensity. This apparent contradiction is explained by the fact that while plastic bands appear and tend to decrease the stress of the system, the axial stress increment provided to the system $E\Delta\epsilon_{xx}$ at each numerical step is high enough to compensate the drop.

When $\Delta\epsilon_{xx}$ is decreased, the three distributions converge quickly, not only in terms of average values but also in shape. Therefore, the model response converges when the displacement increment imposed at the boundary is small enough. The criterion is for $\Delta\epsilon_{xx}$ to be below $0.05\Delta p_{min}$. This analysis demonstrates that the presented model converges with respect to the incremental change in boundary conditions.

## 4.3. Effect of constitutive parameters

The effects of constitutive parameters are studied in this section on dogbone specimens with tetrahedral meshes. For conciseness, only mean values of stress drop $\mu_d^\sigma$, band plastic strain $\mu_d^p$ and band size $\mu_d^w$ will be analyzed. Fig. 11(a), (b) and (c) show the effects of $\Delta p_{min}$, varying from $5 \times 10^{-5}$ to $8 \times 10^{-4}$ on the evolution of $\mu_d^\sigma$, $\mu_d^p$ and $\mu_d^w$ determined from simulations (black dots). In this range, both $\mu_d^\sigma$ and $\mu_d^p$ increase linearly with $\Delta p_{min}$, while the value of $\mu_d^w$ changes marginally (variation well below the standard deviation of each distribution, indicated by a red lines in Fig. 11(c)).

The linear work-hardening modulus $H$ was varied from 0, to 10,000 MPa. The corresponding results for the evolution of $\mu_d^\sigma$, $\mu_d^p$ and $\mu_d^w$ are given in Fig. 11(d), (e) and (f). For $H$ ranging from 0 to 1000 MPa, no characteristic of the bands evolves. Between 1000 and 10,000 MPa, both $\mu_d^\sigma$ and $\mu_d^p$ decrease by a factor 1.7, while $\mu_d^w$ decreases only by a factor of 1.1.

Finally, the last model parameter tested is $v$. It was varied from 0 to 0.49. Evolution of $\mu_d^\sigma$, $\mu_d^p$ and $\mu_d^w$ as $v$ is varied are given in Fig. 11(g), (h) and (i), for two given values of $H$: 10,000 MPa in black, and 1000 MPa in blue. For $H = 10,000$ MPa, varying $v$ from 0 to 0.49 decreases $\mu_d^\sigma$ by a factor of 1.2, $\mu_d^p$ by a factor of 1.8 and increases $\mu_d^w$ by a factor of 1.5. For $H = 1000$ MPa, varying $v$ from 0 to 0.49 decreases $\mu_d^\sigma$ by a factor of 1.4, $\mu_d^p$ by a factor of 2.3 and increases $\mu_d^w$ by a factor of 1.5. Therefore, it is clear that when $H$ is larger, the effect of Poisson's ratio on the band plastic strain and stress drop intensity is less significant. It is also apparent that $v$ has a significant impact on the band width, and that the effect of $H$, while small, is not insignificant, as a constant factor of 1.1 on $\mu_d^w$ is almost always found when comparing results obtained for $H = 10,000$ MPa and $H = 1000$ MPa.

While no explanation has been found yet for the dependence of band width with each parameter, it is possible to understand the evolutions of band strain and stress drop. Assuming that the fraction $\phi$ of the plastic band in the specimen is known and remains constant for each
<figure><img src="image_13.png" /><figcaption>Fig. 10. Effect of varying the ratio of axial strain increment over plastic threshold $\Delta\epsilon_{xx}/\Delta p_{min}$ on (a) the tensile curves, (b) the statistics of stress drops, (c) the statistics of mean plastic strain and (d) the statistics of band widths. In (b), (c) and (d), the dots are representative of the probability density function determined from the simulation and the solid line is the modeled distribution. The color code is indicated in (a). In (a), the tensile curves are shifted in the inset for visibility purposes. In (b), the distributions are indicated in linear scale in the inset.</figcaption></figure>

plastic burst, it is possible to predict the value of the plastic strain in the band:

$$ \|p\| = \frac{3\mu + H}{\phi E + H} \Delta p_{min} \qquad (43) $$

The Eq. (43), analogous to Eqs. (27) and (39), is derived in the same way by replacing $\|\epsilon^P\|$ by $\phi\|\epsilon^P\|$ in Eq. (33). This formula is valid if the stress state in the specimen is tensile and strictly uniform. For a dogbone specimen, the value of $\phi$ for a given band is close to the ratio between the band width and the specimen length $w/L$. As the value of $\phi$ remains a priori unknown in the model, a constant value $\phi^{ref} = 4.5\%$ is chosen for all future calculations, close to the value of $\mu_d^w/L$ found in the simulation.

Using the previous Eqs. (28) and (43) ($\|\sigma_{vM}\| = -E\|p\|$), and assuming a value for $\phi = \phi^{ref}$, it is possible to give an estimate $\mu_d^{\sigma}$ and $\mu_{theo}^{p}$ for all simulations (green dashed curves in Fig. 11). A multiplicative factor between 0.35 and 0.4 has been applied to the obtained values to be closer to the simulations results (as indicated in Fig. 11(b), (e), (h) and (k)). This multiplicative factor can be explained by geometric consideration, making the stress state far from uniform in the dogbone specimen.

Overall, the proportionality factor between $\Delta p_{min}$ and both $\mu_d^{\sigma}$ and $\mu_d^{p}$ is well explained, alongside the slightly more complex evolution of those quantities. The small gap between simulation data and modeled evolution can be explained by the non constant value of $\phi$ in simulation, especially as $v$ change. All model parameters affect the plastic strain in bands and the consecutive stress drop in a theoretically well-understood manner. Most model parameters do not affect the band width, but some parameters, especially $v$, do, for reasons that are not fully understood yet.

## 4.4. Effect of gauge length

To study the impact of geometry on the model response, the gauge length $L$ has been modified from a value $L^{ref}$ (which was the one used

before) to a value of $2L^{ref}$. Fig. 11(j), (k) and (l) display the effects of $L$ on the evolution of $\mu_d^{\sigma}$, $\mu_d^{p}$ and $\mu_d^w$ determined from simulations. The average value of band widths $\mu_d^w$ evolves proportionally to the gauge length, with a factor of proportionality of 4%, which is the average fraction $\phi$ of the specimen in the plastic band in the reference state (see Fig. 10(d)). This means that the ratio $\phi = \mu_d^w/L$ remains constant when $L$ is modified. On the other hand, changing the gauge length has no effect on $\mu_d^{\sigma}$ and $\mu_d^{p}$, which was expected from Eqs. (28) and (43), since $\phi$ remains constant as $L$ changes and other parameters do not depend on $L$.

Simulation with dogbone geometries transformed under homothetic dilatation have been performed. Dilatation factors of 2 and 4 have been tested (results not shown here). Again, all dimensions of the bands evolve proportionally with the corresponding dimensions of the specimen (band width with specimen length, band height with specimen width and band thickness with specimen thickness). Other features of the bands remain unaffected by homothetic transformation (band mean plastic strain and stress drop intensity). The model is thus homothetic. This result is consistent with the absence of an internal length in the model. Only a plastic threshold has been introduced, which is not linked with any length.

# 5. Discussion

## 5.1. Statistics of stress drop amplitude and frequency

In the previous sections, it has been demonstrated that the plastic events occur following two distributions: a Gaussian part (big plastic events) and a power law part (small plastic events). The properties of big plastic events alongside the effect of model parameters on them have been extensively analyzed in the previous sections. Here, we analyze the truncated power law distribution.

In Figs. 8, 9 and 10(b) the distributions of stress drops in logarithmic scales are represented for the simulations testing the effects of mesh
<figure><img src="image_14.png" /><figcaption>**Fig. 11.** Effect of plastic threshold $\Delta p_{\min}$ on (a) $\mu_d^\sigma$, (b) $\mu_d^p$ and (c) $\mu_d^{w}$. Effect of work-hardening modulus $H$ on (d) $\mu_d^\sigma$, (e) $\mu_d^p$ and (f) $\mu_d^{w}$. Effect of Poisson's ratio $\nu$ on (g) $\mu_d^\sigma$, (h) $\mu_d^\sigma$ and (i) $\mu_d^w$. Effect of gauge length $L$ on (j) $\mu_d^\sigma$, (k) $\mu_d^\sigma$ and (l) $\mu_d^w$. In black are indicated the results from simulations. For $\nu$, two different values of $H$ were tested. The values found for $H = 1000$ MPa is indicated in blue, and the values found for $H = 10,000$ MPa is indicated in black. In green dashed line are the theoretical evolutions found using Eqs. (28) and (43) using the reference parameters, as written for each figure.</figcaption></figure>

size, orientation, and prescribed displacement at boundaries, respectively. The simulation results (dots) are displayed with the modeled truncated power law (solid line). When the mesh size is too big, or the axial strain $\Delta\epsilon_{xx}$ too large, deviation from power law is observed for stress drop events smaller than $10^{-1}$ MPa. However, when $h/L$ is lower than 0.5%, and $\Delta\epsilon_{xx}/\Delta p_{\min}$ lower than 1%, then the stress drop distribution follow indeed a power law for its smaller events. For any orientation of the hexahedral mesh, a good match is found between the stress drop distribution and the power law.

The value of the exponent of the power law $\alpha$ is presented for all simulations in Fig. 12. Fig. 12(a), (b) and (c) show the values of

$\alpha$ for non-model parameters ($h/L$, angle $\theta$ of mesh and $\Delta\epsilon_{xx}/\Delta p_{\min}$). The values of $\alpha$ range from 1.0 to 1.5. More precisely, one group of simulations lies around 1.0, while a second group lies between 1.3 and 1.5. The first group corresponds to simulations with a high ratio of $h/L$ and $\Delta\epsilon_{xx}/\Delta p_{\min}$, where the simulation distribution did not follow a power law, or with a mesh consisting of hexahedra oriented at $\theta = 54.74^{\circ}$, where the plastic events showed atypical and unrepresentative behavior. The second group consists of simulations with finer meshes and prescribed displacements. Simulations with hexahedral meshes that are not aligned too favorably with the ideal orientation of the plastic bands have the same exponent value as simulations with fine enough
<figure><img src="image_15.png" /><figcaption>**Fig. 12.** Evolution of power law exponent $\alpha$ as a function of (a) $h/L$, (b) $\theta$ (as defined for hexahedral elements), (c) $\Delta\epsilon_{xx}/\Delta p_{min}$, (d) $L/L^{ref}$, (e) $\Delta p_{min}/\Delta p_{min}^{ref}$, (f) $H$, (g) $v$ with $H = 10,000$ MPa and (h) $v$ with $H = 1000$ MPa.</figcaption></figure>

meshes and small enough displacement steps. It is interesting to note that the power law exponent converges when the mesh is refined and the value of $\Delta\epsilon_{xx}$ is small enough. It is also noteworthy that the values of the simulation parameters required for convergence are similar to those obtained for the convergence of large stress drop distributions. This implies that small plastic events and large ones are closely related.

By varying the model parameters, $\alpha$ varies between 1.3 and 1.6 as shown in Fig. 12(e)–(h). No clear trend can be seen for the effect of $\Delta p_{min}$ on $\alpha$, unlike that found for the large stress drop characteristics. When $\Delta p_{min}$ is small or large, $\alpha$ is close to 1.45. For $H$, values below 1000 MPa lead to $\alpha$ being close to 1.55, while increasing $H$ to 10,000 MPa leads to a decreased value of $\alpha$ = 1.4. When $H$ equals 1000 MPa, having $v$ less than 0.3 leads to a value $\alpha$ close to 1.55, while it decreases to 1.3 when $v$ = 0.49. When $H$ is 10,000 MPa, varying $v$ does not affect $\alpha$, which remains between 1.3 and 1.4. Finally, if the value of the sample length is increased from $L_{ref}$ to $2L_{ref}$, the value of $\alpha$ remains around 1.4, which means that the length has no direct effect on $\alpha$.

In the literature on the plasticity of polycrystalline materials, a power law exponent $\alpha$ between 1 and 1.6 is typically observed for tensile deformation, which is consistent with the values predicted by this model (Csikor et al., 2007; Perchikov and Truskinovsky, 2024; Weiss et al., 2007; Brown, 2012; Patinet et al., 2011; Dimiduk et al., 2006). For example, Borasi et al. (2023) found an exponent $\alpha$ of 1.5 for the tensile deformation of microcast silver crystals with dogbone geometry. This range of exponents is also consistent with those found for metals affected by the PLC effect (Lebyodkin et al., 2000; Ananthakrishna, 2005). Although the system-sized plastic events in the model deviate from the experimentally observed behavior due to their Gaussian nature, the small-scale plastic events exhibit a phenomenology similar to that seen in experiments and physical models (MD, DDD, MTM). Thus, the model demonstrates an ability to capture certain key aspects of the underlying physics behind intermittent plasticity.

## 5.2. Spatiotemporal analysis of strain bursts

The spatiotemporal evolution of plastic strain increment on the $O_x$ axis is given in Fig. 13 for the reference simulation, with a zoom at the start and at the end of the simulation. In those maps, horizontal lines are plastic bands happening in one numerical step.

The first series of bands, from $t = 0$ to 0.02, has a “pseudo-propagative” character in the sense that the nth band will always appear next to the (n-1)-th band. In some simulations, this propagation goes from one end of the gauge to the other. In the current simulation,

a second series of “pseudo-propagative” bands appear at the left end of the gauge length instead of only one, at $t = 0.021$. Then, a third series of correlated bands appear at $t = 0.04$, as a continuation of the first series. The behavior of band nucleation is phenomenologically analogous to type B PLC (Yilmaz, 2011). However, this first propagation occurs at a stress level whose average value remains constant, which is more similar to a Lüders phenomenon. The reason for this correlated nucleation of bands is the overstress appearing at the boundary of each band, the work-hardening rate $H > 0$ and the absence of preliminary defect on the gauge length (no concentration of stress prior to the plastic deformation). When $H$ is less than 100 MPa, bands appear but are confined to the area between the gauge length and the head of the specimen where overstress exists due to the dogbone geometry.

Once the entire tensile specimen has been plastically deformed for the first time ($t = 0.06$), elements with significant local stress concentrations are now randomly distributed throughout the gauge length of the specimen. Between $t = 0.06$ and 0.15, plastic bands appear more random and less correlated. Between $t = 0.85$ and 1 (13 (c)), bands are nucleating purely randomly on the gauge length, and no series of pseudo-propagative bands are found. This regime of nucleation is comparable to type C PLC (Yilmaz, 2011).

The distances between consecutive bands (distance between their center) have been measured during the experiment, and their distribution normalized by the specimen length is shown in Fig. 13(d) in logarithmic scale. It is found that they are distributed following a power law with an exponent $\alpha = 1.0$, which is also typical of type C PLC (Jiang et al., 2005).

The interest of the presented model is to spontaneously induce macroscopic deformation bands within the specimen. These bands, with finite width, can be interpreted as shear bands caused by an avalanche of dislocations. Ultimately, while remaining within an isotropic elasto-plastic framework implementable in finite element codes, the model mimics the same phenomenology as that due to the self-organization of dislocations among themselves in real materials, which is moreover an emerging phenomenon in the model.

## 5.3. Application to structural computations

To evaluate the model under complex multiaxial conditions, simulations were performed on a 3D holed plate geometry. Both the classical continuous plasticity (CCP) and the newly developed time-discontinuous plasticity (TDP) have been tested. Fig. 15 shows the tensile curves of both simulations, the distribution of stress drops from the TDP simulation and the cumulative plastic strain maps for the
<figure><img src="image_16.png" /><figcaption>Fig. 13. (a) Spatiotemporal evolution of plastic strain increment during the reference simulation. X-axis is the relative position on the gauge length (0.15 is the left end of the gauge length, 0.85 its right end). Y-axis is the normalized simulation time. (b) Zoom of (a) at the start of the simulation. (c) Zoom of (a) at the end of the simulation. (d) Distribution of distances between consecutive bands (blue) and the corresponding modeled power law (red).</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>Fig. 14. (a) Tensile curves obtained from the simulations of a holed plate geometry, using either the classical continuous plasticity (black) or the time-discontinuous plasticity (red). (b) Distribution of stress drops from the time-discontinuous simulation (black dots) and the corresponding power law (red solid line). (c) Cumulative plastic strain maps for the continuous plasticity simulation at $ε_1$, (d) $ε_4$ and (e) $ε_5$.</figcaption></figure>

CCP simulation at 3 given deformation step. The boundary conditions consists of mixed conditions on left and right faces (only axial displacement imposed), the other faces being free. Rigid body motion is suitably fixed. Model parameters used are the same as for the reference dogbone geometry simulation (see Table 1). Meshes were very refined around the hole ($h/L = 0.25\%$) and coarser around the edge of the

specimen ($h/L = 1\%$). The simulation was done in 3D and there were 4 tetrahedra in the thickness.

The CCP plate begins to yield macroscopically at $ε_1 = 0.06\%$, when the average stress within the plate becomes greater than $σ_{ys} 100$ MPa, as seen on the corresponding tensile curve at the slope break. Before $ε_1$ only the area around the hole will yield plastically. After $ε_1$ the whole plate deforms plastically. The cumulative plastic strain maps
show that the plasticity is concentrated around the hole at $ε_1$, $ε_4$ and $ε_5$, forming the x-shape expected for this geometry as it corresponds to the maximum values for von Mises stress.

The tensile curve of the TDP plate shows some serrations, as found for the dogbone geometry. The effective yield stress is 125 MPa, which is higher than the CCP tensile curve, as expected from previous analysis on dogbone specimens. The stress drop distribution (Fig. 14(b)) is calculated as indicated in Section 3.4, using the axial stress. Although this value does not have the same meaning as before ($σ_{VM}$ is different from $σ_{xx}$ due to the hole), it is still indicative of plastic activity. The shape of the curves is slightly different, with a first part (below 2.5 MPa) following a power law and the second part following a plateau instead of a Gaussian distribution. This is very similar to the supercritical distributions observed for metallic nanopillars smaller than 500 nm in diameter (Zhang et al., 2020). The first part has been modeled with a truncated power law, defined in Eq. (42). The power law exponent $α$ found is 1.36, which is very close to the values found in the dogbone geometry. This suggests that the rate of nucleation of small plastic events at the boundaries of plastic bands is also geometry, and thus intrinsic to the time-discontinuous model.

The maps of cumulative plastic strain, plastic strain increment and von Mises stress for the time-discontinuous simulations are displayed in Fig. 15. As observed in continuous plasticity for this geometry, stress concentration is found around the hole, resulting in higher cumulative plastic strain fields above and below the hole (step 1). In contrast to classical continuum plasticity, where the strain around the hole appears gradually, the plastic strain in our model appears in a spatially and temporally discontinuous manner.

The first plastic band appearing break the central symmetry of the system (step 2). It is not traversing, with a width of 9 elements. At its edge, von Mises stress concentration is found, which will induce new bands, not necessarily with the same orientation.

The first traversing plastic bands appear in step 3, associated with the first serration on the tensile curve at $ε_3$. It is not contiguous. The average plastic strain of the band is the same as for the equivalent dogbone simulation (≈ 5 $Δp_{min}$). The width of the band increases as it extends away from the central hole. At the boundary of the sample it has a width of 0.04 L, which is the average width of the bands in the corresponding dogbone simulation. The following plastic bands appear next to the previous one due to von Mises stress concentration, with the same orientation (step 4).

Similar to the dogbone geometry, plastic activity consists of small plastic events (occurring in some elements, not always contiguous) and large plastic events. The large plastic events can manifest as traversing plastic bands or non-traversing bands around the hole. At the end of the simulation (step 5), the entire sample is covered by plastic bands with the same orientation of ≈ 54.74°. This phenomenon can be understood by the shape of the residual von Mises stress, which follows the same direction and favors this direction.

Finally, the comparison of the cumulative plastic strain map at $ε_5$ of the CCP plate (Fig. 14(e)) and the TDP plate (Fig. 15) reveals a key characteristic of the model. While both plates exhibit large-scale plastic deformation, the map for the time-discontinuous model reveals significant heterogeneity in plastic strain distribution, whereas the CCP map shows a more homogeneous deformation away from the hole. Therefore, this demonstrates that the model successfully achieves the critical objective of obtaining intense localization of plastic deformation into discrete bands, as found in intermittent plasticity (Charpagne et al., 2021). Unlike previous models, this result is achieved without the introduction of probabilistic elements into the plasticity model itself, thus eliminating the need for prior knowledge of the laws governing plastic flow correlations (Gélébart, 2021; Marano et al., 2019).

## 6. Conclusions

A new time-discontinuous plasticity model was developed in this work and implemented independently in the FEniCSX (Scroggs et al., 2022) and Zset (Besson and Foerch, 1997) solvers. By introducing a single new parameter to conventional J2-plasticity, called the plastic threshold $Δp_{min}$, the model is able to capture the spatial and temporal localization of plastic flow.

Simulations under homogeneous deformation conditions demonstrate that plastic deformation occurs through abrupt increases in plastic strain, comparable to dislocation avalanches, and these events are always associated with sharp stress drops. The model is interpreted through the existence of two yield surfaces: the upper surface controls the initiation of plasticity, while the lower surface defines the stress reached after plastic relaxation. The results of the new model under homogeneous fields are fully understood under any triaxiality.

When the geometry of the specimen is complex and induces stress concentration, plastic strain bands spontaneously emerge in a random manner, without introducing any stochasticity as an input to the model. Stochasticity of plastic events is an outcome of the model, not an input in contrast to many existing models in this field. Temporal intermittency of plastic activity and consecutive stress drops obey a power law for small events and a Gaussian distribution for large events. Spatial intermittency of plastic activity follows a power law. The effect of model and simulation parameters, such as $Δp_{min}$, elastic constants, work-hardening modulus and mesh size, have been extensively studied. Notably the band widths are proportional to the length of the tensile specimen, and the plastic strain carried by the band is proportional to $Δp_{min}$. The model is homothetic in nature as it contains no internal length. It has been proven to be usable under any geometry, no matter its complexity, while remaining efficient in terms of calculation time. One limitation of the model is that, to obtain serrations and reliable statistics on stress drops and plastic events, the applied strain in the specimen must be at least below to 5% of $Δp_{min}$. No spurious mesh dependency is observed in contrast to softening models for localization.

This model developed in the von Mises plasticity framework is the first step towards the development of a time-discontinuous model in crystal plasticity. It is believed that this next step is necessary for the current work to be more physically relevant and thus comparable to experimental data. To control band size, a statistical spatial distribution of $Δp_{min}$ across the specimen should be introduced into the model. The statistic of $Δp_{min}$ could be calibrated using data from MD, MTM, or DDD models, which represent the most physically relevant approaches.

## Funding sources

This work was funded by the French ANR (Agence Nationale de la Recherche) under the MESOCRYSP project (ANR-21-CE08-0030).

## CRediT authorship contribution statement

**M. Lamari:** Writing – original draft, Visualization, Software, Methodology, Investigation, Formal analysis, Conceptualization.
**P. Kerfriden:** Writing – original draft, Supervision, Software, Methodology, Formal analysis, Conceptualization.
**O.U. Salman:** Writing – review & editing, Supervision, Methodology, Formal analysis.
**V. Yastrebov:** Writing – review & editing, Visualization, Supervision, Methodology, Conceptualization.
**K. Ammar:** Validation, Software.
**S. Forest:** Writing – review & editing, Supervision, Project administration, Methodology, Funding acquisition, Conceptualization.

## Declaration of generative AI in scientific writing

During the preparation of this work the author(s) used ChatGPT 3.5 and DeepL in order to improve language and readability. After using this tool/service, the author(s) reviewed and edited the content as needed and take(s) full responsibility for the content of the published article.
<figure><img src="image_18.png" /><figcaption>**Fig. 15.** Map of cumulative plastic strain, plastic strain increment and von Mises stress at five different levels of deformation of the time-discontinuous simulation. The strain levels are shown in Fig. 14. At stage 1, the fields are close to those obtained with classical plasticity. At stage 2 the first significant band appears. At stage 3, the first traversing band appears. At stage 4, several traversing bands have appeared. Stage 5 is the state at the end of the simulation.</figcaption></figure>

## Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

## Appendix A. Finite element solver and time discretization of the time-discontinuous elasto-plasticity model

The time interval $\tau = [0, T]$ is divided into $N$ time steps. The equation of continuum mechanics is enforced at discrete times $\bar{T} = \{t_0, t_1, \dots, t_N\}$, using the finite element method. In the following, $n$ will be an integer between 0 and $N$.
**Finite element solver.** The domain $\Omega$ is discretized into the set $\mathcal{E}$ of finite elements such that $\Omega \approx \Omega^h = \bigcup_{v \in \mathcal{E}} \Omega^{(v)}$. The space of linear functions over each element $\Omega^{(v)}$ is denoted $\mathcal{P}^1(\Omega^{(v)})$. The finite element space $\mathcal{V}_{n+1}$ in which the displacement field $u_{n+1}$ at time $t_{n+1}$ will be looked for is defined by:

$$ \mathcal{V}_{n+1} = \{u \in C^0(\Omega) | \forall v \in \mathcal{E}, u_{|\Omega^{(v)}} \in \mathcal{P}^1(\Omega^{(v)}), u = u_d|_{t=t_{n+1}} \text{ over } \partial_u \Omega \} \quad (44) $$

and the associated vector space of displacement fields that vanish over $\partial_u \Omega$:

$$ \mathcal{V}_0 = \{u \in C^0(\Omega) | \forall v \in \mathcal{E}, u_{|\Omega^{(v)}} \in \mathcal{P}^1(\Omega^{(v)}), u = 0 \text{ over } \partial_u \Omega\} \quad (45) $$

The weak form of the time-discontinuous elasto-plasticity equations at time $t_{n+1} \in \tilde{T}$, discretized by the finite element method, is now introduced. For any virtual displacement $v \in \mathcal{V}_0$, the displacement field $u_{n+1} \in \mathcal{V}_{n+1}$ is searched such that the finite element stress field $\sigma_{n+1}$ at time $t_{n+1}$ satisfies (in absence of volumetric forces and with trivial Neumann boundary conditions):

$$ R_{n+1}(v; u_{n+1}) = \int_{\Omega} \nabla_s v \cdot \sigma (\nabla_s u_{n+1}; \varepsilon_n^p, p_n) \, d\Omega = 0 \quad (46) $$

where $\nabla_s u_{n+1} = \varepsilon_{n+1}$ is the finite element strain field at time $t_{n+1}$. The dependency of the stress on the local history of the strain tensor is introduced via the functional dependency of $\sigma$ on the internal variables $\varepsilon_n^p$ and $p_n$ computed at time step $t^n$.

**Time discretization of the constitutive relation.** The time-discontinuous constitutive equations introduced in the previous section need to be appropriately discretized in time. More precisely, Eqs. (13), (14), (15) and (16) must be replaced by their time-discrete counterparts so that $\sigma (\nabla_s u_{n+1}; \varepsilon_n^p, p_n)$ can be evaluated numerically. A fully implicit time-stepping scheme is used. Eq. (6) is enforced at time $t_{n+1}$:

$$ \sigma_{n+1} = \mathbb{C} : (\varepsilon_{n+1} - \varepsilon_{n+1}^p) \quad (47) $$

The normality rule Eq. (13) is discretized as follows:

$$ \varepsilon_{n+1}^p - \varepsilon_n^p = (p_{n+1} - p_n) \mathbf{n}_{n+1}^{-} \quad (48) $$

Then, introducing the set $\mathcal{P}_{n+1}$ of cumulative plastic strain increments $\Delta p^*$ satisfying

$$ \Delta p^* f_{n+1}^* = 0 \quad (49) $$

and

$$ \Delta p^* \langle \Delta p_{\min} - \Delta p^* \rangle = 0 \quad (50) $$

the maximum value in set $\mathcal{P}_{n+1}$ is chosen as the solution of the constitutive update, i.e.

$$ p_{n+1} = \max_{\Delta p^* \in \mathcal{P}_{n+1}} \Delta p^* + p_n \quad (51) $$

In the previous set of equations, $\mathbf{n}_{n+1}^{-}$ denotes the normal to the yield surface at time $t_{n+1}$, assuming that no plastic jump took place between $t_n$ and $t_{n+1}$. Mathematically, this means that $\mathbf{n}_{n+1}^{-} = \frac{\partial f}{\partial \sigma}|_{\sigma_{n+1}=C:(\varepsilon_{n+1}-\varepsilon_n^p), p_{n+1}=p_n}$. Notation $f_{n+1}^* \in \mathcal{P}_{n+1}$ stands for the value of yield function $f$ at the end of the time step, i.e. $f_{n+1}^* = f(\mathbb{C} : (\varepsilon_{n+1} - \varepsilon_n^p - \Delta p^* \mathbf{n}_{n+1}^{-}; p_n + \Delta p^*))$. Knowing $\varepsilon_{n+1}$, the equations of the time-discrete constitutive law are solved algorithmically, and pointwise, to obtain the values of $p_{n+1}$, $\varepsilon_{n+1}^p$, and $\sigma_{n+1}$. This local algorithmic procedure (the so-called constitutive update) will be detailed in the next section.

**Global Newton solver.** The Newton-Raphson solver used to solve the nonlinear finite element system of equations at time $t_{n+1}$ requires the linearization of residual (46). The Gateau-derivative of $R_{n+1}$ in direction $\delta u \in \mathcal{V}_0$ is the following bilinear form (linear both in $\delta u$ and in $v$):

$$ \delta R_{n+1} (\delta u, v; u_{n+1}) = \int_{\Omega} \nabla_s v \cdot \left. \frac{\partial \sigma_{n+1}}{\partial \varepsilon_{n+1}} \right|_{\varepsilon_{n+1} = \nabla_s u_{n+1}} : \nabla_s \delta u \, d\Omega \quad (52) $$

The discrete equations of time-discontinuous elasto-plasticity may be formally solved using a Newton algorithm. At iteration $k$ of the Newton procedure, we look for a finite element displacement correction $\Delta u_{n+1,k+1} = u_{n+1,k+1} - u_{n+1,k} \in \mathcal{V}_0$ that satisfies

$$ \delta R_{n+1} (\Delta u_{n+1,k+1}, v; u_{n+1,k}) = -R_{n+1} (v; u_{n+1,k}) \quad (53) $$

for all virtual displacements $v \in \mathcal{V}_0$. Newton initialization $u_{n+1,0}$ must be chosen such that it satisfies the Dirichlet conditions at time $t_{n+1}$, i.e. $u_{n+1,0} \in \mathcal{V}_{n+1}$.

When solving the finite element system of equations at time $t_{n+1}$, the stress field $\sigma_{n+1}$ needs to be evaluated at the quadrature points of the finite element mesh. Internal variables from $t_n$, namely $\varepsilon_n^p$ and $p_n$, are stored at quadrature points, following the algorithmic procedure detailed in Section 2.4.

## Appendix B. Supplementary data

Supplementary material related to this article can be found online at https://doi.org/10.1016/j.ijsolstr.2024.113171.

### Data availability

The data that support the findings of this study are openly available in the Zenodo repository titled 'Time-discontinuous plasticity, data and videos' at https://zenodo.org/records/14266823, reference number 10. 5281/zenodo.14266823. The code used for FENiCSX simulations can be found at https://github.com/Mathias-Lamari/Time-discontinuous-plasticity.

### References

Alstott, J., Bullmore, E., Plenz, D., 2014. Powerlaw: A Python package for analysis of heavy-tailed distributions. In: Rapallo, F. (Ed.), PLoS ONE (ISSN: 1932-6203) 9 (1), e85777. http://dx.doi.org/10.1371/journal.pone.0085777.

Ananthakrishna, G., 2005. Spatio-temporal features of the Portevin-Le Chatelier effect. Mater. Sci. Eng. A (ISSN: 0921-5093) 400-401, 210-213. http://dx.doi.org/10. 1016/j.msea.2005.03.037, Dislocations 2004, URL https://www.sciencedirect.com/ science/article/pii/S0921509305002650.

Baggio, R., Salman, O.U., Truskinovsky, L., 2023a. Homogeneous nucleation of dislocations as a pattern formation phenomenon. Eur. J. Mech. A Solids (ISSN: 0997-7538) 99, 104897. http://dx.doi.org/10.1016/j.euromechsol.2022.104897, URL https:// www.sciencedirect.com/science/article/pii/S0997753822003278.

Baggio, R., Salman, O.U., Truskinovsky, L., 2023b. Inelastic rotations and pseudoturbulent plastic avalanches in crystals. Phys. Rev. E 107, 005000-1-005000-19.

Besson, J., Cailletaud, G., Chaboche, J.-L., Forest, S., 2009. Non-linear Mechanics of Materials, vol. 167. Springer Science & Business Media.

Besson, J., Foerch, R., 1997. Large scale object-oriented finite element code design. Comput. Methods Appl. Mech. Engrg. 142, 165-187.

Borasi, L., Frasca, S., Charbon, E., Mortensen, A., 2023. The effect of size, orientation and temperature on the deformation of microcast silver crystals. Acta Mater. (ISSN: 1359-6454) 249, 118817. http://dx.doi.org/10.1016/j.actamat.2023.118817, URL https://www.sciencedirect.com/science/article/pii/S1359645423001489.

Brown, L., 2012. Constant intermittent flow of dislocations: central problems in plasticity. Mater. Sci. Technol. 28 (11), 1209-1232. http://dx.doi.org/10.1179/ 174328412X13409726212768.

Charpagne, M.A., Hestroffer, J.M., Polonsky, A.T., Echlin, M.P., Texier, D., Valle, V., Beyerlein, I.J., Pollock, T.M., Stinville, J.C., 2021. Slip localization in Inconel 718: A three-dimensional and statistical perspective. Acta Mater. (ISSN: 1359- 6454) 215, 117037. http://dx.doi.org/10.1016/j.actamat.2021.117037, URL https://www.sciencedirect.com/science/article/pii/S1359645421004171.

Chevalier, J., Brassart, L., Lani, F., Bailly, C., Pardoen, T., Morelle, X.P., 2018. Unveiling the nanoscale heterogeneity controlled deformation of thermosets. J. Mech. Phys. Solids (ISSN: 0022-5096) 121, 432-446. http://dx.doi.org/10.1016/j.jmps.2018.08. 014, URL https://www.sciencedirect.com/science/article/pii/S0022509618303909.

Clauset, A., Shalizi, C.R., Newman, M.E.J., 2009. Power-law distributions in empirical data. SIAM Rev. 51 (4), 661-703. http://dx.doi.org/10.1137/070710111.

Colas, D., Finot, E., Flouriot, S., Forest, S., Mazière, M., Paris, T., 2014. Investigation and modeling of the anomalous yield point phenomenon in pure tantalum. Mater. Sci. Eng. A (ISSN: 0921-5093) 615, 283-295. http://dx.doi.org/ 10.1016/j.msea.2014.07.028, URL https://www.sciencedirect.com/science/article/ pii/S0921509314008934.
Csikor, F.F., Motz, C., Weygand, D., Zaiser, M., Zapperi, S., 2007. Dislocation avalanches, strain bursts, and the problem of plastic forming at the micrometer scale. Science 318 (5848), 251–254. http://dx.doi.org/10.1126/science.1143719, URL https://www.science.org/doi/abs/10.1126/science.1143719.

de Geus, T.W.J., Wyart, M., 2022. Scaling theory for the statistics of slip at frictional interfaces. Phys. Rev. E 106, 065001. http://dx.doi.org/10.1103/PhysRevE.106.065001.

de Souza Neto, E.A., Peric, D., Owen, D.R.J., 2011. Computational Methods for Plasticity: Theory and Applications. John Wiley & Sons.

Dimiduk, D.M., Woodward, C., LeSar, R., Uchic, M.D., 2006. Scale-free intermittent flow in crystal plasticity. Science 312 (5777), 1188–1190, (ISSN: 00368075, 10959203) URL http://www.jstor.org/stable/3846259.

Fernández-Castellanos, D., Roux, S., Patinet, S., 2021. Insights from the quantitative calibration of an elasto-plastic band from a lennard-jones atomic glass. Comptes Rendus. Physique 22, 135–162.

Gélébart, L., 2021. Grain size effects and weakest link theory in 3D crystal plasticity simulations of polycrystals. Comptes Rendus. Physique 22 (S3), 313–330. http://dx.doi.org/10.5802/crphys.53.

Gómez-García, D., Devincre, B., Kubin, L.P., 2006. Dislocation patterns and the similitude principle: 2.5D mesoscale simulations. Phys. Rev. Lett. 96, 125503–1–125503–6.

Greer, A.L., Cheng, Y.Q., Ma, E., 2013. Shear bands in metallic glasses. Mater. Sci. Eng. R (ISSN: 0927-796) 74 (4), 71–132. http://dx.doi.org/10.1016/j.mser.2013.04.001, URL https://www.sciencedirect.com/science/article/pii/S0927796X13000259.

Guillermín, N., Besson, J., Koster, A., Lacour, L., Maziere, M., Chalons, H., Forest, S., 2023. Experimental and numerical analysis of the portevin-le Chatelier effect in a nickel-base superalloy for turbine disks application. Int. J. Solids Struct. (ISSN: 0020-7683) 264, 112076. http://dx.doi.org/10.1016/j.ijsolstr.2022.112076, URL https://www.sciencedirect.com/science/article/pii/S0020768322005297.

Jiang, Z., Zhang, Q., Jiang, H., Chen, Z., Wu, X., 2005. Spatial characteristics of the Portevin-Le Chatelier deformation bands in Al-4at. Mater. Sci. Eng. A (ISSN: 0921-5093) 403 (1), 154–164. http://dx.doi.org/10.1016/j.msea.2005.05.059, URL https://www.sciencedirect.com/science/article/pii/S0921509305004934.

Kerfriden, P., 2024. Quantized plasticity 10cm, pierre kerfriden. https://www.youtube.com/watch?v=nW0KEG-W4jk.

Kiener, D., Grosinger, W., Dehm, G., Pippan, R., 2008. A further step towards an understanding of size-dependent crystal plasticity: In situ tension experiments of miniaturized single-crystal copper samples. Acta Mater. (ISSN: 1359-6454) 56 (3), 580–592. http://dx.doi.org/10.1016/j.actamat.2007.10.015, URL https://www.sciencedirect.com/science/article/pii/S1359645407006969.

Lamari, M., Allain, S.Y.P., Geandier, G., Ponçot, M., Perlade, A., Zhu, K., 2024. Behavior of TRIP-aided medium mn steels investigated by in situ synchrotron X-ray diffraction experiments and microstructure-based micromechanical modelling. Int. J. Plast. (ISSN: 0749-6419) 173, 103866. http://dx.doi.org/10.1016/j.jiplas.2023.103866, URL https://www.sciencedirect.com/science/article/pii/S0749641923003509.

Lebyodkin, M., Dunin-Barkowskii, L., Bréchet, Y., Estrin, Y., Kubin, L., 2000. Spatio-temporal dynamics of the Portevin–Le Chatelier effect: experiment and modelling. Acta Mater. (ISSN: 1359-6454) 48 (10), 2529–2541. http://dx.doi.org/10.1016/S1359-6454(00)000667-7, URL https://www.sciencedirect.com/science/article/pii/S1359645400000677.

Marano, A., Gelebart, L., Forest, S., 2019. Intragranular localization induced by softening crystal plasticity: Analysis of slip and kink bands localization modes from high resolution FFT-simulations results. Acta Mater. (ISSN: 1359-6454) 175, 262–275. http://dx.doi.org/10.1016/j.actamat.2019.06.010, URL https://www.sciencedirect.com/science/article/pii/S1359645419303696.

Mazière, M., Forest, S., 2015. Strain gradient plasticity modeling and finite element simulation of Lüders band formation and propagation. Contin. Mech. Thermodyn. (ISSN: 1432-0959) 27 (1), 83–104. http://dx.doi.org/10.1007/s00161-013-0331-8.

Oh, S.H., Legros, M., Kiener, D., Dehm, G., 2009. In situ observation of dislocation nucleation and escape in a submicrometre aluminium single crystal. Nature Mater. (ISSN: 1476-4660) 8 (2), 95–100. http://dx.doi.org/10.1038/nmat2370.

Okuyucu, C., Ulucan, T.H., Abboud, M., Motallebzadeh, A., Özerinç, S., Kalay, İ., Kalay, Y.E., 2023. Nanomechanical properties of Al-Tb marginal metallic glass. Mater. Sci. Eng. A (ISSN: 0921-5093) 888, 145809. http://dx.doi.org/10.1016/j.msea.2023.145809, URL https://www.sciencedirect.com/science/article/pii/S0921509323012339.

Patinet, S., Bonamy, D., Proville, L., 2011. Atomic-scale avalanche along a dislocation in a random alloy. Phys. Rev. B 84, 174101. http://dx.doi.org/10.1103/PhysRevB.84.174101.

Perchikov, N., Truskinovsky, L., 2024. Quantized plastic deformation. J. Mech. Phys. Solids (ISSN: 0022-5096) 190, 105704. http://dx.doi.org/10.1016/j.jmps.2024.105704, URL https://www.sciencedirect.com/science/article/pii/S0022509624001704.

Perilyev, V., Moshkovich, A., Lapsker, I., Laikhtman, A., Rapoport, L., 2013. The effect of vanadium content and temperature on stick-slip phenomena under friction of CrV(x)N coatings. Wear (ISSN: 0043-1648) 307 (1), 44–51. http://dx.doi.org/10.1016/j.wear.2013.08.012, URL https://www.sciencedirect.com/science/article/pii/S0043164813004651.

Richeton, T., Weiss, J., Louchet, F., 2005. Breakdown of avalanche critical behaviour in polycrystalline plasticity. Nature Mater. (ISSN: 1476-4660) 4 (6), 465–469. http://dx.doi.org/10.1038/nmat1393.

Ruestes, C.J., Segurado, J., 2024. A stochastic discrete slip approach to microplasticity: Application to submicron W pillars. Int. J. Plast. (ISSN: 0749-6419) 176, 103965. http://dx.doi.org/10.1016/j.ijplas.2024.103965, URL https://www.sciencedirect.com/science/article/pii/S0749641924000925.

Rýš, M., Kursa, M., Petryk, H., 2024. Spontaneous emergence of deformation bands in single-crystal plasticity simulations at small strain. Comput. Mech. (ISSN: 1432-0924) http://dx.doi.org/10.1007/s00466-024-02519-8.

Salman, O.U., Truskinovsky, L., 2011. Minimal integer automaton behind crystal plasticity. Phys. Rev. Lett. 106, 175503. http://dx.doi.org/10.1103/PhysRevLett.106.175503.

Salman, O.U., Truskinovsky, L., 2012. On the critical nature of plastic flow: One and two dimensional models. Internat. J. Engrg. Sci. (ISSN: 0020-7225) 59, 219–254. http://dx.doi.org/10.1016/j.ijengsci.2012.03.012, The Special Issue in honor of VICTOR L. BERDICHEVSKY, URL https://www.sciencedirect.com/science/article/pii/S0020722512000559.

Scroggs, M.W., Dokken, J.S., Richardson, C.N., Wells, G.N., 2022. Construction of arbitrary order finite element degree-of-freedom maps on polygonal and polyhedral cell meshes. ACM Trans. Math. Software 48, 18–1–18–23.

Simo, J., 1998. Numerical analysis and simulation of plasticity. In: Numerical Methods for Solids (Part 3) Numerical Methods for Fluids (Part 1). In: Handbook of Numerical Analysis, vol. 6, Elsevier, (ISSN: 1570-8659) pp. 183–499. http://dx.doi.org/10.1016/S1570-8659(98)80009-4, URL https://www.sciencedirect.com/science/article/pii/S1570865998800094.

Sornette, D., Ouillon, G., 2012. Dragon-kings: Mechanisms, statistical methods and empirical evidence. Eur. Phys. J. Spec. Top. (ISSN: 1951-6401) 205 (1), 1–26. http://dx.doi.org/10.1140/epjst/e2012-01559-5.

Uchic, M.D., Shade, P.A., Dimiduk, D.M., 2009. Micro-compression testing of fcc metals: A selected overview of experiments and simulations. JOM (ISSN: 1543-1851) 61 (3), 36–41. http://dx.doi.org/10.1007/s11837-009-0038-2.

Vermeij, T., Wijnen, J., Peerlings, R.H.J., Geers, M.G.D., Hoefnagels, J.P.M., 2024. A quasi-2D integrated experimental-numerical approach to high-fidelity mechanical analysis of metallic microstructures. Acta Mater. (ISSN: 1359-6454) 264, 119551. http://dx.doi.org/10.1016/j.actamat.2023.119551, URL https://www.sciencedirect.com/science/article/pii/S1359645423008807.

Wang, H., Berdin, C., Maziere, M., Forest, S., Prioul, C., Parrot, A., Le-Delliou, P., 2011. Portevin–Le Chatelier (PLC) instabilities and slant fracture in C–Mn steel round tensile specimens. Scr. Mater. (ISSN: 1359-6462) 64 (5), 430–433. http://dx.doi.org/10.1016/j scriptamat.2010.11.005, URL https://www.sciencedirect.com/science/article/pii/S13596464210007530.

Weiss, J., Rhouma, W.B., Richeton, T., Dechanel, S., Louchet, F., Truskinovsky, L., 2015. From mild to wild fluctuations in crystal plasticity. Phys. Rev. Lett. 114, 105504–1–105504–6.

Weiss, J., Richeton, T., Louchet, F., Chmelik, F., Dobron, P., Entemeyer, D., Lebyodkin, M., Lebedkina, T., Fressengeas, C., McDonald, R.J., 2007. Evidence for universal intermittent crystal plasticity from acoustic emission and high-resolution extensometry experiments. Phys. Rev. B 76, 224110–1–224110–8.

Wijnen, J., Peerlings, R.H.J., Hoefnagels, J.P.M., Geers, M.G.D., 2021. A discrete slip plane model for simulating heterogeneous plastic deformation in single crystals. Int. J. Solids Struct. 228, 111094–1–111094–12.

Xu, X., Cheng, L., Zhao, X., Wang, J., Chen, X., 2022a. Micro/nano periodic surface structures and performance of stainless steel machined using femtosecond lasers. Micromachines (ISSN: 2072-666X) 13 (6), http://dx.doi.org/10.3390/mi13060976, URL https://www.mdpi.com/2072-666X/13/6/976.

Xu, M., Zhao, Z., Wang, P., Duan, S., Lei, H., Fang, D., 2022b. Mechanical performance of bio-inspired hierarchical honeycomb metamaterials. Int. J. Solids Struct. (ISSN: 0020-7683) 254–255, 111866. http://dx.doi.org/10.1016/j.ijsolstr.2022.111866, URL https://www.sciencedirect.com/science/article/pii/S0020768322003341.

Yan, Z., Yan, J., Hu, Y., Dang, S., 2010. Crystallization in Zr60Al15Ni25 bulk metallic glass subjected to rolling at room temperature. Sci. China Ser. E: Technol. Sci. (ISSN: 1869-1900) 53 (1), 278–283. http://dx.doi.org/10.1007/s11431-009-0315-x.

Yilmaz, A., 2011. The Portevin–Le Chatelier effect: a review of experimental findings. Sci. Technol. Adv. Mater. 12, 063001–1–063001–16.

Yu, Q., Chatterjee, S., Roche, K.J., Po, G., Marian, J., 2021a. Coupling crystal plasticity and stochastic cluster dynamics models of irradiation damage in tungsten. Modelling Simul. Mater. Sci. Eng. 29 (5), 055021. http://dx.doi.org/10.1088/1361-651X/ac01ba.

Yu, Q., Legros, M., Minor, A.M., 2015. In situ TEM nanomechanics. MRS Bull. 40 (1), 62–70. http://dx.doi.org/10.1557/mrs.2014.306.

Yu, Q., Martínez, E., Segurado, J., Marian, J., 2021b. A stochastic solver based on the residence time algorithm for crystal plasticity models. Comput. Mech. 68, 1369–1384.

Zepeda-Ruiz, L.A., Stukowski, A., Oppelstrup, T., Bertin, N., Barton, N.R., Freitas, R., Bulatov, V.V., 2021. Atomistic insights into metal hardening. Nature Mater. 20, 315–320.
Zhang, P., Salman, O.U., Weiss, J., Truskinovsky, L., 2020. Variety of scaling behaviors in nanocrystalline plasticity. Phys. Rev. E 102, 023006. http://dx.doi.org/10.1103/PhysRevE.102.023006.

Zhang, P., Salman, O.U., Zhang, J.Y., Liu, G., Weiss, J., Truskinovsky, L., Sun, J., 2017. Taming intermittent plasticity at small scales. Acta Mater. 128, 351–364.

Zhu, Z., To, S., Ehmann, K.F., Xiao, G., Zhu, W., 2016. A novel diamond micro-/nano-machining process for the generation of hierarchical micro-/nano-structures. J. Micromech. Microeng. 26 (3), 035009. http://dx.doi.org/10.1088/0960-1317/26/3/035009.