# A plastic correction algorithm for full-field elasto-plastic finite element simulations: critical assessment of predictive capabilities and improvement by machine learning

Abhishek Palchoudhary, Simone Peter, Vincent Maurel, Cristian Ovalle, Pierre Kerfriden

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

## ► To cite this version:

Abhishek Palchoudhary, Simone Peter, Vincent Maurel, Cristian Ovalle, Pierre Kerfriden. A plastic correction algorithm for full-field elasto-plastic finite element simulations: critical assessment of predictive capabilities and improvement by machine learning. Computational Mechanics, 2024, 10.1007/s00466-024-02561-6. hal-04943408

HAL Id: hal-04943408

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

Submitted on 12 Feb 2025

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published 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_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# A plastic correction algorithm for full-field elasto-plastic finite element simulations: critical assessment of predictive capabilities and improvement by machine learning

Abhishek Palchoudhary<fnref n="1" />,<fnref n="2" /> · Simone Peter<fnref n="1" />,<fnref n="2" /> · Vincent Maurel<fnref n="1" /> · Cristian Ovalle<fnref n="1" /> · Pierre Kerfriden<fnref n="1" />

Received: 1 February 2024 / Accepted: 17 September 2024

© The Author(s) 2024

## Abstract

This paper introduces a new local plastic correction algorithm that is aimed at accelerating elasto-plastic finite element (FE) simulations for structural problems exhibiting localised plasticity (around e.g. notches, geometrical defects). The proposed method belongs to the category of generalised multi-axial Neuber-type methods, which process the results of an elastic prediction point-wise in order to calculate an approximation of the full elasto-plastic solution. The proposed algorithm relies on a rule of local proportionality, which, in the context of J2 plasticity, allows us to express the plastic correction problem in terms of the amplitude of the full mechanical tensors only. This lightweight correction problem can be solved for numerically using a fully implicit time integrator that shares similarities with the radial return algorithm. The numerical capabilities of the proposed algorithm are demonstrated for a notched structure and a specimen containing a distribution of spherical pores, subjected to monotonic and cyclic loading. As a second point of innovation, we show that the proposed local plastic correction algorithm can be further accelerated by employing a simple meta-modelling strategy, with virtually no added errors. At last, we develop and investigate the merits of a deep-learning-based corrective layer designed to reduce the approximation error of the plastic corrector. A convolutional architecture is used to analyse the neighbourhoods of material points and outputs a scalar correction to the point-wise Neuber-type predictions. This optional brick of the proposed plastic correction methodology relies on the availability of a set of full elasto-plastic finite element solutions to be used as a training data-set.

**Keywords** Model reduction · Plastic correction algorithm · Neuber's rule · Nonlinear finite element analysis · Convolutional neural networks

## List of Symbols

### Material behavior constants

$\lambda, \mu$ Lamé coefficients

$\sigma_y$ Yield stress

$C, D$ Kinematic hardening parameters

$Q, b$ Isotropic hardening parameters

Simone Peter, Vincent Maurel, Cristian Ovalle and Pierre Kerfriden
these authors contributed equally to this work

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

Abhishek Palchoudhary
abhishek.palchoudhary@minesparis.psl.eu

Pierre Kerfriden
pierre.kerfriden@minesparis.psl.eu

## Quantities from the elasto-static computation

$\bar{\sigma}_{\text{VM}}^{\#}$ Von Mises stress coming from the elasto-static computation at a fixed loading

$\zeta^{\#}, \xi^{\#}$ Stress and strain tensors coming from the elasto-static computation

$\zeta_d^{\#}, \xi_d^{\#}$ Deviatoric stress and strain tensors coming from the elasto-static computation

$\bar{\zeta}_d^{\#}, \bar{\xi}_d^{\#}$ Deviatoric stress and strain tensors coming from the elasto-static computation at a fixed loading

$\zeta_{d,o}^{\#}, \xi_{d,o}^{\#}$ Deviatoric stress and strain tensors coming from the elasto-static computation at the last peak of cyclic loading

## Quantities from a reference elasto-plastic computation

$\zeta, \xi$ Reference stress and strain tensors

$\zeta_d, \xi_d$ Reference deviatoric stress and strain tensors

<footnote n="1">¹ Centre for Material Sciences (CMAT), CNRS UMR 7633, BP 87, Mines Paris, PSL University, 91003 Evry, France</footnote>

<footnote n="2">² Technical University of Munich, Arcisstr. 21, 80333 Munich, Germany</footnote>
$\mathcal{E}^p$ Reference plastic strain tensor

$\mathcal{E}_d^p$ Reference deviatoric plastic strain tensor

$f_y$ Reference yield surface function

$\mathcal{J}$ Reference von Mises stress

$\mathcal{X}$ Reference nonlinear kinematic hardening tensor

$R$ Reference nonlinear isotropic hardening function

$p$ Reference cumulative plastic strain

## Quantities from an elasto-plastic computation by the plastic corrector

$\hat{\mathcal{E}}_d, \hat{\mathcal{E}}_d$ Deviatoric stress and strain tensors computed by the plastic corrector

$\hat{\mathcal{E}}_{d,o}, \hat{\mathcal{E}}_{d,o}$ Deviatoric stress and strain tensors at the last peak of cyclic loading computed by the plastic corrector

$\hat{\mathcal{E}}^p$ Plastic strain tensor computed by the plastic corrector

$\hat{\mathcal{E}}_d^p$ Deviatoric plastic strain tensor computed by the plastic corrector

$\hat{f}_y$ Yield surface function computed by the plastic corrector

$\hat{\mathcal{J}}$ Von Mises stress computed by the plastic corrector

$\hat{\mathcal{X}}$ Nonlinear kinematic hardening tensor computed by the plastic corrector

$\hat{p}$ Cumulative plastic strain computed by the plastic corrector

## Scalar variables for proportional evolution rule

$s$ Ratio of the approximated deviatoric stress tensor to the deviatoric stress tensor from the elasto-static computation

$e$ Ratio of the approximated deviatoric strain tensor to the deviatoric strain tensor from the elasto-static computation

$e^p$ Ratio of the approximated deviatoric plastic strain tensor to the deviatoric strain tensor from the elasto-static computation

$x$ Ratio of the approximated kinematic hardening tensor to the deviatoric strain tensor from the elasto-static computation

$s_o$ Ratio of the approximated deviatoric stress tensor to the deviatoric stress tensor from the elasto-static computation at the last peak of cyclic loading

$e_o$ Ratio of the approximated deviatoric strain tensor to the deviatoric strain tensor from the elasto-static computation at the last peak of cyclic loading

$e_o^p$ Ratio of the approximated deviatoric plastic strain tensor to the deviatoric strain tensor from the elasto-static computation at the last peak of cyclic loading

$f$ Loading function

$f_o$ Loading function at the last peak of cyclic loading

## Other miscellaneous notation

$\hat{\mathcal{Q}}_d$ Projected deviatoric stress tensor using the local proportionality rule

$\underline{u}$ Displacement vector

$\underline{\mathcal{I}}$ Identity tensor

$\Delta p$ Cumulative plastic strain range in a cycle

$\phi$ Intrinsic dissipation in a cycle

# 1 Introduction

The computer simulation of industrial components is often based on plasticity analysis around critical areas with high stress concentrations such as notches or defects like pores [1–8]. These analyses usually require finite element simulations. However, plasticity simulations are computationally expensive due to the number of elements required to accurately represent the geometry of stress concentrators, and due to the number of time increments needed to integrate nonlinear material laws over time. In this context, we propose a new methodology, which belongs to the class of plastic correction approaches, to rapidly obtain an approximation of the full-field elasto-plastic response of structures subjected to proportional loading sequences, from a single 3D elasto-static finite element solution.

Several types of plastic correction methodologies have been developed in the past to post-process the elasto-plastic response of structures from elastic finite element solutions. One family of methods uses homogenization theory, which consists in viewing the plastic zone as an inclusion in an elastic matrix and deriving local constraints to simulate the evolution of plastic quantities at the notch tip under load, starting from an elastic finite element solution [9–11]. Another family of methods relies on the use of Neuber-type rules. These rules are heuristics that relate the stresses and strains in an elasto-plastic body to those in a geometrically similar elastic body undergoing similar loading conditions. While originally developed for a uni-axial stress state at a notch subjected to monotonic loading [12] and then cyclic loading [13, 14], the Neuber rule (and other Neuber-type methods like the Equivalent Strain Energy Density approach [15]) have been generalized to multi-axial loading states at notch tips. One sub-family of Neuber-type methods for multi-axiality is based upon independently employing Neuber-type rules for every scalar component of the stress and strain tensors [16, 17]. The other sub-family of Neuber-type approaches reduces the complexity of the
previous approach by employing variations of proportional evolution rules for the stress tensor, strain tensor, or for a combination of stress and strain [18–21]. As far as we are aware, these existing pieces of work concentrate specifically on the development of Neuber-type plastic correction methods for the prediction of multiaxial stress and strain states at notch tips. One specific piece of work by Desmorat et al. [19] suggests using locally proportional evolutions together with Neuber-type rules to predict the evolution of elasto-plastic fields. Yet we could not find any detailed analysis of the accuracy of the suggested methodology away from notches and free boundaries, nor could we find the derivations of general-purpose algorithms to integrate the resulting plastic correction equations numerically, i.e. under arbitrary (proportional) load histories.

The first point of innovation proposed in this paper is the development of a multiaxial Neuber-type plastic correction method that may be used to produce approximations of elasto-plastic quantities at every quadrature point of a finite element mesh. As such, the developed methodology will not rely on particular strain state assumptions regarding the recovery of elasto-plastic solutions at traction-free boundaries. Our proposal is to use a local proportionality rule for the deviatoric stress and strain tensors, which will therefore be linearly related to their counterparts as calculated using an elasto-static analysis. Remarkably, in the context of J2 plasticity (we use a Chaboche model with kinematic and isotropic hardening), the local proportionality rule results in scalar constitutive equations for the deviatoric strain and stress amplitudes, without further assumption on plastic flow, which was already identified in [19]. Complemented by a scalar Neuber-type rule, the constitutive equations may be solved locally, the load stemming directly from that heuristic rule. We use the standard change of peaks method [14] to account for load cycles with non-zero mean stresses. We will analyse the accuracy of full 3D elasto-plastic solutions computed using the Neuber-type approximation, i.e. at and away from notches and free boundaries. This full-field aspect is particularly relevant to fields whereby fracture criteria are based on full-field elasto-plastic solutions, for example non-local fatigue models that use elasto-plastic fields around critical points [22], or weakest-link models that compute probabilities of failure via a weighted average of mechanical stresses over the entire computational domain [23, 24].

The second point of novelty is a machine-learning-based acceleration of the proposed plastic corrector, which is particularly useful for finite element models exhibiting a large number of degrees of freedom and long time analyses. This is because the local time integration of the elasto-plastic equations, using the plastic correction methodology mentioned previously, remains computationally expensive. We show that under the previously stated rules (J2 plasticity, proportionality of deviatoric strain and stress tensors, scalar Neuber

rule), and for a given load history, any output of the plastic correction algorithm exhibits a (nontrivial) scalar dependency to the von Mises stress of the elastic finite element simulation. Therefore, we suggest a meta-modelling strategy in the form of a Gaussian process regression [25], that will be trained to produce the elasto-plastic quantity of interest given a von Mises stress stemming from the elasto-static simulation as input. A dataset is generated by populating the real positive axis and computing the corresponding outputs using the plastic corrector. Few (30 to 150) such datapoints are necessary for a one-dimensional regression. We show that training the Gaussian process regression on such a small dataset is sufficient to act as a virtually cost-free surrogate for all remaining local plastic corrections, without sacrificing accuracy. As a consequence, the cost of acquiring approximated elasto-plastic quantities for specimens with single or multiple stress concentrators reduces to the cost of the elastic finite element simulation, as the number of degrees of freedom of the finite element model increases.

As a third and more exploratory part of the development, we propose a methodology of plastic correction using neural networks, which we coin Neural plastic corrector (NPC). The Neural Network developed here will be used to correct the output of the proposed Neuber-type methodology in order to better reproduce the mechanical fields delivered by a full elasto-plastic finite element analysis. While the Neuber-type method described thus far never requires elasto-plastic finite element simulations to be performed, the NPC is based on training with examples (supervised learning), and is therefore based on the availability of reference elasto-plastic solutions that will be used as a dataset. The Neuber corrector being local by nature, we aim to correct it by using information about the local topology of material point neighbourhoods, which will be analysed in an end-to-end fashion by a Convolutional Neural Network (CNN). More precisely, we will first voxelise the neighbourhood of any point of interest. We will then project the result of the previously described Neuber-type methodology to be improved onto the voxelised neighbourhood. Finally, we will provide this gridded data as input to a Neural network and train it to predict a correction to the elasto-plastic quantity of interest obtained approximately using the Neuber-type plastic corrector. Whilst the previously introduced Neuber-type method may be applied to general structures exhibiting stress concentrations, the Neural plastic corrector will be dedicated to the analysis of specimens with randomly placed pores, such as those that may be of interest when simulating the failure of porous alloys obtained by casting, welding or additive manufacturing [26–28]. In this context, the necessity to analyse specimens, or batches of specimens, containing large quantities of defects may justify the deployment of an AI that learns from a small to medium quantity of reference elasto-plastic simulations in order to predict the elasto-plastic output of subsequent simulations.
We do not anticipate a strong generalisation ability of the NPC, consistently with observations and analyses made in our previous work [29, 30], albeit in a different mechanical context. However, we will show that the output of the Neuber-type methodology may be improved upon by the proposed AI, using a reasonable amount of data for training, and that using the output of our Neuber-type methodology as full-field input to the CNN is indeed beneficial as compared to using inputs of the NPC with a lesser mechanical content.

This paper is divided into four sections: the first section introduces the proposed Neuber-type plastic corrector scheme. The second section is devoted to results and error analysis of the plastic correction algorithm for two different problems: A structure with a notch undergoing monotonic loading, and a specimen made of a porous alloy with a spherical pore population undergoing monotonic and cyclic loading. The third section introduces the meta-modelling strategy proposed to accelerate the plastic correction algorithm, and the fourth section presents the Neural Plastic Correction strategy. The three elements introduced in the paper are separately validated in the section where they are respectively presented. The paper is concluded by a discussion and perspectives of future studies.

# 2 Full field Neuber-type plastic corrector

## 2.1 Linear elasticity problem

Neuber-type methods approximate plasticity by locally processing the stress and strain fields stemming from a single elasto-static finite element analysis. We set up a linear elastic problem in an isotropic material, with stiffness tensor C and Lamé coefficients λ, μ, over a time interval [0, T], whereby proportional loading conditions are assumed at the global level, i.e. that of the structure. The time-dependent displacement, stress and strain tensors for the elasto-static simulation are denoted by $\underline{u}^{\#}(t)$, $\underline{\sigma}^{\#}(t)$ and $\underline{\varepsilon}^{\#}(t)$, respectively. The identity tensor is denoted by $\mathbb{I}$. The boundary $\partial\Omega$ of computational domain $\Omega$ is additively split into a Dirichlet part $\partial\Omega_u$ and a Neumann part $\partial\Omega_t$. The equations of linearised elasticity are introduced as follows:

$$ \mathrm{div} \, \underline{\sigma}^{\#}(t) + f(t) \underline{\xi} = 0 \qquad (1) $$

$$ \underline{\sigma}^{\#}(t) = C : \underline{\varepsilon}^{\#}(t) = \lambda \mathrm{Tr}(\underline{\varepsilon}^{\#}(t))\mathbb{I} + 2\mu \underline{\varepsilon}^{\#}(t) \qquad (2) $$

$$ \underline{\varepsilon}^{\#}(t) = \frac{1}{2}(\nabla \underline{u}^{\#}(t) + \nabla \underline{u}^{\#}(t)^T) \qquad (3) $$

$$ \underline{u}^{\#}(t) = f(t) \bar{u}_{\alpha} \quad \text{over } \partial\Omega_{u} \qquad (4) $$

$$ \underline{\sigma}^{\#}(t) \cdot \underline{n} = f(t) \bar{\underline{t}}_{\alpha} \quad \text{over } \partial\Omega_{t} \qquad (5) $$

In the previous set of equations, $f: [0, T] \mapsto \mathbb{R}$ is an arbitrary function of time that may be set to simulate, for example, proportionally monotonic and cyclic loading. $\bar{\underline{\xi}}$ is a vector-valued field of volume forces, $\bar{u}_\alpha$ is a vector-valued field of applied displacements, and $\bar{\underline{t}}_\alpha$ is a vector-valued field of applied traction loads. The previous elastic time-dependent problem is proportional, in the sense that all prescribed loading conditions, be it Neumann, Dirichlet conditions or volume sources, are introduced as fixed vector-valued fields multiplied by a function of time $f$.

Therefore, $\underline{u}^{\#}(t)$, $\underline{\sigma}^{\#}(t)$ and $\underline{\varepsilon}^{\#}(t)$ may be obtained by solving the equation of elasto-statics for $f=1$, and processed in the following way to recover the history of elastic solutions over $[0, T]$:

$$ \underline{u}^{\#}(t) = f(t) \bar{u}^{\#} \qquad (6) $$

$$ \underline{\sigma}^{\#}(t) = f(t) \bar{\underline{\alpha}}^{\#} \qquad (7) $$

$$ \underline{\varepsilon}^{\#}(t) = f(t) \bar{\underline{\xi}}^{\#} \qquad (8) $$

where $\bar{\underline{u}}^{\#}$, $\bar{\underline{\alpha}}^{\#}$ and $\bar{\underline{\xi}}^{\#}$ are, respectively, the displacement field, the stress tensor and the strain tensor fields obtained by setting $f=1$. (note: a bar symbol ($\bar{\cdot}$) is used to denote quantities that are obtained by solving the equations of elasto-statics with $f=1$).

If the elastic deviatoric stress and strain tensors obtained for $f=1$ are denoted by $\bar{\underline{\alpha}}_d^{\#}$ and $\bar{\underline{\xi}}_d^{\#}$ respectively, then the elastic deviatoric tensors at any time $t$ (denoted by $\underline{\alpha}_d^{\#}(t)$ and $\underline{\xi}_d^{\#}(t)$) are obtained by the following scaling:

$$ \underline{\alpha}_d^{\#}(t) = f(t) \bar{\underline{\alpha}}_d^{\#} \qquad (9) $$

$$ \underline{\varepsilon}_d^{\#}(t) = f(t) \bar{\underline{\varepsilon}}_d^{\#} \qquad (10) $$

## 2.2 Von Mises plasticity constitutive model

A von Mises plasticity model with a non-linear kinematic hardening and a non-linear isotropic hardening is chosen for the description of the evolution of the yield surface [31]. The deviator of the stress tensor $\underline{\sigma}$ is henceforth denoted by $\underline{\sigma}_d$, the total strain by $\underline{\varepsilon}$ and the plastic strain tensor by $\underline{\varepsilon}^p$. Elasticity is given by

$$ \underline{\sigma} = C : (\underline{\varepsilon} - \underline{\varepsilon}^p) \qquad (11) $$
The von Mises stress is given as:

$$ \mathcal{J}(\sigma_d - \mathcal{X}) = \sqrt{\frac{3}{2}(\sigma_d - \mathcal{X}) : (\sigma_d - \mathcal{X})} \qquad (12) $$

where $\mathcal{X}$ is the non-linear kinematic hardening tensor. The expression for the plastic strain rate $\dot{\varepsilon}^p$ is given as:

$$ \dot{\varepsilon}^p = \dot{p} \left( \frac{3}{2} \frac{\sigma_d - \mathcal{X}}{\mathcal{J}(\sigma_d - \mathcal{X})} \right) \qquad (13) $$

The evolution of $\mathcal{X}$ is given by the following expression:

$$ \dot{\mathcal{X}} = \frac{2}{3} C \dot{\varepsilon}^p - D \mathcal{X} \dot{p} \qquad (14) $$

where $\dot{p}$ is the time evolution of the cumulative plastic strain $p$, and $C$ and $D$ are kinematic hardening material parameters. The isotropic hardening $R(p)$ is given by the following expression:

$$ R(p) = Q(1 - \exp(-bp)) \qquad (15) $$

where $Q$ and $b$ are isotropic hardening material parameters. The deviator of the plastic strain tensor is used to define the cumulative plastic strain rate:

$$ \dot{p} = \sqrt{\frac{2}{3} \dot{\varepsilon}_d^p : \dot{\varepsilon}_d^p} \qquad (16) $$

The evolution of the yield surface $f_y(\sigma_d; \mathcal{X}, p)$ is required to satisfy the following two constraints:

$$ f_y(\sigma_d; \mathcal{X}, p) = \mathcal{J}(\sigma_d - \mathcal{X}) - \sigma_y - R(p) \le 0 \qquad (17) $$

$$ f_y \dot{p} = 0 \qquad (18) $$

## 2.3 Modified Neuber rule

### Neuber rule for deviatoric stress and strain tensors

The Neuber-type rule proposed in this paper operates on the deviatoric parts of the stress and strain tensors [16]:

$$ \hat{\sigma}_d : \hat{\varepsilon}_d = \sigma_d^{\#} : \varepsilon_d^{\#} \qquad (19) $$

where $\hat{\sigma}_d$, $\hat{\varepsilon}_d$ stand for the approximated deviatoric stress and strain tensors, respectively (note: a hat symbol ($\hat{\cdot}$) is used to denote quantities that are approximated by the plastic corrector).

Tensors $\sigma_d^{\#}$, $\varepsilon_d^{\#}$, which are obtained by solving elasto-statics instead of elasto-plasticity, are time-dependent, but the ($t$) notation has been dropped for conciseness.

Cyclic loading is handled with the classical change of origin at every peak, as proposed by Chaudonneret [14]:

$$ (\hat{\sigma}_d - \hat{\sigma}_{d,o}) : (\hat{\varepsilon}_d - \hat{\varepsilon}_{d,o}) = (\sigma_d^{\#} - \sigma_{d,o}^{\#}) : (\varepsilon_d^{\#} - \varepsilon_{d,o}^{\#}) \qquad (20) $$

In the previous equation, quantities $\hat{\sigma}_{d,o}$, $\hat{\varepsilon}_{d,o}$ stand for the approximated deviatoric stress and strain tensors at the last peak of loading. The tensors $\sigma_{d,o}^{\#}$, $\sigma_{d,o}^{\#}$ stand for the deviatoric stress and strain tensors coming from the elastic finite element solution evaluated at the last peak.

### Proportional evolution rule for deviatoric stress and strain tensors

It is important to note that stresses and strains in a structure depend on the geometry, the material behaviour and the boundary conditions. When there is a change in the ratio of any two components of the stress tensor at a given material point, these local stresses become non-proportional, by definition. If there is plastic flow anywhere in the structure, local non-proportionality may arise even in the case where the external loading is proportional.

The Neuber rule developed so far and the constitutive relations operate on symmetric second order tensors. We assume local proportionality [19], that is, while the local evolution of the actual deviatoric stress and strain tensors at a point may be arbitrarily complex, we postulate that the evolution is well approximated by assuming that it remains in the direction of the deviatoric stress and strain tensors obtained from elasto-statics, i.e. there is no shift in the direction in which plasticity develops.

Therefore, the deviatoric stress tensor may be written as a scaling of the stress tensor stemming from elasto-statics, which reads as

$$ \hat{\sigma}_d(t) = s(t) \bar{\sigma}_d^{\#} \qquad (21) $$

and similarly for the deviatoric strain tensor,

$$ \hat{\varepsilon}_d(t) = e(t) \bar{\varepsilon}_d^{\#} \qquad (22) $$

where the scalar variables $s$ and $e$ are introduced as scaling factors to be determined.

### Neuber rule in terms of proportionality ratios

Taking the local proportionality rules into account, the Neuber rule reduces to a constraint on scalar variables $s$ and $e$. By substituting equations (9) and (10), and (21) and (22) into equation (19), the monotonic Neuber rule, in the case of monotonic loading functions, becomes:

$$ s e \bar{\sigma}_d^{\#} : \bar{\varepsilon}_d^{\#} = f^2 \bar{\sigma}_d^{\#} : \bar{\varepsilon}_d^{\#} \qquad (23) $$

which yields

$$ s e = f^2 \qquad (24) $$
<figure><img src="image_4.png" /><figcaption>**Fig. 1** An illustration of the plastic correction algorithm during the (a) first branch of loading (b) second branch of loading, with $s_o$ and $e_o$ updated to their respective values at the last peak</figcaption></figure>

where we remind the reader that $s$, $e$ and $f$ are time-dependent variables (the (t) notation has been dropped for conciseness).

For cyclic loading, we introduce $s_o$ and $e_o$, which are the values of the scalar variables $s$ and $e$ at the last peak. Substituting equations (9) and (10), and (21) and (22) into equation (20):

$$ (s\bar{\bar{Q}}_d^\# - s_o\bar{\bar{Q}}_d^\#) : (e\bar{\bar{\xi}}_d^\# - e_o\bar{\bar{\xi}}_d^\#) = (f\bar{\bar{Q}}_d^\# - f_o\bar{\bar{Q}}_d^\#) : (f\bar{\bar{\xi}}_d^\# - f_o\bar{\bar{\xi}}_d^\#) \quad (25) $$

$$ \bar{\bar{Q}}_d^\# : \bar{\bar{\xi}}_d^\# (s - s_o)(e - e_o) = \bar{\bar{Q}}_d^\# : \bar{\bar{\xi}}_d^\# (f - f_o)(f - f_o) \quad (26) $$

$$ (s - s_o)(e - e_o) = (f - f_o)^2 \quad (27) $$

The scalar constraint between (deviatoric) stress and strain amplitudes is illustrated in Fig. 1 in the context of cycling loading. During the first branch of loading (denoted by the superscript 1), the quantities $s_o$ and $e_o$ take on the initial values of $s_o^1$ and $e_o^1$, i.e. 0, and $s$ and $e$ evolve as $s^1$ and $e^1$. When the peak of loading is reached, and a load reversal is made for a second branch of loading (denoted by the superscript 2), the values of $s_o$ and $e_o$ are updated to the previous peak values, shown in the figure as $s_o^2$ and $e_o^2$. Next, $s$ and $e$ evolve as $s^2$ and $e^2$ according to the updated origin. This process is repeated, with the quantities $s_o$ and $e_o$ successively taking on the values of $s$ and $e$ at every peak of loading following all the load reversals defined in the function $f(t)$. Equation (27) will be complemented by the elasto-plastic constitutive equations to yield a constitutive update algorithm.

## 2.4 Elasto-plastic constitutive equations under proportional tensor evolutions

The aim of this section is to reduce the tensorial equations of the constitutive law (presented in section 2.2) to scalar equations, using the rules of proportionality stated in equations (21) and (22), leading to a set of reduced constitutive equations in a set of scalar variables including $s$ and $e$. The reader is reminded that a hat symbol ($\hat{\cdot}$) is used to denote quantities approximated by the plastic corrector.

### Stress-strain relation

The approximated stress tensor is given by:

$$ \hat{\hat{\alpha}} = 2\mu(\hat{\xi} - \hat{\xi}^p) + \lambda \mathrm{Tr}(\hat{\xi} - \hat{\xi}^p)\mathbb{I} \quad (28) $$

The approximated stress tensor $\hat{\hat{\alpha}}$ can be split into its approximated deviatoric $\hat{\hat{\alpha}}_d$ and hydrostatic $\hat{\hat{\alpha}}_h$ parts:

$$ \hat{\hat{\alpha}} = \hat{\hat{\alpha}}_d + \hat{\hat{\alpha}}_h \quad (29) $$

where the approximated deviatoric stress tensor $\hat{\hat{\alpha}}_d$, under isochoric plastic flow, becomes:

$$ \hat{\hat{\alpha}}_d = 2\mu(\hat{\xi}_d - \hat{\xi}_d^p) \quad (30) $$

The deviatoric strain is approximated by assuming proportional evolution to the deviatoric strain coming from the elasto-static simulation (refer to Sect. 2.3). By virtue of equations (13) and (30) the plastic deviatoric strain evolves proportionally to the deviatoric strain coming from the elasto-static simulation, $\bar{\bar{\xi}}_d^\#$. A scalar variable $e^p$ can be introduced to encode this:

$$ \hat{\hat{\xi}}_d^p(t) = e^p(t)\bar{\bar{\xi}}_d^\# \quad (31) $$
<table><caption>**Table 1** Plastic correction equations written in terms of tensor variables ($\hat{\alpha}_d, \hat{\varepsilon}_d, \hat{\varepsilon}^p, \hat{p}, \hat{\mathbf{X}}$) to the left, and in terms of scalar variables ($s, e, e^p, \hat{p}, x$) to the right. *f* is the global loading function, which is time dependent,</caption><thead><tr><th>Type</th><th>Tensorial variables</th><th>Scalar variables</th></tr></thead><tbody><tr><td>Neuber rule</td><td>$\hat{\alpha}_d : \hat{\varepsilon}_d = \alpha_d^# : \varepsilon_d^#$</td><td>$s = f^2$</td></tr><tr><td>Elasticity</td><td>$\hat{\alpha} = 2\mu(\hat{\varepsilon} - \hat{\varepsilon}^p) + \lambda \text{Tr}(\hat{\varepsilon} - \hat{\varepsilon}^p)\mathbb{I}$</td><td>$s = (e - e^p)$</td></tr><tr><td>Yield function</td><td>$f_y(\hat{\alpha}; \hat{\mathbf{X}}, \hat{p}) = \mathcal{J}(\hat{\alpha} - \hat{\mathbf{X}}) - \sigma_y - R(\hat{p})$</td><td>$\hat{f}_y(s; x, \hat{p}) = \mathcal{J}(s, x) - \sigma_y - R(\hat{p})$</td></tr><tr><td>Evolution of yield function</td><td>$f_y(\alpha_d; \mathbf{X}, \hat{p}) \dot{\hat{p}} = 0$ and $f_y(\alpha_d; \mathbf{X}, \hat{p}) \le 0$</td><td>$\hat{f}_y(s; x, \hat{p}) \dot{\hat{p}} = 0$ and $\hat{f}_y(s; x, \hat{p}) \le 0$</td></tr><tr><td>von Mises stress</td><td>$\mathcal{J}(\hat{\alpha}_d - \hat{\mathbf{X}}) = \sqrt{\frac{3}{2}(\hat{\alpha}_d - \hat{\mathbf{X}})} : (\hat{\alpha}_d - \hat{\mathbf{X}})$</td><td>$\hat{\mathcal{J}}(s, x) = \left|s - \frac{x}{2\mu}\right| \bar{\sigma}_{\text{VM}}^{\#}$</td></tr><tr><td>Isotropic hardening</td><td>$R(\hat{p}) = Q(1 - \exp(-b\hat{p}))$</td><td>$R(\hat{p}) = Q(1 - \exp(-b\hat{p}))$</td></tr><tr><td>Kinematic hardening</td><td>$\dot{\hat{\mathbf{X}}} = \frac{2}{3} C \dot{\hat{\varepsilon}}^p - D \hat{\mathbf{X}} \dot{\hat{p}}$</td><td>$\dot{x} = \frac{2}{3} C \dot{e}^p - D x \dot{\hat{p}}$</td></tr><tr><td>Flow rule</td><td>$\dot{\hat{\varepsilon}}^p = \dot{\hat{p}} \left( \frac{3}{2} \frac{\hat{\alpha}_d - \hat{\mathbf{X}}}{\mathcal{J}(\hat{\alpha}_d - \hat{\mathbf{X}})} \right)$</td><td>$\dot{e}^p = \dot{\hat{p}} \left( \frac{3}{2} \frac{2\mu s - x}{\mathcal{J}(s, x)} \right)$</td></tr><tr><td>Cumulative plastic strain</td><td>$\dot{\hat{p}} = \sqrt{\frac{2}{3} \dot{\hat{\varepsilon}}_d^p : \dot{\hat{\varepsilon}}_d^p}$</td><td>$\dot{\hat{p}} = \frac{1}{3\mu} |\dot{e}^p| \bar{\sigma}_{\text{VM}}^{\#}$</td></tr></tbody></table>

and $\bar{\sigma}_{\text{VM}}^{\#}$ denotes the von Mises stress stemming from the elastic finite element simulation performed with $f = 1$

<table><caption>**Table 2** Parameters of the reference elasto-plastic equations detailed in Eqs. (11)-(18) and the plastic corrector equations detailed in Table 1, used for test cases 1 and 2</caption><thead><tr><th>E<br>MPa</th><th>σ<sub>y</sub><br>MPa</th><th>b</th><th>Q<br>MPa</th><th>C<br>MPa</th><th>D</th></tr></thead><tbody><tr><td>200000</td><td>100</td><td>10</td><td>100</td><td>40000</td><td>400</td></tr></tbody></table>

For arbitrary (proportional) load functions $f(t)$, this is done by computing the value of $e^p(t)$ at a series of time steps, using the fully implicit time stepping scheme described in the appendix A.

We can finally reconstruct the approximated elasto-plastic tensor variables ($\hat{\alpha}(t), \hat{\varepsilon}(t), \hat{\varepsilon}^p(t), \hat{p}(t), \hat{\mathbf{X}}(t)$) from the scalar quantities computed using the plastic corrector, i.e. ($s(t), e(t), e^p(t), \hat{p}(t), x(t)$). For example, to compute a component of the stress tensor, we can first reconstruct $\hat{\alpha}(t) = s(t)\hat{\alpha}_d + f(t)\text{Tr}(\hat{\alpha}^{\#})\mathbb{I}$ and then report the time evolution of desired component of this approximated stress tensor.

## 3 Numerical investigations and results

The aim of this section is to evaluate the accuracy of the plastic corrector with respect to reference computations, which consist of full elasto-plastic FE simulations using the Z-Set suite [32]. Two test cases, a notched structure and a specimen with spherical pores, will be presented with detailed full-field and point-wise comparisons of variables like the cumulative plastic strain. A third test case, that of a specimen with pores of realistic morphology, meshed based on images from X-ray tomography, is presented in Appendix B. The meshes for the test cases were created with GMSH [33],

and elasto-static simulations needed for the plastic corrector methodology were performed with FENiCS [34].

The parameters used for the reference elasto-plastic equations (equations (11)-(18)) and the plastic corrector equations (detailed in Table 1) chosen for this investigation, for test cases 1 and 2, are summarized in Table 2.

### 3.1 Test cases

#### Notched structure

Test case 1 is a notched structure, illustrated in Fig. 2a. For the plastic corrector, an elasto-static FEA computation is computed with prescribed displacements $\bar{u}_a = [u_x, 0, 0]$ for $x \ge L_c$ and $\bar{u}_a = [-u_x, 0, 0]$ for $x \le -L_c$ such that $\bar{\sigma}_{\text{VM}}^{\#}$ at the notch area is at the yield stress of the material (this computation corresponds to $f = 1$). The function $f(t)$ is then chosen to monotonically increase such that $f(t)\bar{\sigma}_{\text{VM}}^{\#}$ reaches 155% of the yield stress of the material at the notch tip. This monotonic load function is shown in Fig. 2b. The full-field $\bar{\sigma}_{\text{VM}}^{\#}$ in the structure is used for the plastic corrector approximations. The reference elasto-plastic computation is performed in such a way as to achieve the same loading level as in the plastic corrector computation.

#### Specimen with spherical pores

Test case number 2 is a specimen with a spatial distribution of spherical pores in the gauge section shown in Fig. 3. The spherical pores are of a fixed size, and their centers were generated using random sampling. For the plastic corrector, an elasto-static FEA computation is computed with prescribed displacements $\bar{u}_a = [u_x, 0, 0]$ for $x \ge L_c$ and $\bar{u}_a = [-u_x, 0, 0]$ for $x \le -L_c$ such that $\bar{\sigma}_{\text{VM}}^{\#}$ in the gauge section, away from pores, is at the yield stress of the material (this computation corresponds to $f = 1$). The function $f(t)$
<figure><img src="image_5.png" /><figcaption>**Fig. 2** Test case 1: Geometry, boundary conditions and applied loading</figcaption></figure>

(a) Notched structure, with displacement $\bar{u}_a$ at $f = 1$ scaled by $f(t)$ on highlighted red regions, in opposite directions

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

(b) At the notch tip, the von Mises stress from the elastic computation $\bar{\sigma}_{\text{VM}}^{\#}$ at $f = 1$ is scaled by $f(t)$. $f(t)\bar{\sigma}_{\text{VM}}^{\#}$ goes up to 155% of the yield stress at the peak of loading. The red points indicate the specific time-steps at which the cumulative plastic strain results from the full-field plastic corrector will be compared against the reference cumulative plastic strain.

is evaluated twice, generating loading sequences of 2 cycles and 20 cycles. At the peak of these cyclic loads, $f(t)\bar{\sigma}_{\text{VM}}^{\#}$ is chosen to reach 80% of the yield stress of the material in the gauge section away from pores. These cyclic load functions are shown in Fig. 3(b-c). Due to the presence of pores in the specimen, the local loading exceeds the yield stress in several regions. The full-field $\bar{\sigma}_{\text{VM}}^{\#}$ in the specimen is used for the plastic corrector approximations. For the reference elasto-plastic computation, the same type of loading as in the plastic corrector is applied.

## 3.2 Accuracy of the plastic corrector predictions

The accuracy of some quantities of interest computed using the plastic corrector will now be examined for the previously presented test cases. For test case 1, the evolution of the cumulative plastic strain, at a point of the notch tip and for the full 3D geometry will be shown, for monotonic loading. For test case 2, the evolution of the cumulative plastic strain will be shown in the 3D geometry encompassing the pores, for both monotonic and cyclic loading. The accuracy of the

deviatoric stresses and strains will also be evaluated for a few points around the pores. Finally, the accuracy of the cumulative plastic strain range (denoted by $\Delta p = p_{\text{max}}^{\text{cycle}} - p_{\text{min}}^{\text{cycle}}$) over the 20th cycle will be evaluated. The 20th cycle is chosen here because in subsequent cycles, the stress strain responses do not change significantly from one cycle to the next, i.e. the cycling response is stabilised. Quantities of interest extracted from such stabilised cycles are notably suitable for the computation of fatigue criteria [35].

### Time-evolution of cumulative plastic strain

The field of cumulative plastic strain approximated by the plastic corrector, at three time steps of the loading sequence, is shown for test case 1 in Fig. 4 and compared with reference computations at the same levels of applied loading. The time steps correspond to the notch von Mises stress (of the elastic FEA computation) reaching 102%, 127% and 152% of the elastic yield limit. The plastic corrector predicts overall higher cumulative plastic strain in the notch region, as compared to the reference elasto-plastic solution. Scatter plots of the cumulative plastic strain in all integration points in the mesh are shown in Fig. 5. The approximated solution of
<figure><img src="image_7.png" /><figcaption>(a) Specimen with spherical pores in gauge section, with displacement $\bar{u}_a$ at $f = 1$ scaled by $f(t)$ on highlighted red regions, in opposite directions</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>(b) Load function 1: In the gauge section (c) Load function 2: In the gauge section away from pores, the von Mises stress away from pores, the von Mises stress from the elastic computation $\bar{\sigma}_{\text{VM}}^{\#}$ at $f = 1$ is scaled by $f(t)$. The plastic corrector 1 is scaled by $f(t)$. The plastic corrector and reference cumulative plastic strain and the reference cumulative plastic strain will be compared at the specific time-steps range over the 20th cycle (red), will be indicated in red compared</figcaption></figure>

the plastic corrector algorithm aligns closely with the reference solution. However, there exists some scatter, which is expected due to the cumulative errors linked to the rules of the plastic correction algorithm. The reference computation has points that undergo plastic accommodation, which is not captured in the Neuber-type computation, therefore leading to over-estimation by the plastic correction algorithm. This is a known phenomenon and has been reported before [10]. Secondly, the Neuber approach also under-estimates the cumulative plastic strain at some points as it is localised and does not take into account redistribution of stresses that occur in the reference computation [15].

The time-evolution of a maximally loaded point on the notch region is shown in Fig. 6. The plastic corrector over-estimates the cumulative plastic strain throughout the loading sequence because stress redistribution was not considered. The relative error is very high at first, due to very low values of plasticity. During loading at the notch between 127% – 152%$\sigma_y$, the relative error lowers to around 15–20%. These over-estimations match with the literature [15, 36].

For test case 2, the cumulative plastic strains obtained using the plastic corrector at a couple of time steps of the applied loading (corresponding to 40% and 80% of the yield stress of the material) are shown in Fig. 7 and compared with reference computations at the same levels of applied loading. Plasticity develops around the pores, and the plastic corrector

once again predicts overall higher cumulative plastic strain in a slightly larger region, as compared to the reference. Scatter plots of the cumulative plastic strain in all integration points in the mesh are shown in Fig. 8. The full-field predictions match well with the reference, especially for the time instants corresponding to lower levels of loading. This is because the plastic corrector’s accuracy improves when the plasticity is relatively confined. However, a slight over-estimation by the plastic correction algorithm is observed.

The time-evolution of a maximally loaded point near a pore is shown in Fig. 9a–b. Similar to case 1, the plastic corrector over-estimates the cumulative plastic strain throughout the loading sequence. Interestingly, the relative error is lower at the beginning of loading and rises afterwards, contrary to case 1. The relative error in the prediction at the highest stressed points is around 35 – 40% during the cyclic loading. The error oscillates with load reversals, with the oscillations reducing with successive cycles. The same analysis is done with a point seeing lesser loading, shown in Fig. 9c–d. Conversely, the plastic corrector under-estimates the cumulative plastic strain throughout the loading sequence. Under-estimation at points away from highly loaded regions is also a result of stress re-distribution, and has rarely been reported in the literature as Neuber-type methods have primarily been studied for highly loaded points or surfaces.
<figure><img src="image_9.png" /><figcaption>**Fig. 4** Test case 1 (notched geometry): Time-evolution of cumulative plastic strain</figcaption></figure>

## Time-evolution of stresses and strains

We limit the study of stresses and strains to the more complex loading case, i.e. the test case 2 with spherical pores. We evaluate the accuracy of the plastic corrector on components of the deviatoric stress and strain tensors. The same two points are selected as the previous section, i.e. a highly loaded one with high relative error in cumulative plastic strain (Fig. 10a) and another point with lesser loading and lower relative error (Fig. 11a).

Figure 10b–c shows the evolution of the axial and shear components of the stress and strain tensors for the point under high local loading.

While the stresses are relatively well approximated, the strains have a much higher error–this is a known result [10,

19, 36]. The axial component is overestimated by the plastic corrector. As aforementioned, the Neuber-type method used does not take into account redistribution of stresses. The shear component is not well predicted due to the proportional evolution rule in the plastic corrector–this will be explained in detail in a later section.

Figure 11b–c shows the evolution of the axial and shear components of the stress and strain tensors for the point under low local loading.

Both the axial and shear components are reasonably well-estimated by the plastic corrector. This is because this point, experiencing lesser loading than the first one, is not as affected by the proportional evolution rule–again, this is proved later on.
<figure><img src="image_10.png" /><figcaption>**Fig. 6** Test case 1 (notched geometry): Comparison of the time-evolution of the cumulative plastic strain at the notch tip obtained by the plastic corrector and the reference computation during monotonic loading</figcaption></figure>

<figure><img src="image_11.png" /></figure>
<figure><img src="image_12.png" /><figcaption>**Fig. 7** Full-field cumulative plastic strain results computed by the plastic correction algorithm compared to the reference, for two time steps corresponding to successively higher levels of loading</figcaption></figure>

Finally, the von Mises stress approximated by the plastic corrector is compared to the reference computation in Fig. 12. The results show an excellent match between the plastic corrector and the reference, for progressively higher loading.

**Cumulative plastic strain range in the 20th cycle**

The field of cumulative plastic strain range $\Delta p$ over the 20th cycle, obtained by the plastic correction algorithm and the reference elasto-plastic FE computation, is shown respectively in Fig. 13a and b for the test case number 2, specimen with spherical pores. A comparison of the values at all the quadrature points in the mesh is shown in Fig. 13c. A good

overall match is found in the values of the cumulative plastic strain range obtained by the plastic correction algorithm and the reference computation. Like the previously presented cases, there is over and under-estimation by the plastic corrector – the number of over-estimated points is significantly higher than the under-estimated ones. The order of relative error in the loaded zones is around 30 – 40%. However, we acknowledge that the errors may be much higher if the material presents ratcheting behaviour, which was not the case here [10].
<figure><img src="image_13.png" /><figcaption>**Fig. 9** Comparison of the time-evolution cumulative plastic strain obtained by the plastic corrector and a reference computation at two points in the mesh of test case 2, showing an over-estimation and under-estimation of the plastic corrector</figcaption></figure>

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

(a) Selected point in the mesh with high loading and high relative error in $p$

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

(b) Evolution of cumulative plastic strain at the shown point with high loading and high error

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

(c) Selected point in the mesh with low loading and low relative error in $p$

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

(d) Evolution of cumulative plastic strain at the shown point with low loading and low error
<figure><img src="image_18.png" /><figcaption>(a) Selected point in the mesh with high error</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>(b) Axial component</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>(c) Shear component</figcaption></figure>

## 3.3 Computational time needed for a full-field plastic corrector computation

A breakdown of the computational costs required for obtaining the elasto-plastic solution by the two methods, i.e. a full reference computation and a plastic corrector computation, is compared in Table 3. The mesh considered is the same porous mesh as test case 2, and has around 661k quadrature points (one per linear element). The loading considered here comprises 1000 time-steps.

## 3.4 Errors due to the rule of local proportionality

We will now analyze the error that arises solely from the local proportionality rule. Our objectives are to evaluate (i) the contribution of this proportionality rule to the error in the plastic corrector method and (ii) the characteristics of this error. To achieve this, we will construct a projected elasto-plastic solution based exclusively on the local proportionality rule, omitting the Neuber-type rule. Specifically, we will project

the reference stress history at a given point onto the direction of the stress derived from the elastic computation. This projected stress will be compared to the stress calculated using the plastic corrector, which incorporates both the Neuber-type rule and the local proportionality rule. This comparison will allow us to isolate the part of the error that is due to the local proportionality rule. Additionally, we will compare the projected stress computed using only the proportionality rule to the reference stress to assess the acceptability of the proportionality rule in terms of error. Deviations between the stresses computed with these two methods will indicate significant local non-proportionality, which our algorithm may not adequately capture. We will show the characteristics of this error, i.e. the conditions and reasons for which usage of the proportionality rule could result in high relative error.

The projected stress tensor resulting from the local proportionality rule is computed using the following projection of the reference stress $\sigma_d$:

$$ \hat{\bar{\sigma}}_d = \frac{(\bar{\sigma}_d : \sigma_d)}{(\bar{\sigma}_d : \bar{\sigma}_d)} \bar{\sigma}_d \qquad (46) $$
<figure><img src="image_21.png" /><figcaption>**Fig. 11** Point in the mesh of test case 2 with a low amount of plasticity and low error in $p$ (9% relative error between reference $p$ and plastic corrector $\hat{p}$ at $t=0.12$ s): evolution of the individual components of the approximated deviatoric stress and strain tensors by the plastic corrector (denoted by $\hat{\sigma}_d^{ij}$, $\hat{\varepsilon}_d^{ij}$) compared to the respective reference curves obtained via a complete elasto-plastic computation (denoted by $\sigma_d^{ij}$, $\varepsilon_d^{ij}$)</figcaption></figure>

Quantities denoted by $\hat{\cdot}$ represent projected quantities assuming local proportionality; this approximation is thus different from the plastic corrector approximations developed in the previous sections as it directly operates on the reference stress (which is, of course, not available in practice). The absolute error $\xi$ and the relative error $\xi_{rel}$ between the projected deviatoric stress tensor using the hypothesis of proportionality ($\hat{\mathcal{Q}}_d$), and the reference deviatoric stress tensor via a full non-linear FE computation ($\sigma_d$) can be calculated as follows:

$$ \xi = \sigma_d - \hat{\mathcal{Q}}_d \quad (47) $$

The relative error in stress is computed by taking the Frobenius norm of the tensors:

$$ \xi_{rel} = \frac{\|\xi\|_F}{\|\sigma_d\|_F} \quad (48) $$

### Contribution of the local proportionality rule to the error in the plastic corrector solution

We will now analyze the contribution of the proportional evolution rule in terms of error in the plastic corrector solution, at a few points of the computational domain of test case 2. The same two points previously analyzed in Fig. 9a,c are considered. Figure 14a–b shows a comparison of the evolution of the axial and shear components of the deviatoric stress tensors for the first point near a pore, with high plasticity and high error, computed using (i) purely the proportional evolution rule (denoted by $\hat{\sigma}_d^{ij}$) (ii) the plastic corrector (denoted by $\hat{\sigma}_d^{ij}$) (iii) reference (denoted by $\sigma_d^{ij}$). For the components computed using purely the local proportional evolution rule, we observe that the axial component is very well predicted. The shear component is badly predicted, but the maximum stress value is acceptable. This shows that in the plastic corrector solution, the proportional evolution rule does not contribute significantly to the overall solution, and it is the Neuber-type rule that causes most of the error. A similar comparison is done for the second point with low
<figure><img src="image_22.png" /><figcaption>**Fig. 12** Full-field von Mises stress results computed by the plastic correction algorithm compared to the reference, for two time steps corresponding to successively higher levels of loading</figcaption></figure>

plasticity and low error (Fig. 15a–b). For this point, both the axial and shear components computed with the local proportional evolution rule are reasonably well predicted, and also do not differ from the plastic corrector solution.

The relative error $\xi_{rel}$ due to the local proportionality rule is shown for the two points in Figs. 14c and 15c, alongside their cumulative plastic strain. The load reversal naturally causes $\xi_{rel}$ to increase sharply. $\xi_{rel}$ becomes close to 1, which indicates a 100% relative error. The explosion of relative errors at these load inversion points is because the stress values are very close to zero (see Figs. 14a,b and 15a,b).

The absolute error at these inversion points is very low. This shows that even if the relative error is significant, in reality, its effect is negligible.

Furthermore, the cumulative plastic strain does not increase during these load reversals. In other words, the error during load reversal does not contribute significantly to the error in plasticity computations. For the first point with a higher *p* error, cumulative plastic strain starts increasing again despite a high $\xi_{rel}$. For the second point, with lower *p* error, we observe that the $\xi_{rel}$ always stays low during increase in plasticity. Thus, some points of loading, which are closer
<figure><img src="image_23.png" /><figcaption>**Fig. 13** Full-field comparison of $\Delta p$ in the 20th cycle in the test case 2 mesh</figcaption></figure>

to a proportional loading sequence, are better predicted than others.

**Relative error due to the local proportionality rule evaluated in the full mesh**

Only test case number 2, i.e. the specimen with spherical pores, is considered. The relative error $\xi_{rel}$ due to the local proportionality rule is calculated for all the loading history. The percentage of points in the mesh that stay below 15%

relative error during the loading sequence is shown in Fig. 16. During the first branch of monotonic loading, around 70-75% of all the points in the mesh remain below 15% relative error, and during the load reversal, this number goes down to 20%. As previously discussed, this occurs because the near-zero stress values at these load inversion points lead to a sharp increase in relative errors. The percentage of points below 15% relative error goes back up to 70% after the load reversal, when plasticity starts developing again.
<table><caption>**Table 3** CPU times (in seconds) for full-field *p* computation for the test case 2: mesh of a specimen with spherical pores (661771 quadrature points, 1000 time-steps)</caption><thead><tr><th>Operation</th><th>Reference FEA (s)</th><th>Plastic corrector (s)</th></tr></thead><tbody><tr><td>Elastic FEA for σ&#773;<sup>#</sup><sub>VM</sub></td><td>-</td><td>236</td></tr><tr><td>Elasto-plastic computation</td><td>52850</td><td>156</td></tr><tr><td>Total CPU time</td><td>52850</td><td>392</td></tr></tbody></table>

<figure><img src="image_24.png" /><figcaption>**Fig. 14** For a point with high plasticity (shown in Fig. 9a): Evolution of the individual components of the projected deviatoric stress tensor (denoted by $\hat{\sigma}_{d}^{ij}$) compared to the respective reference curves obtained via a complete elasto-plastic computation without any reduction (denoted by $\sigma_{d}^{ij}$). The plastic corrector solution is recalled in red. The relative error fraction for the stresses is shown alongside the cumulative plastic strain</figcaption></figure>

## 4 Machine learning-based acceleration of plastic corrector computations

We wish to compute a scalar quantity of interest (QoI) from the elasto-plastic solution at arbitrary time $t \in [0, T]$. More precisely, we wish to extract a scalar value $Q \in \mathbb{R}$ from $(\hat{\mathcal{X}}(t), \hat{\mathcal{E}}(t), \hat{\mathcal{E}}^p(t), \hat{p}(t), \hat{\mathcal{X}}^*(t))$ at an arbitrary number of quadrature points of the finite element mesh.

During the process of acquiring the QoI, the computation of $(s(t), e(t), e_p(t), \hat{p}(t), x(t))$ requires implicit time integration, as described in section 2.4 and the appendix A, and needs to be done for every quadrature point of the finite element mesh, which may lead to significant computational expense, especially for long time analysis. Remarkably, variables $(s(t), e(t), e_p(t), \hat{p}(t), x(t))$ are dependent on the solution of the elastic finite element prediction through the von Mises stress only, which can be seen by inspection of the system of equations in Table 1, right column. Therefore,

any elasto-plastic quantity of interest, which depends solely on these proportionality ratios also depends on the sole von Mises stress $\bar{\sigma}_{\text{VM}}^{\#}$ stemming from the elastic finite element computation, for a given $f(t)$.

We propose to acquire the scalar QoI with a one-dimensional Gaussian process regression (GP) algorithm [25] trained on plastic corrector computations. One may use several independent GP regression models for several quantities of interest.

### Gaussian process regression algorithm

The GP needs to capture the dependence of the scalar QoI on $\bar{\sigma}_{\text{VM}}^{\#}$. To this end, we need to define a training interval for $\bar{\sigma}_{\text{VM}}^{\#}$. A uniformly distributed set of $n_s$ von Mises stress values, between 0 and $s^+$ times the $\sigma_y$, is created in logarithmic space. As stress values can span several orders of magnitude, usage of logarithmic space ensures that the GP has good interpolation ability. These values are input to the plastic corrector to calculate the QoI. $\bar{\sigma}_{\text{VM}}^{\#}$ and QoI constitute the training data
<figure><img src="image_25.png" /><figcaption>**Fig. 15** For a point with low plasticity (shown in Fig. 9c): Evolution of the individual components of the projected deviatoric stress tensor calculated by assuming a hypothesis of proportionality (denoted by $\hat{\sigma}_d^{ij}$) compared to the respective reference curves obtained via a complete elasto-plastic computation without any reduction (denoted by $\sigma_d^{ij}$). The plastic corrector solution is recalled in red. The relative error fraction for the stresses is shown alongside the cumulative plastic strain</figcaption></figure>

<figure><img src="image_26.png" /><figcaption>**Fig. 16** Time evolution of the percentage of elements in the FE computation (test case 2: spherical pores) below 15% relative error, alongside the applied loading history in the gauge section away from pores</figcaption></figure>

for the GP. For numerical stability, we applied a logarithmic transformation to both the input and target data. Any zero values in the target data were replaced with a small positive value to avoid undefined logarithmic values. The model was then fitted to the log-transformed data.

Once trained, the GP can be applied to all integration points within a given mesh of a structure without requiring any time integration.

**Numerical example on a QoI** As an example we train a GP to predict $e^p$ at the 1000th time-step of a given $f(t)$ (shown in Fig. 17a) for the mesh of a specimen with spherical pores (test case 2). The training data set for learning the scalar relation between $\bar{\sigma}_{VM}^\#$ and $e^p$ is created using $n_s = 150$ points and $s^+ = 12$. This is shown in Fig. 17b. The choice of 150 points ensures sufficient resolution for accurately capturing the relationship between the two quantities, which is either zero or a simple monotonically increasing function. The upper limit, set at 12 times the yield stress, is chosen to be slightly higher than the maximum stress concentration factor induced by the pores in the mesh, which is approximately 8.

A comparison between the results obtained via the plastic correction algorithm and the predicted results via the Gaussian process regression is shown in Fig. 18 for all integration points. The results show virtually no difference, indicating that the 1D meta-model can be used to further accelerate computation of scalar QoIs of the plastic corrector with no added error.

The Gaussian process regression only needs a few seconds in total to use for the prediction of a scalar QoI, including the cost of computing the training set via the plastic corrector, training time (computation and factorisation of data covariance matrix), and prediction time (shown in Table 4).
<figure><img src="image_27.png" /><figcaption>**Fig. 17** a The load function chosen, with 1000 time-steps, with $e^p$ being extracted at the last time-step at the peak of loading (highlighted in red) b Training data created for the Gaussian Process (GP) using the plastic corrector, along with the GP's predicted $e^p$ for the last time-step</figcaption></figure>

<figure><img src="image_28.png" /><figcaption>**Fig. 18** Validation of the Gaussian process regression: Comparison between the results obtained via the plastic correction algorithm (pure Neuber-type algorithm) and the predicted results via the Gaussian process regression</figcaption></figure>

## 5 Neural plastic corrector: learning from examples

The aim of this section is to investigate the relative merits of the previously described Neuber-type plastic corrector that corrects elasto-static simulations and that of an AI-

based alternative that learns plastic corrections from full elasto-plastic simulations, before attempting to merge the two approaches.

The application of machine learning as a surrogate or correction technique for finite element analysis (FEA) has been successful across various fields, including the acceleration of nonlinear computations in porous specimens [30], structural analysis [37, 38] and other nonlinear mechanics problems [39, 40]. However, these methods tend to perform poorly when applied to cases that are geometrically very different from the cases in the training set [29]. Taking this knowledge on-board, we restrict our AI approach using Convolutional Neural Networks (CNNs) to the correction of elasto-plastic fields in specimens with given nominal geometries but with random geometrical defects, the morphology and positioning statistics of which are known in advance. We will concentrate here on a tensile specimen including uniformly distributed spherical pores of fixed radius.

The AI strategy is the following: we voxelise the neighbourhood of points for which plastic correction is to be performed. Next, we provide either (i) the von Mises stress from the elastic computation or (ii) the result of the proposed Neuber-type plastic corrector as input to a CNN designed to predict a pointwise output of interest stemming from an elasto-plastic simulation. Training is done by acquiring reference quantities of interests (QoI) from a series of full elasto-plastic simulations using the Z-Set suite [32], for a sufficiently large number of specimens with random pore placement.

<table><caption>**Table 4** CPU times (in seconds) for QoI computation for the test case 2: mesh of a specimen with spherical pores (661771 quadrature points, 1000$^{th}$ time-step). The training set comprises 150 points integrated over all 1000 time-steps</caption><thead><tr><td>Operation</td><td>Computing training set</td><td>Training</td><td>QoI computation</td><td>Total CPU time</td></tr></thead><tbody><tr><td>Plastic corrector</td><td>-</td><td>-</td><td>156</td><td>156</td></tr><tr><td>Gaussian process regression</td><td>0.6</td><td>0.3</td><td>1.7</td><td>2.6</td></tr></tbody></table>
<figure><img src="image_29.png" /><figcaption>**Fig. 19** Obtaining training data from a mesh for the Convolutional Neural Network</figcaption></figure>

## 5.1 Dataset generation

The objective here is to generate a robust dataset that contains sufficient mechanical information about pore-surface and pore-pore interactions, so that the QoI in specimens containing a new distribution of spherical random pores is well predicted.

The Convolutional Neural Network (CNN) used here is designed to work with 3D images. A procedure was developed in order to convert finite element results into uniform grids of voxels (see Fig. 19a–b). CNNs are resource-heavy when working with 3D images. Therefore, a size restriction is necessary for the images to fit in memory. We choose the CNN prediction of the QoI at a given voxel to be conditioned on the surrounding sub-volume of 16x16x16 voxels of the low-fidelity simulation. This sub-volume is assumed to contain enough information about the mechanical behaviour of the voxel’s surroundings. In other words, the influence of the mechanical state beyond this sub-volume on the center point is assumed to be sufficiently well-represented by the mechanical fields that will be provided as input channel for the 16x16x16 input volume. The QoI in the subvolumes is normalised by the maximum value of the QoI across the training set. Empty space due to the pores or due to regions beyond the specimen’s surfaces are encoded with a negative value, as input to our CNN is positive. If a negative input is to

be given, an additional binary channel can be used to encode the geometric information of the specimens [29].

Next, the QoI computed by the plastic corrector is extracted specifically for the points that experience plastic deformation (i.e. all points with the von Mises stress exceeding $\sigma_y$). A K-medoids clustering algorithm (scikit-learn) is used to extract clusters from these points. K-medoids is preferred because it chooses actual data points as cluster centers, rather than creating new, synthetic points like K-means. This approach ensures that each cluster center is a real sample from the dataset, and is necessary as the mechanical field around the cluster center is needed. The number of clusters per realisation of a random porous specimen was chosen to be 200, to extract a reasonable amount of data per specimen. An example of the cluster points in a mesh is shown in Fig. 19c.

Finally, an input-output pair consists of a 16x16x16 volume of real values (the von Mises stress from the elastic prediction, or the actual QoI computed using the Neuber-type plastic corrector), associated with a reference value for the QoI corresponding to the center of this volume.

## 5.2 Architecture

The input of the CNN (Fig. 20) consists of 16 x 16 x 16 subvolumes of low-fidelity information around a voxel of the QoI to be corrected. The convolutional blocks used consist of
<figure><img src="image_30.png" /><figcaption>**Fig. 20** Architecture of the CNN developed for multi-fidelity corrections of a quantity of interest (QoI)</figcaption></figure>

3D convolutional operations with zero padding and a stride of two voxels, which reduces the size of the volumes after each block, followed by ReLU activation functions. After two blocks, the volumes are flattened and reduced through a fully connected layer followed by a ReLU activation function, and a last fully connected layer reduces the dimensions of the pseudo-output to unity. An exponential function can be used on the pseudo-output, to ensure that the CNN predictions remain positive, which is activated if the QoI is required to be strictly positive. The model parameters consist of the weights of the kernel and the fully connected layers. The output of the network is the value of the predicted QoI at the centre of the given subvolume. A Mean Squared Error between the CNN predictions and the reference data summed over the batch (size of 32) is used as the loss function for training.

## 5.3 Numerical results

For the numerical study, the intrinsic dissipated energy $\phi$ integrated over a loading cycle from time $t_i$ to $t_f$ [41] is considered as a QoI to be corrected via a CNN:

$$ \phi = \int_{t_i}^{t_f} \left( f_y(\mathcal{X}; \mathcal{X}, p) + \sigma_y + \frac{(R(p))^2}{2Q} + \frac{D}{2C} (\mathcal{J}(\mathcal{X}))^2 \right) \dot{p} dt \quad (49) $$

$$ \text{with } \mathcal{J}(\mathcal{X}) = \sqrt{\frac{3}{2} \mathcal{X} : \mathcal{X}}. $$

### Impact of mechanical inputs on the performance of the CNN

The CNN is trained using the von Mises stress from the elasto-static simulation ($\bar{\sigma}_{\text{VM}}^{\#}$) as input. The CNN learns to predict the reference dissipation $\phi$. We evaluate the accuracy of the CNN on only the clusters of the K-medoids algorithm, as we expect that this reflects well the CNN accuracy on full fields due to the good representative nature of the clusters. Figure 21 shows the CNN predictions on the clusters of a mesh that is not in the training set after being trained on a progressively increasing number (20 to 100) of meshes. The results show that the accuracy of the prediction increases as a function of the amount of training data used. These results

highlight that the CNN is able to predict the QoI using only $\bar{\sigma}_{\text{VM}}^{\#}$ as input.

Next, the CNN is trained using the plastic corrector dissipation ($\hat{\phi}$) as input to the CNN. The CNN learns to predict the difference between the plastic corrector dissipation $\hat{\phi}_i$ and the reference dissipation $\phi$. After the prediction, this difference is then subtracted from the plastic corrector dissipation $\hat{\phi}_i$ to get the final CNN-predicted dissipation. We focus on learning the difference because it typically yields better results than directly learning the QoI [42].

Figure 22a shows the plastic corrector dissipation as a function of the reference dissipation, and Fig. 22b-e show the CNN predictions on the K-medoid clusters of a mesh that is not in the training set after being trained on a progressively increasing number (20 to 380) of meshes. The results show the accuracy of the predictions improved as a function of the amount of training data used. The CNN trained with 380 meshes has 89% of cluster points of a mesh, not included in the training set, falling in a ±20% error cone, and the highest values of dissipation, which are generally of the most interest, have less than a percent of error.

Figure 22f shows the mean square error (MSE) calculated on all the clusters of the mesh that is not in the training set, as a function of the number of meshes used for training. The downward trend of the CNNs trained with the two inputs is a clear indicator that the CNN is able to perform better if given more data during training.

We observe that training the CNN with the plastic corrector is beneficial when there is a low amount of available training data. The percentage of cluster points of a mesh not included in the training set, falling in a ±20% error cone, for the CNN trained on different amounts of training data, is shown in the second column of Table 5. The CNN trained with $\bar{\sigma}_{\text{VM}}^{\#}$ as input, on 20 meshes, performs significantly worse than the CNN trained with $\hat{\phi}$ as input. The difference of the accuracy of the predictions between the CNNs, however, becomes smaller and smaller as the amount of training data increases. In our opinion, this is because the corrections provided by the proposed Neuber-type approach are purely local. The plastic corrections do not contain any additional information as compared to elastic results about the source
<figure><img src="image_31.png" /><figcaption>**Fig. 21** Using von Mises stress coming from elasto-static simulations ($\bar{\sigma}_{\text{VM}}^{\#}$) as input: CNN predictions of dissipation for the clusters of a new, unseen mesh, after training for 2000 epochs on a varying number of meshes</figcaption></figure>

of the discrepancy between the input and output of the CNN model from the topology of the neighbourhood of a point.

### Accuracy of the CNN predictions for a non-spherical defect

We now illustrate the behaviour of the CNN, trained successfully on a data-set consisting of 380 specimens with spherical pores, when making predictions on a specimen containing a defect of non-spherical morphology. To this hand, we will generate a new specimen with a single non spherical defect whose geometry is illustrated in Fig. 23a. The geometric morphology of this defect comes from a computed tomography scan of a cast aluminium alloy. We ask the CNN to perform the plastic correction on all the points of the mesh that plastify (shown by the red points in Fig. 23a). After being corrected by CNN, the number of predictions that are within the 20% error cone reaches 40.4%, starting from 40.8% for correction-less Neuber-type predictions, as shown in Fig. 23. This behavior is in agreement with previous scientific findings [29] that report a poor performance of neural networks when used to make predictions for cases outside the distribution of the training set.

## 6 Conclusion

A plastic correction algorithm has been proposed that rapidly post-processes elasto-statics simulations to approximate the full-field elasto-plastic response of structures subjected to proportional loading. The classical Neuber's rule has been modified to operate on the deviatoric part of stress and strain tensors. Cyclic loading is handled via a change of origin at every loading peak, using the method in [14]. A pointwise rule of proportional evolution of the deviatoric stress and strain tensors is used to reduce J2 elasto-plastic constitutive laws to a set of equations operating on scalar representations of stresses and strains. A fully implicit time integration algorithm for this system of equations has been developed, which leads to approximations of elasto-plastic solutions for arbitrary (proportional) loading sequences.

Numerical investigations on two test cases, one being a notched structure and the other one being a specimen with spherical pores, have shown good overall approximations of cumulative plastic strain fields, for both monotonic and cyclic loading, with errors in the notched structure matching previous studies [15, 36]. The overly stiff behavior observed in these plastic corrector predictions arises because plastic accommodation and stress redistribution effects have not been taken into account [10, 15]. For cyclic loading sequences, quantities computed over later cycles, like the cumulative plastic strain range, are reasonably well approximated. The plastic corrector's accuracy improves when the plasticity is relatively confined.

Furthermore, we have attempted to isolate the error due to the local proportional evolution rule in the plastic corrector algorithm. Investigations reveal that the proportionality rule contributes significantly less error compared to the Neuber-type rule. We have also shown that the major part of the error owed to the rule of local proportional evolution occurs during load reversal, when there is no accumulating plasticity.

Despite the plastic corrector's low computational cost compared to a full nonlinear finite element simulation, time-integration over long time histories may remain expensive. However, as a consequence of the rules underlying the proposed plastic corrector (J2 plasticity, pointwise proportionality of deviatoric strain and stress tensors, scalar Neuber rule), any output of the plastic correction algorithm depends on the elastic finite element simulation via the von Mises stress only. We have shown that this property can be exploited to build 1D meta-models that accelerate the plastic correction method even further, for any quantity of interest, by sampling the output of the plastic corrector for a relatively small number of von Mises stresses, training the meta-model and using it at almost no cost in place of the time integration algorithm.

Lastly, we have investigated the possibility of performing plastic correction using deep learning. These neural networks need reference elasto-plastic computations on a set of geometries for training, and do not generalise well to other types of geometries. We have shown that the output of the Neuber-
<figure><img src="image_32.png" /><figcaption>**Fig. 22** Using plastic corrector dissipation ($\hat{\phi}$) as input: CNN predictions of dissipation in the clusters of a new, unseen mesh: **a** Plastic corrector, before passing into CNN, **b–e** CNN predictions after training for 2000 epochs on a varying number of meshes, and **f** Sensitivity analysis showing decrease of the MSE between the CNN predictions and reference data with increase in training data</figcaption></figure>

<table><caption>**Table 5** Percentage of cluster points in a mesh not included in the training set falling in a ±20% error cone for the networks trained with different inputs ($\bar{\sigma}_{VM}^{#}$ and $\hat{\phi}$) at different amounts of training data used</caption><thead><tr><td>Training data used (number of meshes)</td><td>Input σ&#773;<sup>#</sup><sub>VM</sub> (% points)</td><td>Input φ&#773; (% points)</td></tr></thead><tbody><tr><td>20</td><td>62%</td><td>80.5%</td></tr><tr><td>40</td><td>73%</td><td>77.5%</td></tr><tr><td>60</td><td>69%</td><td>76.0%</td></tr><tr><td>80</td><td>76.5%</td><td>79.0%</td></tr><tr><td>100</td><td>89%</td><td>87.5%</td></tr><tr><td>380</td><td>-</td><td>89.0%</td></tr></tbody></table>

type plastic corrector can be improved by using a CNN that leverages information about the local topology of material point neighbourhoods. We have also shown that the output of our Neuber-type methodology, as full-field input to the CNN, is indeed beneficial compared to using inputs of the NPC with a lesser mechanical content, i.e. the von Mises stress from elasto-static simulations. We have shown that the benefit rapidly becomes negligible with an increase in the training data. Thus, the plastic corrector offers an advantage

for deep learning-based corrections, particularly when the amount of training data is limited.

To further advance the research presented in this paper, several directions can be pursued. First, the accuracy of the plastic corrector could be evaluated on cyclic loading sequences with non-zero mean values and irregular shapes. Second, the plastic corrector algorithm could be extended to accommodate non-proportional applied loading conditions, with an assessment of the algorithm's robustness and limitations under these scenarios. Finally, the potential of
<figure><img src="image_33.png" /><figcaption>**Fig. 23** Accuracy of the CNN trained on 380 specimens with spherical pores on a specimen containing a defect of non-spherical morphology. The highlighted red points indicate points that undergo plastic deformation, that are selected for correction, in a mesh not included in the training set with a defect of complex morphology</figcaption></figure>

the plastic corrector in predicting elasto-viscoplastic material behavior and creep could be explored by incorporating time-dependent plasticity and viscous flow effects into the algorithm.

$$se = f^2 \qquad (50)$$

The plastic corrector equations derived in Sect. 2.4 are reminded here for the sake of readability.

## Appendix A: Algorithm for solving the plastic corrector equations

This section details the algorithm used to solve numerically the equations of the local plastic corrector, for arbitrary (proportional) loading functions $f(t)$.

### Monotonic loading

In that case, the proposed Neuber-type rule (from Sect. 2.3) reads as:

$$\text{Stress strain relation } s = (e - e^p) \qquad (51)$$

$$\text{Yield surface } \hat{f}_y(s; x, \hat{p}) = \hat{\mathcal{J}}(s, x) - \sigma_y - R(\hat{p}) \qquad (52)$$

$$\text{Yield surface evolution } \hat{f}_y \dot{\hat{p}} = 0 \text{ and } \hat{f}_y(s; x, \hat{p}) \le 0 \qquad (53)$$

$$\text{von Mises stress } \hat{\mathcal{J}}(s, x) = \left| s - \frac{x}{2\mu} \right| \bar{\sigma}_{\text{VM}}^{\#} \qquad (54)$$

$$\text{Kinematic hardening } \dot{x} = \frac{2}{3} C \dot{e}^p - D x \dot{\hat{p}} \qquad (55)$$

$$\text{Isotropic hardening } R(\hat{p}) = Q(1 - e^{-b\hat{p}}) \qquad (56)$$
$$ \text{Cumulative plastic strain } \hat{p} = \frac{1}{3\mu} |\dot{e}^p| \bar{\sigma}_{\text{VM}}^{\#} \quad (57) $$

Implicit time integration of this set of equations is done by introducing the following set of discretised equations for $s_{i+1}$, $e_{i+1}$, $e_{i+1}^p$, $\hat{p}_{i+1}$ and $x_{i+1}$ at time $t_{i+1} \in [0, T]$, scalar quantities $e_i^p$, $\hat{p}_i$ and $x_i$ at time $0 \le t_i < t_{i+1}$ being known.

$$ s_{i+1}e_{i+1} = f_{i+1}^2 \quad (58) $$

$$ s_{i+1} = e_{i+1} - e_{i+1}^p \quad (59) $$

$$ \hat{p}_{i+1} = \hat{p}_i + |e_{i+1}^p - e_i^p| \frac{1}{3\mu} \bar{\sigma}_{\text{VM}}^{\#} \quad (60) $$

$$ x_{i+1} = \frac{x_i + \frac{2}{3}C(e_{i+1}^p - e_i^p)}{1 + D(e_{i+1}^p - e_i^p)} \quad (61) $$

$$ \hat{f}_y(s_{i+1}, x_{i+1}, \hat{p}_{i+1}) \le 0 \\ \hat{f}_y(s_{i+1}, x_{i+1}, \hat{p}_{i+1})(\hat{p}_{i+1} - \hat{p}_i) = 0 \quad (62) $$

To solve this system of equations, we first need to express $s_{i+1}$ as a function of $e_{i+1}^p$ by making use of equations (58), and (59):

$$ s_{i+1} = \begin{cases} \frac{-e_{i+1}^p + \sqrt{(e_{i+1}^p)^2 + 4f_{i+1}^2}}{2} & \text{if } f_{i+1} > f_i \\ \frac{-e_{i+1}^p + \sqrt{(e_{i+1}^p)^2 + 4f_{i+1}^2}}{2} & \text{if } f_{i+1} < f_i \end{cases} \quad (63) $$

There are two roots for $s$. Substituting the expression for $e$ from the scalar stress-strain relation of the material law in the scalar Neuber-type equation, a second order polynomial equation is obtained: The positive sign in the root is taken for increasing $f$ (tensile loading), the negative sign is taken in the case of decreasing $f$ (compression loading).

Now, the following solution algorithm is proposed

* compute

$$ s_{i+1}^* = \begin{cases} \frac{-e_i^p + \sqrt{(e_i^p)^2 + 4f_{i+1}^2}}{2} & \text{if } f_{i+1} > f_i \\ \frac{-e_i^p - \sqrt{(e_i^p)^2 + 4f_{i+1}^2}}{2} & \text{if } f_{i+1} < f_i \end{cases} \quad (64) $$

* compute $f_{y,i+1}^* = \hat{f}_y(s_{i+1}^*, x_i, \hat{p}_i)$, i.e. the value of the yield function assuming that no plastic flow takes place between $t_i$ and $t_{i+1}$.

* if $f_{y,i+1}^* \le 0$, set $e_{i+1}^p = e_i^p$

* if $f_{y,i+1}^* > 0$, find $e_{i+1}^p$ such that

$$ \hat{f}_y(s_{i+1}, x_{i+1}, \hat{p}_{i+1}) = 0 \quad (65) $$

with $s_{i+1}$ given as a function of $e_{i+1}^p$ in equation (63), $x_{i+1}$ given as a function of $e_{i+1}^p$ in equation (61) and $\hat{p}_{i+1}$ given as a function of $e_{i+1}^p$ in equation (60). The

root of this equation is found by a Newton algorithm, which is initialised by setting $e_{i+1}^p = e_i^p$, $x_{i+1} = x_i$ and $\hat{p}_{i+1} = p_i$. In our implementation, the derivative of $\hat{f}_y$ with respect to $e_{i+1}^p$ is computed by finite differences.

* set $x_{i+1} = \frac{x_i + \frac{2}{3}C(e_{i+1}^p - e_i^p)}{1+D(e_{i+1}^p - e_i^p)}$ and $\hat{p}_{i+1} = \hat{p}_i + |e_{i+1}^p - e_i^p| \frac{1}{3\mu} \bar{\sigma}_{\text{VM}}^{\#}$

### Cyclic loading

In that case, the proposed Neuber-type rule (from Sect. 2.3) reads as:

$$ (s - s_o)(e - e_o) = (f - f_o)^2 \quad (66) $$

From the plastic corrector equations derived in Sect. 2.4, the scalar stress-strain relation extended for cyclic loading is:

$$ (s - s_o) = (e - e_o) - (e^p - e_o^p) \quad (67) $$

The rest of the plastic corrector equations are the same as equations (52)-(57). The equations (66) and (67) involve the quantities $s_o$, $e_o$, $e_o^p$ and $f_o$ which are updated with $s_i$, $e_i$, $e_i^p$ and $f_i$ respectively, each time a load reversal occurs. The time discretisation of these equations is given here:

$$ (s_{i+1} - s_o)(e_{i+1} - e_o) = (f_{i+1} - f_o)^2 \quad (68) $$

$$ s_{i+1} - s_o = (e_{i+1} - e_o) - (e_{i+1}^p - e_o^p) \quad (69) $$

And the expression for $s_{i+1}$ as a function of $e_{i+1}^p$ by making use of equations (68), and (69):

$$ s_{i+1} - s_o = \begin{cases} \frac{-(e_{i+1}^p - e_o^p) + \sqrt{((e_{i+1}^p - e_o^p)^2 + 4(f_{i+1} - f_o)^2)}}{2} & \text{if } f_{i+1} > f_i \\ \frac{-(e_{i+1}^p - e_o^p) - \sqrt{((e_{i+1}^p - e_o^p)^2 + 4(f_{i+1} - f_o)^2)}}{2} & \text{if } f_{i+1} < f_i \end{cases} \quad (70) $$

The same implicit time integrator used for the monotonic case presented in the previous paragraph is used to find the value of $e^p$ point-wise at each time-step.

### Python implementation

A Python implementation of this algorithm is made available under GPL licence. The code and detailed usage instructions can be found in this GitHub repository (https://github.com/AbhishekPalchoudhary/PlasticCorrector/tree/main). The repository includes scripts and instructions for point-wise plastic correction and correction of full elastic FEA computations. Examples are also included. This will enable reproduction of our results and enable further exploration of the algorithm's capabilities.
# Appendix B: Accuracy of the plastic corrector on a specimen with pores from X-ray tomography loaded in the high-cycle fatigue regime

As a final test case, we present here a specimen containing pores (arising due to the manufacturing process) that were meshed using information from tomography of a porous AlSi7Mg0.3 alloy. This case was chosen to show the particular suitability of application of the plastic corrector for high-cycle fatigue models that require elasto-plastic fields around pores of non-spherical geometry [26, 28, 43]. The parameters of the Chaboche law are chosen according to reference [44], and are summarised in Table 6.

For the plastic corrector, an elasto-static FEA computation is computed with prescribed displacements $\bar{u}_a = [u_x, 0, 0]$ for $x \ge L_c$ and $\bar{u}_a = [-u_x, 0, 0]$ for $x \le -L_c$ such that $\bar{\sigma}_{VM}^{\#}$ in the gauge section away from the pores is at the yield stress of the material (this computation corresponds to $f = 1$). The function $f(t)$ is chosen to have 20 cycles. At the peak of the cyclic load, $f(t)\bar{\sigma}_{VM}^{\#}$ is chosen to reach 47% of the yield

stress of the material in the gauge section away from pores. The boundary conditions are shown in Fig. 24.

For the reference elasto-plastic computation, prescribed displacements are applied to both ends of the porous specimen in the same way. The von Mises stress in the gauge section away from pores reaches 47% of the material's yield stress at peak loads in the cycles, assuming the body behaves elastically.

The plastic corrector and reference computations for the cumulative plastic strain range $\Delta p$ in the 20th cycle for this mesh are respectively shown in Fig. 25a and b. A scatter plot comparing $\Delta p$ in all the integration points in the mesh is shown in Fig. 25c. The accuracy of the full-field approximation indicates that the plastic corrector performs well for varying material parameters and on different sizes and geometries of pores. Despite the nominal loading being $0.47\sigma_y$, there are regions with higher stress concentrations arising due to a maximum stress concentration factor of $k_t \sim 4.3$. The plastic corrector approximates the solution well in these regions despite the high $k_t$ as

<table><caption>**Table 6** Parameters of the elasto-plastic model detailed in Eqs. (11)–(18) [31, 45]</caption><thead><tr><th>E<br>MPa</th><th>σ<sub>y</sub><br>MPa</th><th>b</th><th>Q<br>MPa</th><th>C<br>MPa</th><th>D</th></tr></thead><tbody><tr><td>75500</td><td>170</td><td>19</td><td>20</td><td>127499</td><td>1334</td></tr></tbody></table>

<figure><img src="image_34.png" /><figcaption>**Fig. 24** Boundary conditions (shown in red) for a specimen containing a sub-volume of tomography-informed pores, showing where displacement is applied to get a cyclic loading in the gauge section of the specimen (away from pores) with peak equal to 47% of the yield stress of the material. The 20th cycle is chosen for the computation of $\Delta p$</figcaption></figure>
<figure><img src="image_35.png" /><figcaption>**Fig. 25** **a–b** A comparison between $\Delta p$ in the 20th cycle calculated via the plastic corrector and a reference computation via Z-Set [32] in a few pores of a specimen containing a subvolume of pores (with the maximum stress concentration factor being $k_t \sim 4.3$). The loading corresponds to 80 MPa in the gauge section at the peak of cyclic loading, away from pores (around 47% of $\sigma_y$). **c** A scatter plot comparing $\Delta p$ calculated via the plastic corrector and a reference computation via Z-Set [32] in all the integration points of the specimen containing the subvolume of pores</figcaption></figure>

the plasticity remains confined due to the small size of the pores.

**Acknowledgements** The authors would like to thank the French National Research Agency (ANR-20-THIA-0022) for the partial funding of this project. The direct financial support of Mines Paris - PSL is also acknowledged.

**Funding** Open access funding provided by Mines Paris - PSL.

**Declarations**

**Conflict of 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.

**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.
References

1. Kuguel R (1961) A relation between theoretical stress concentration factor and fatigue notch factor deduced from the concept of highly stressed volume. Proc ASTM 61:732-748

2. Morel F, Palin-Luc T (2002) A non-local theory applied to high-cycle multiaxial fatigue. Fatigue Fracture Eng Mater Struct 25:649-665. https://doi.org/10.1046/j.1460-2695.2002.00527.x

3. Van Hooreweder B, Moens D, Boonen R, Sas P (2012) Fatigue strength analysis of notched aluminium specimens using the highly stressed volume method. Fatigue Fracture Eng Mater Struct 35(2):154-159. https://doi.org/10.1111/j.1460-2695.2011.01602.x

4. Krzyżak D, Łagoda T (2014) Fatigue life estimation of notched elements with use of non-local volumetric method. Int J Fatigue 61:59-66. https://doi.org/10.1016/j.ijfatigue.2013.12.004

5. El Khoukhi D, Morel F, Saintier N, Bellett D, Osmond P, Le V-D, Adrien J (2019) Experimental investigation of the size effect in high cycle fatigue: role of the defect population in cast aluminium alloys. Int J Fatigue 129:105222. https://doi.org/10.1016/j.ijfatigue.2019.105222

6. Dezecot S, Maurel V, Buffiere J-Y, Szmytka F, Koster A (2017) 3d characterization and modeling of low cycle fatigue damage mechanisms at high temperature in a cast aluminum alloy. Acta Mater 123:24-34. https://doi.org/10.1016/j.actamat.2016.10.028

7. Paul SK, Majumdar S, Kundu S (2014) Low cycle fatigue behavior of thermo-mechanically treated rebar. Mater Des 58:402-411. https://doi.org/10.1016/j.matdes.2014.01.079

8. Branco R, Prates PA, Costa JD, Berto F, Kotousov A (2018) New methodology of fatigue life evaluation for multiaxially loaded notched components based on two uniaxial strain-controlled tests. Int J Fatigue 111:308-320. https://doi.org/10.1016/j.ijfatigue.2018.02.027

9. Herbland, T (2009) Une méthode de correction élastoplastique pour le calcul en fatigue des zones de concentration de contraintes sous chargement cyclique multiaxial non proportionnel. Theses, École Nationale Supérieure des Mines de Paris (December). https://pastel.hal.science/tel-00479991

10. Chouman M, Gaubert A, Chaboche JL, Kanoute P, Cailletaud G, Quilici S (2014) Elastic-viscoplastic notch correction methods. Int J Solids Struct 51(18):3025-3041. https://doi.org/10.1016/j.ijsolstr.2014.04.017

11. Levieil B, Doudard C, Thevenet D, Bridier F, Ezanno A, Calloch S (2019) An original simplified method based on the use of an adjustable localization operator for low-cycle fatigue life predictions in the case of confined plasticity. Theoret Appl Fract Mech 104:102383. https://doi.org/10.1016/j.tafmec.2019.102383

12. Neuber H (1961) Theory of stress concentration for shear-strained prismatical bodies with arbitrary nonlinear stress-strain law. J Appl Mech 28(4):544-550. https://doi.org/10.1115/1.3641780 (https://asmedigitalcollection.asme.org/appliedmechanics/article-pdf/28/4/544/5444285/544_1.pdf)

13. Topper T, Wetzel R, Morrow J (1967) Neuber's rule applied to fatigue of notched specimens. J Mater ASTM 4:21

14. M, C, JP, C (1985) Adaptation of neuber's theory to stress concentration in viscoplasticity. La Recherche Aérospatiale 4, 33-40

15. Molski K, Glinka G (1981) A method of elastic-plastic stress and strain calculation at a notch root. Mater Sci Eng 50:93-100. https://doi.org/10.1016/0025-5416(81)90089-6

16. Buczynski A, Glinka G (2003): An analysis of elasto-plastic strains and stresses in notched bodies subjected to cyclic non-proportional loading paths. In: Carpinteri, A., de Freitas, M., Spagnoli, A. (eds.) Biaxial/Multiaxial Fatigue and Fracture. European Structural Integrity Society, vol. 31, pp. 265-283. Elsevier, Amsterdam, Netherlands. https://doi.org/10.1016/S1566-1369(03)80015-8

17. Ince A, Glinka G (2013) A numerical method for elasto-plastic notch-root stress-strain analysis. J Strain Anal Eng Des 48(4):229-244. https://doi.org/10.1177/0309324713477638

18. Moftakhar A, Buczynski A, Glinka G (1994) Calculation of elasto-plastic strains and stresses in notches under multiaxial loading. Int J Fract 70:357-373. https://doi.org/10.1007/BF00032453

19. Desmorat R (2002) Fast estimation of localized plasticity and damage by energetic methods. Int J Solids Struct 39(12):3289-3310. https://doi.org/10.1016/S0020-7683(02)00002-1

20. Ye D, Hertel O, Vormwald M (2008) A unified expression of elastic-plastic notch stress-strain calculation in bodies subjected to multiaxial cyclic loading. Int J Solids Struct 45(24):6177-6189. https://doi.org/10.1016/j.ijsolstr.2008.07.012

21. McDonald RJ, Socie DF (2011) A technique to estimate the local multiaxial elastic-plastic behavior from a purely elastic solution. Eng Fract Mech 78(8):1696-1704. https://doi.org/10.1016/j.engfracmech.2010.12.001

22. Taylor D (2005) Analysis of fatigue failures in components using the theory of critical distances. Eng Fail Anal 12(6):906-914. https://doi.org/10.1016/j.engfailanal.2004.12.007

23. Beremin FM, Pineau A, Mudry F et al (1983) A local criterion for cleavage fracture of a nuclear pressure vessel steel. Metall Trans A 14:2277-2287. https://doi.org/10.1007/BF02663302

24. Zok FW (2017) On weakest link theory and weibull statistics. J Am Ceram Soc 100(4):1265-1268. https://doi.org/10.1111/jace.14665

25. Rasmussen CE, Williams CK (2006) Gaussian processes for machine learning. MIT Press, Cambridge

26. Bercelli L, Moyne S, Dhondt M, Doudard C, Calloch S, Beaudet J (2021) A probabilistic approach for high cycle fatigue of wire and arc additive manufactured parts taking into account process-induced pores. Addit Manuf 42:101989. https://doi.org/10.1016/j.addma.2021.101989

27. Merot P, Morel F, Robert C, Pessard E, Gallegos Mayorga L, Buttin P (2024) Non local multiaxial fatigue modeling of defects: a unified approach to interpret size and shape effects. Theoret Appl Fract Mech 131:104378. https://doi.org/10.1016/j.tafmec.2024.104378

28. Palchoudhary A, Maurel V, Ovalle C, Kerfriden P (2024) A multiscale model for fast HCF characterization of alloys in the presence of process-induced defects. In: 16ème Colloque National en Calcul de Structures, Giens, France. CNRS, CSMA, ENS Paris-Saclay, CentraleSupélec. https://hal.science/hal-04610882

29. Krokos V, Bui Xuan V, Bordas SPA, Young P, Kerfriden P (2021) A bayesian multiscale cnn framework to predict local stress fields in structures with microscale features. Comput Mech 69(3):733-766. https://doi.org/10.1007/s00466-021-02112-3

30. Krokos V, Bordas SPA, Kerfriden P (2024) A graph-based probabilistic geometric deep learning framework with online enforcement of physical constraints to predict the criticality of defects in porous materials. Int J Solids Struct 286:112545. https://doi.org/10.1016/j.ijsolstr.2023.112545

31. Chaboche JL (1989) Constitutive equations for cyclic plasticity and cyclic viscoplasticity. Int J Plast 5(3):247-302. https://doi.org/10.1016/0749-6419(89)90015-6

32. Besson J, Leriche R, Foerch R, Cailletaud G (2012) Object-oriented programming applied to the finite element method part ii. Application to material behaviors. Revue Européenne des Éléments 7:567-588. https://doi.org/10.1080/12506559.1998.10511322

33. Geuzaine, C., Remacle, J.-F.: Gmsh. http://gmsh.info/

34. Baratta IA, Dean JP, Dokken JS, Habera M, Hale JS, Richardson CN, Rognes ME, Scroggs MW, Sime N, Wells GN (2023) DOLFINx: The next generation FEniCS problem solving environment. preprint. https://doi.org/10.5281/zenodo.10447666

35. Schijve J (2009) Fatigue of structures and materials, vol 623, 2nd edn. Springer, Dordrecht

36. Jones R, Science D, (Australia) TO, Aeronautical (Australia), MRL (1998) Stress and strain estimation at notches in aircraft structures.
DSTO Aeronautical and Maritime Research Laboratory, Melbourne (1998). http://www.dsto.defence.gov.au/corporate/reports/DSTO-GD-0196.pdf

37. Nie Z, Jiang H, Kara LB (2019) Stress field prediction in cantilevered structures using convolutional neural networks. J Comput Inf Sci Eng 20(1):630. https://doi.org/10.1115/1.4044097

38. Deshpande S, Bordas SPA, Lengiewicz J (2024) Magnet: a graph u-net architecture for mesh-based simulations. Eng Appl Artif Intell 133:108055. https://doi.org/10.1016/j.engappai.2024.108055

39. Raissi M, Perdikaris P, Karniadakis GE (2019) Physics-informed neural networks: a deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations. J Comput Phys 378:686-707. https://doi.org/10.1016/j.jcp.2018.10.045

40. Bai J, Liu G-R, Rabczuk T, Wang Y, Feng X-Q, Gu Y (2024) A robust radial point interpolation method empowered with neural network solvers (rpim-nns) for nonlinear solid mechanics. Comput Methods Appl Mech Eng 429:117159. https://doi.org/10.1016/j.cma.2024.117159

41. Besson J, Cailletaud G, Chaboche J-L, Forest S (2010) Non-linear mechanics of materials. Springer, Dordrecht

42. Chinesta F, Cueto E, Abisset-Chavanne E, Duval J-L, El Khaldi F (2020) Virtual, digital and hybrid twins: a new paradigm in data-based engineering and engineered data. Arch Comput Methods Eng 27:105-134. https://doi.org/10.1007/s11831-018-9301-4

43. Lacourt, L (2019) Étude numérique de la nocivité des défauts dans les soudures. Theses, Université Paris sciences et lettres (October). https://pastel.hal.science/tel-02512870

44. Le, V.D.: Etude de l'influence des hétérogénéités microstructurales sur la tenue en fatigue à grand nombre de cycles des alliages d'aluminium de fonderie. PhD thesis, ENSAM (2016). http://www.theses.fr/2016ENAM0012/document

45. Le V-D, Saintier N, Morel F, Bellett D, Osmond P (2018) Investigation of the effect of porosity on the high cycle fatigue behaviour of cast al-si alloy by x-ray micro-tomography. Int J Fatigue 106:24-37. https://doi.org/10.1016/j.ijfatigue.2017.09.012

**Publisher's Note** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.