# Fast mesoscopic model of plasticity in polycrystals to compute probabilistic S–N curves in high cycle fatigue

Insaf Echerradi, Daniel Weisz-Patrault, Michael Peigney

► To cite this version:

Insaf Echerradi, Daniel Weisz-Patrault, Michael Peigney. Fast mesoscopic model of plasticity in polycrystals to compute probabilistic S–N curves in high cycle fatigue. International Journal of Solids and Structures, 2025, 315 (1), pp.113348. 10.1016/j.ijsolstr.2025.113348. hal-05382373

HAL Id: hal-05382373

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

Submitted on 25 Nov 2025

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
# Fast mesoscopic model of plasticity in polycrystals to compute probabilistic S-N curves in high cycle fatigue.

Insaf Echerradi<sup>a,b</sup>, Daniel Weisz-Patrault<sup>a,*</sup>, Michael Peigney<sup>b</sup>

## Abstract

High cycle fatigue in polycrystals is mostly governed by deterministic laws such as crystal plasticity, but also depends on probabilistic properties, such as random defects and crystallographic and morphological textures, which result in significant scatter of fatigue lifetime at the macroscopic scale. Thus, modeling fatigue phenomena so that the probabilistic density function of failure is anticipated, would be useful especially for very high cycle fatigue involving up to $10^9$ cycles. To do so, the grain structure with crystal orientations should be considered in full field computations, which usually involve prohibitive computation cost therefore hindering numerical exploration of statistical distribution of fatigue failures.

This paper therefore consists of developing a very fast full field mesoscopic model of polycrystals subjected to crystal plasticity during cyclic loading based on energy minimization techniques. As a result, the uniform plastic slip in each grain is obtained in the form of a relatively simple recursive formula, which guarantees short computation time even for very high cycle fatigue. The proposed approach has been validated against a classical crystal plasticity finite element model in 2D, and satisfying agreement is observed. In addition the model has been applied in combination with classical fatigue criteria to rapidly compute the fatigue lifetime and then derive probabilistic S-N curves, hence creating a substantial link between crystallographic and morphological textures on the one hand, and fatigue lifetime estimations on the other hand.

**Keywords:** High cycle fatigue, Full field computation, Mesoscopic model, Crystal plasticity, Energetic principle

## 1. Introduction

Safety analysis of many industrial processes relies on fatigue properties of materials and structures. As shown in figure 1, several fatigue regimes are classically distinguished. Ratcheting is a regime for which the elastic-plastic behavior does not stabilize leading to an accumulation of net plastic strain. Plastic shakedown is a regime for which the behavior stabilizes into a repetitive elasto-plastic loop so that the cumulative plastic strain increases although the net plastic strain

<footnote><sup>a</sup>LMS, CNRS, École Polytechnique, Institut Polytechnique de Paris, F-91128 Palaiseau, France</footnote>

<footnote><sup>b</sup>Lab Navier, Univ Gustave Eiffel, ENPC, CNRS, F-77447 Marne la Vallée, France</footnote>

<footnote>*Corresponding author: daniel.weisz-patrault@cnrs.fr</footnote>
remains constant. Elastic shakedown is similar to plastic shakedown but the stabilized steady state is purely elastic. If the applied load is sufficiently small then the behavior remains purely elastic. This paper focuses on high or very high cycle fatigue and therefore only involves plastic and elastic shakedown. While microstructure has little influence on low-cycle fatigue properties, it does play a role in fatigue properties at high or very high cycle counts, which has been observed by optical microscope measurements [1] and synchrotron micro-tomography [2, 3]. This paper aims at developing numerical tools to investigate this link between microstructure and very high cycle fatigue properties.

<figure><img src="image_2.png" /><figcaption>Figure 1: Definition of fatigue regimes (adapted from [4]).</figcaption></figure>

High cycle fatigue in polycrystals is mostly governed by deterministic laws such as crystal plasticity, but also depends on probabilistic properties, such as random defects in the material, which result in significant scatter of fatigue lifetime at the macroscopic scale. Since crystal plasticity depends on local features of polycrystals (e.g. crystal lattice orientation, grain size and shape etc.) very high cycle fatigue properties depend on the detailed grain structure (i.e., mesoscopic scale). However, even though this mesoscopic scale is meaningful for understanding fatigue phenomena [5], design and risk analysis are usually conducted for very large parts, which necessitates macroscopic fatigue criteria and lifetime estimation. However, statistically identical polycrystals at the macroscopic scale have different grain structures at the mesoscopic scale that evolve differently with respect to cyclic loading. This is responsible for a part of the material variability resulting in scatter of fatigue lifetime, which is a major issue for design and risk analysis. Indeed, classical S-N curves [6] necessitate to repeat several times fatigue experiments in order to
explore lifetime statistical distributions for all loads considered. For design purposes, S-N curves are generally defined as the average lifetime minus two standard deviations [7]. This approach not only implies significant experimental efforts, especially for very high cycle fatigue, but is also insufficient for proof of safety as a non negligible proportion of failures occur for a number of cycles below the threshold of two standard deviations (e.g., 2.5% for a normal distribution). In addition, probabilistic approaches [8] fully relying on experimental results and identification of probability density functions are difficult to perform due to high experimental costs. Indeed, the probability density function for rare events, which is essential for proof of safety, would be poorly estimated if experiments cannot be repeated a sufficient number of times.

Thus, modeling fatigue phenomena so that the probabilistic density function of failure is anticipated, would be useful especially for very high cycle fatigue involving up to $10^8$ cycles. To develop such models, the grain structure with crystal orientations should be considered in full field computations. However, the exact grain structure is unknown in most industrial applications, only statistical information is known such as grain size and shape distributions or crystal orientation distribution. Thus, to capture the probability density function of rare events, it would be necessary to perform numerous full field computations corresponding to different draws of random grain structures, which seems very difficult to achieve for 3D polycrystals due to high computation costs [9].

Within this framework, this paper aims at contributing to a general upscaling strategy of high cycle fatigue in order to estimate probabilistic fatigue lifetime at the macroscopic scale as a function of statistical descriptors of the grain structure. Such an upscaling strategy has already been applied within the context of grain growth [10, 11]. The proposed methodology consists of developing a very fast full field mesoscopic model of polycrystals subjected to crystal plasticity during cyclic loading, and then applying classical fatigue criteria to derive probabilistic S-N curves, hence creating a substantial link between crystallographic and morphological textures on the one hand, and fatigue lifetime estimations on the other hand.

Deterministic approaches to characterize the endurance limit (i.e., the stress level below which an infinite number of cycles can be applied without causing fatigue failure), such as the well known Dang Van criterion [5, 12–16], have been established at the macroscopic scale, but rely on mesoscopic theoretical bases with the concept of misoriented grain or critical plane. Furthermore, it may be noted that [17] shows that the mechanisms envisaged to explain fatigue phenomena at the microstructure scale [18, 19] involve the Dang Van criterion, thus giving it a solid micromechanical basis. These approaches present the advantage to involve a reduced number of parameters that can be easily identified with some experiments. However, at the macroscopic scale, even though these classical fatigue criteria explicitly refer to mesoscopic properties (e.g., misoriented grain or critical plane), there is no quantitative link between crystallographic and morphological textures and fatigue lifetime prediction. The main drawbacks of these models are to rely only on experimental data to consider the scatter of fatigue lifetime, and to involve new experimental investigations for each type of texture (obtained from various fabrication and forming processes).

To investigate in more details fatigue mechanisms and the relationship between the grain structure and fatigue lifetime, numerical simulations at the mesoscopic scale have been developed [20–23]. These models involve the computation of crystal plasticity in polycrystals subjected to cyclic loads with hundred of grains whose crystal lattice orientation is specified. Slip systems are
therefore identified and the cumulative plastic slip can be computed in each grain. Crystal plasticity finite elements models (CPFEM) are usually utilized, but fast Fourier transform (FFT) techniques can also be considered. For instance, $10^6$ cycles have been computed in a 2D polycrystal of two thousand grains. These full field computations only enable to consider high cycle fatigue for a small number of specific 2D polycrystalline structures due to their computational cost. Moreover, some interesting studies [22] attempted to introduce probabilistic features on the basis of full field computations, but the survival probability at *n* cycles was identified on the basis of very few different mesostructures, and for a unique number of cycles *n*.

To overcome this difficulty of large computation costs, simplified fast methodologies have been proposed such as *direct cyclic methods* which predict the asymptotic response of an elasto-plastic body subjected to cyclic loading [24–27]. Such approaches have been used to develop multiscale fatigue lifetime models [4], which are based on homogenization techniques. However, there is no reference to the detailed microstructure and the multiscale model is identified on macroscopic S-N curves. Therefore, despite their great interest, these approaches cannot be used as a replacing tool of CPFEM of FFT full field computations within the context of the envisaged upscaling strategy.

To significantly reduce computation time of full field mesoscopic computations, a general semi-analytical approach is proposed both in 3D or 2D in this contribution. A finite polycrystal composed of *N* grains of arbitrary shape is considered with uniform applied stress at the boundary. Each grain is subjected to an unknown and uniform plastic strain. On this basis, the total energy, composed of the potential energy (i.e., the elastic energy minus the work of external forces) and the dissipated energy, can be calculated analytically as a function of the unknown plastic strain in each grain. Within the framework of rate independent plastic mechanisms [28, 29], the unknown plastic strain can be determined eventually by minimizing the total incremental energy (see e.g., [30–37]). Lagrange multipliers enable to determine the plastic slip analytically in each grain. The obtained analytical solution is general and takes into account the real shape of the grains and the interactions between grains. However, this solution is inconvenient because 1) numerical integrals over arbitrary shapes are involved, which is computationally costly, and 2) the plastic strain in each grain depends in a very non-linear way of the plastic strain of all the other grains.

To overcome these difficulties and to reach very short computation time, several assumptions are introduced. In particular, integrals computation is facilitated by approximating the entire mesostructure by ellipsoids fully embedded in the grains [38]. In addition, a weakly coupled solution is derived naturally from the general fully coupled formula by introducing assumptions that are reasonable within the context of very high cycle fatigue involving low applied stress. More precisely the weak coupling consists of considering the fully decoupled solution (i.e., independent grains) as a first estimation, and then computing a correction corresponding to the effect of each grain on its neighbors, and the reciprocal effect of the corrected neighbors on each grain. As a result, the plastic slip in each grain and for each cycle is obtained as a simple recursive formula, which guarantees very short computation time. Furthermore, even though the plastic response is highly anisotropic as slip systems significantly depend on crystal orientation, it is reasonable to consider isotropic stiffness tensor to compute the elastic response despite the lattice symmetry (i.e., cubic in the proposed examples). Therefore, the stiffness tensor does not depend on crystal orientations, which facilitates the calculation of the elastic energy. Moreover, this assumption is also useful to determine which slip system is activated in each grain. Indeed, an elastic test
is easily computed under the assumption of elastic isotropy. The obtained trial resolved shear stresses according to the slip planes of each grain can be compared to the critical shear stress (i.e., a threshold above which the grain undergo plastic slip). For each grain, we assumed that only one slip system is activated. In general, this assumption is not verified as several slip mechanisms are usually activated simultaneously, but for small applied stress as in very high cycle fatigue, this assumption is well verified. For example, we present in the following a CPFEM computation for which several plastic slip systems can be activated simultaneously. However, for low applied stress a single slip system is activated in each grain. The proposed semi-analytical strategy to compute very large number of loading cycles in polycrystals undergoing crystal plasticity can be interpreted as a major extension of the well known Eshelby elliptical inclusion problem subjected to uniform eigenstrain in a infinite matrix [39, 40].

Eventually, the fatigue life of each polycrystalline structure is determined using different failure criteria, depending on whether an infinite endurance regime or regimes leading to a finite lifetime are reached.

It should be noted that stress gradient is known to have a significant influence on fatigue lifetime [17, 20, 41]. Indeed, high stress gradients correspond to small area of maximum stress, which may contain only grains with higher fatigue resistance (i.e., favorable crystal orientation, size or shape). This effect is taken into account in the proposed methodology by considering that the area where the stress field is more or less uniform depends on the stress gradient. The higher the stress gradient is, and the smaller the RVE should be, hence the smaller the number of grains N should be set.

This paper is organized as follows. First theoretical developments are detailed in section 2, which leads to the determination of the general fully coupled solution. The weakly coupled solution is then detailed. In addition, to efficiently capture infinite endurance, asymptotic formula are derived to provide the stabilized plastic slip and shear stress fields without computing successive cycles. The model validity and limitations are investigated by a comparison with CPFEM performed on FreeFEM++ [42] and presented in section 3.1. Results are presented in section 3 to analyze the model sensitivity to different aspects of the microstructure and derive probabilistic S-N curves. Conclusive remarks are given in section 4.

## 2. Model derivation

### 2.1. Energy optimization strategy

In this section, theoretical developments underlying the proposed mesoscopic model are presented. The macroscopic load is denoted by $\bar{\sigma}^j$ for the j-th loading step. Since we assumed that a single slip direction may be activated, the cyclic macroscopic stress is assumed to be proportional to a constant macroscopic stress $\bar{\sigma}$:

$$ \bar{\sigma}^j = f(j) \times \bar{\sigma} \quad (1) $$

Where $j \in \mathbb{N} \mapsto f(j) \in \mathbb{R}$ is a scalar function (e.g., oscillating around its average value). It should be noted that $f(j)$ is not necessarily a periodic function, but can represent cycles of different amplitudes.

Consider a representative volume element (RVE) denoted by $\Omega$ composed of N grains of arbitrary shape denoted by $\Omega_k$ where $1 \le k \le N$. In practice we use Voronoi-Laguerre tessellations
but other techniques could be used as well. Each grain is associated with a crystal orientation. The model is implemented for slip systems corresponding to face-centered cubic system (FCC), which covers several common metals such as copper, aluminum, nickel or austenitic steels, but the proposed methodology could be easily adapted for other lattices such as body-centered cubic systems (BCC). Thus, in each grain there are $N_n$ slip planes (e.g., $N_n = 4$ for FCC) and for each slip plane there are $N_s$ slip directions (e.g., $N_s = 3$ for FCC), hence $N_m = N_n \times N_s$ slip mechanisms. Normal vectors and slip directions are respectively denoted by $\mathbf{n}_k^m$ and $\mathbf{s}_k^m$ (where $1 \le k \le N$ and $1 \le m \le N_m$). For each grain $k$ and slip mechanism $m$ the following maximum resolved elastic trial shear stress is computed as

$$ \bar{\tau}_k^m(j) = \bar{\sigma}^j : (\mathbf{n}_k^m \otimes \mathbf{s}_k^m) \qquad (2) $$

To determine the activated slip mechanism we consider the first slip mechanism whose resolved trial shear stress given by (2) reaches the critical shear stress of grain $k$ denoted by $\tau_k^C$. Otherwise, if there is no slip mechanism so that the critical shear stress is reached, the grain $k$ remains elastic. In the following, the exponent $m$ is discarded for the sake of readability.

Since grains undergo plastic slips, the problem to solve is non-linear. To overcome this difficulty, the *unknown* plastic slip in grain $k$ denoted by $\gamma_k^p$, so that the plastic strain reads $\varepsilon_k^p = (\gamma_k^p/2) (\mathbf{n}_k \otimes \mathbf{s}_k + \mathbf{s}_k \otimes \mathbf{n}_k)$, is replaced by a *trial* and *known* eigenstrain denoted by

$$ \varepsilon_k^* = \frac{\gamma_k^*}{2} (\mathbf{n}_k \otimes \mathbf{s}_k + \mathbf{s}_k \otimes \mathbf{n}_k) \qquad (3) $$

where in the following the star exponent always relate to *known* trial quantities. The problem becomes linear for any given eigenstrain. The determination of the plastic slips $\gamma_k^p$ is then carried out through an energetic principle. Indeed, the total energy composed of the elastic potential energy $\Phi[\gamma_1^*, \dots, \gamma_N^*]$ associated with the dissipated energy $\mathcal{D}[\gamma_1^*, \dots, \gamma_N^*]$, is minimized to determine the plastic slips $(\gamma_1^p, \dots, \gamma_N^p)$, which reads:

$$ (\gamma_1^p, \dots, \gamma_N^p) = \arg \min_{(\gamma_1^*, \dots, \gamma_N^*) \in \mathbb{R}^N} \left( \Phi[\gamma_1^*, \dots, \gamma_N^*] + \mathcal{D}[\gamma_1^*, \dots, \gamma_N^*] \right) \qquad (4) $$

It should be noted that the non-linearity associated to the elastic-plastic behavior is entirely contained in the minimization problem (4).

In the following, the plastic strain increment during the $j$-th loading step is derived as an explicit analytical formula, which enables us to compute easily the cumulative plastic strain as a function of cycles. Following calculations are derived for the $j$-th loading step, knowing that the $j-1$ previous loading step has already been computed. Thus, the plastic slip increment $\Delta\gamma_k^p$ in the $k$-th grain and taking place during the $j$-th loading step has to be computed to update the total plastic slip such that

$$ \gamma_k^{p,j} = \Delta\gamma_k^p + \gamma_k^{p,j-1} \qquad (5) $$

where $\gamma_k^{p,j-1}$ is known from previous loading steps. Similarly, the trial additional plastic slip $\Delta\gamma_k^*$ during the $j$-th loading step is involved in the following calculation, and the total trial plastic slip reads as

$$ \gamma_k^* = \Delta\gamma_k^* + \gamma_k^{p,j-1} \qquad (6) $$
In addition it should be noted that the proposed minimization procedure is performed over the set of piece-wise uniform plastic strains according to the grain structure (i.e., plastic strains are uniform in each grain). This assumption is essential to enable the derivation of an analytical solution and avoid time consuming FE computation by meshing the grains.

## 2.2. Elastic potential energy

The elastic bulk energy $w_k$ in the k-th grain subjected to uniform eigenstrain $\varepsilon_k^*$ is given by

$$w_k = \frac{1}{2} \int_{\Omega_k} (\varepsilon - \varepsilon_k^*) : C : (\varepsilon - \varepsilon_k^*) d\Omega \quad (7)$$

where $\varepsilon$ is the total strain field defined in $\Omega$, $\varepsilon_k^*$ is the uniform eigenstrain defined in the grain k, and the elasticity tensor C is the same for all the grains (i.e., isotropy assumption). The elastic potential energy $\Phi$ of the entire polycrystal is defined as

$$\Phi = \sum_{k=1}^{N} w_k - \int_{\partial\Omega} (\bar{\sigma}^j \cdot n) \cdot u \, dS + \frac{1}{2} \sum_{k=1}^{N} \int_{\Omega_k} X : \varepsilon_k^* \, d\Omega \quad (8)$$

where the second term in (8) corresponds to the work of external forces and the third term corresponds to the free energy due to kinematic hardening where $X$ is the back-stress 2nd order tensor. For each grain $X(\varepsilon_k^*) = \zeta_X \varepsilon_k^*$ where $\zeta_X$ (MPa) is a kinematic hardening parameter. The displacement field $u$ and the corresponding strain field $\varepsilon$ in (8) solve the elasticity problem

$$\left\{
\begin{array}{l}
\sigma = C : (\varepsilon - \varepsilon_k^*) \text{ in } \Omega_k, (\forall k, 1 \le k \le N) \\
\text{div}(\sigma) = 0 \text{ in } \Omega \\
\sigma \cdot n = \bar{\sigma}^j \cdot n \text{ on } \partial\Omega \\
\varepsilon = \frac{1}{2} (\nabla u + \nabla u^\top)
\end{array}
\right.
\quad (9)$$

The following developments detail how to rewrite $\Phi$ as an explicit function of $\bar{\sigma}^j$ and $(\varepsilon_1^*, \dots, \varepsilon_N^*)$. Let $\chi_k$ be the characteristic function of $\Omega_k$ ($1 \le k \le N$), i.e. $\chi_k(x) = 1$ if $x \in \Omega_k$ and $\chi_k(x) = 0$ otherwise. Setting

$$\tau = C^{-1} : \bar{\sigma}^j + \sum_{k=1}^{N} \chi_k \varepsilon_k^* \quad (10)$$

and noting that $\int_{\partial\Omega}(\bar{\sigma}^j \cdot n) \cdot u \, dS = \int_{\Omega} \bar{\sigma}^j : \varepsilon \, d\Omega$, the elastic potential energy can be rewritten as

$$\Phi = \frac{1}{2} \int_{\Omega} \varepsilon : C : \varepsilon \, d\Omega - \int_{\Omega} \varepsilon : C : \tau \, d\Omega + \frac{1}{2} \sum_{k=1}^{N} |\Omega_k| \varepsilon_k^* : C : \varepsilon_k^* + \frac{1}{2} \sum_{k=1}^{N} \int_{\Omega_k} X : \varepsilon_k^* \, d\Omega \quad (11)$$

As already mentioned, for a given known eigenstrain $\varepsilon_k^*$, the considered problem is linear, and it is well known that the solution $u$ to the elasticity problem (9) minimizes the potential energy (11) among all kinematically admissible displacement fields. This implies that for any strain field $\varepsilon'$ corresponding to a kinematically admissible displacement, the following relation holds:

$$\int_{\Omega} \varepsilon : C : \varepsilon' \, d\Omega = \int_{\Omega} \tau : C : \varepsilon' \, d\Omega \quad (12)$$
Choosing the uniform strain $\varepsilon' = \int_{\Omega}(\varepsilon - \tau) d\Omega$ in (12) gives the average value of the strain $\varepsilon$:

$$ \frac{1}{|\Omega|} \int_{\Omega} \varepsilon d\Omega = \frac{1}{|\Omega|} \int_{\Omega} \tau d\Omega = C^{-1} : \bar{\sigma}^j + \sum_{k=1}^{N} c_k \varepsilon_k^* \quad (13) $$

where $c_k = |\Omega_k|/|\Omega|$ is the volume fraction of grain $k$. In addition, using (12) with $\varepsilon' = \varepsilon$ leads to

$$ \int_{\Omega} \varepsilon : C : \varepsilon d\Omega = \int_{\Omega} \tau : C : \varepsilon d\Omega \quad (14) $$

Hence by plugging in (11), we obtain

$$ \Phi = -\frac{1}{2} \int_{\Omega} \tau : C : \varepsilon d\Omega + \frac{1}{2} \sum_{k=1}^{N} |\Omega_k| \varepsilon_k^* : C : \varepsilon_k^* + \frac{1}{2} \sum_{k=1}^{N} \int_{\Omega_k} X : \varepsilon_k^* d\Omega \quad (15) $$

The first integral in (15) can be expressed as a linear function of $\bar{\sigma}^j$ and $\varepsilon_1^*, \dots, \varepsilon_N^*$. Indeed, by using definition (10) of $\tau$ and relation (13), we find:

$$ \begin{aligned} \int_{\Omega} \tau : C : \varepsilon d\Omega &= \bar{\sigma}^j : \int_{\Omega} \varepsilon d\Omega + \sum_{k=1}^{N} \varepsilon_k^* : C : \left( \int_{\Omega_k} \varepsilon d\Omega \right) \\ &= |\Omega| \bar{\sigma}^j : C^{-1} : \bar{\sigma}^j + \sum_{k=1}^{N} |\Omega_k| \bar{\sigma}^j : \varepsilon_k^* + \sum_{k=1}^{N} \varepsilon_k^* : C : \left( \int_{\Omega_k} \varepsilon d\Omega \right) \end{aligned} \quad (16) $$

The superposition principle implies that

$$ \varepsilon = C^{-1} : \bar{\sigma}^j + \sum_{k=1}^{N} \varepsilon_k \quad (17) $$

where $\varepsilon_k$ is the strain field corresponding to the solution of the problem where only a single grain $k$ is subjected to the eigenstrain $\varepsilon_k^*$ with free boundary conditions. This problem is defined for any given $k = 1, \dots, N$ as

$$ \left\{ \begin{array}{l} \sigma = C : (\varepsilon - \varepsilon_k) \text{ in } \Omega_k \\ \sigma = C : \varepsilon \text{ in } \Omega - \Omega_k \\ \operatorname{div}(\sigma) = 0 \text{ in } \Omega \\ \sigma \cdot n = 0 \text{ on } \partial\Omega \end{array} \right. \quad (18) $$

Problem (18) is linear in $\varepsilon_k^*$ so that $\varepsilon_k(x)$ can be written as

$$ \varepsilon_k(x) = S_k(x) : \varepsilon_k^* \quad (19) $$

where the Eshelby-like fourth-order tensors $S_k$ do not depend on $\varepsilon_k^*$ and are entirely determined from $\Omega, \Omega_k$ and $C$. Explicit formulae can be derived and are detailed in section 2.5 It follows that

$$ \int_{\Omega_k} \varepsilon d\Omega = |\Omega_k| \left( C^{-1} : \bar{\sigma}^j + \sum_{l=1}^{N} S_{kl} : \varepsilon_l^* \right) \quad (20) $$
where the uniform fourth-order tensors

$$ S_{kl} = \frac{1}{|\Omega_k|} \int_{\Omega_k} S_l(x) \, d\Omega \qquad (21) $$

have been introduced. The fourth order stiffness tensor $C$ is assumed to be isotropic, i.e.

$$ C = 2\mu J + k_0 \mathbf{1} \otimes \mathbf{1} \qquad (22) $$

where $\mu$ (MPa) and $k_0$ (MPa) are the shear and bulk moduli respectively. In (22) $\mathbf{1}$ is the second order identity tensor, $I$ denotes the fourth order identity tensor and $J = I - (1/3)\mathbf{1}\otimes\mathbf{1}$. Substituting (20) in (16) yields

$$ \begin{aligned} \Phi = -\frac{1}{2} |\Omega| \bar{\sigma}^j : C^{-1} : \bar{\sigma}^j - \sum_{k=1}^{N} |\Omega_k| \bar{\sigma}^j : \varepsilon_k^* - \frac{1}{2} \sum_{k=1}^{N} \sum_{l=1}^{N} |\Omega_k| \varepsilon_k^* : C : S_{kl} : \varepsilon_l^* & \\ + \frac{1}{2} \sum_{k=1}^{N} |\Omega_k| \varepsilon_k^* : C : \varepsilon_k^* + \frac{1}{2} \sum_{k=1}^{N} \int_{\Omega_k} \mathbf{X} : \varepsilon_k^* \, d\Omega & \end{aligned} \qquad (23) $$

Expression (23) is quadratic in $(\bar{\sigma}^j, \varepsilon_1^*, \dots, \varepsilon_N^*)$ as could be expected from the linearity of the elasticity problem (9). All information on the microstructure (shape and distribution of the grains) is contained in the Eshelby like tensors $S_{kl}$.

The solution to the single inclusion problem (18) satisfies

$$ \int_{\Omega} \varepsilon' : C : (\varepsilon_k - \chi_k(x) \varepsilon_k^*) \, d\Omega = 0 \qquad (24) $$

for any strain field $\varepsilon'$ deriving from a kinematically admissible displacement. Choosing in particular $\varepsilon' = \varepsilon_l$ (for $1 \le l \le N$) we obtain

$$ \int_{\Omega} \varepsilon_l : C : (\varepsilon_k - \chi_k(x) \varepsilon_k^*) \, d\Omega = 0 \qquad (25) $$

Swapping the role of $k$ and $l$ in the reasoning, we obtain in a similar fashion that

$$ \int_{\Omega} \varepsilon_k : C : (\varepsilon_l - \chi_l(x) \varepsilon_l^*) \, d\Omega = 0 \qquad (26) $$

Hence

$$ \int_{\Omega_r} \varepsilon_k^* : C : \varepsilon_l d\Omega = \int_{\Omega_l} \varepsilon_l^* : C : \varepsilon_k d\Omega \qquad (27) $$

Using (20) the following relation is obtained:

$$ |\Omega_k| \varepsilon_k^* : C : S_{kl} : \varepsilon_l^* = |\Omega_l| \varepsilon_l^* : C : S_{lk} : \varepsilon_k^* \qquad (28) $$

Recalling that only one slip system is activated and using (3), the potential energy (23) reduces to:

$$ \begin{aligned} \Phi[\Delta\gamma_1^*, \dots, \Delta\gamma_N^*] = & -\frac{1}{2} |\Omega| \bar{\sigma}^j : C^{-1} : \bar{\sigma}^j - \sum_{k=1}^{N} |\Omega_k| \bar{\tau}_k^j (\Delta\gamma_k^* + \gamma_k^{p,j-1}) \\ & + \frac{\mu}{2} \sum_{k=1}^{N} \sum_{l=1}^{N} |\Omega_k| \alpha_{kl} (\Delta\gamma_k^* + \gamma_k^{p,j-1}) (\Delta\gamma_l^* + \gamma_l^{p,j-1}) \end{aligned} \qquad (29) $$
where $\bar{\tau}_k^j$ is the resolved macroscopic shear stress and $\alpha_{kl}$ are the following dimensionless coefficients:

$$ \left\{ \begin{array}{l} \bar{\tau}_k^j = \mathbf{n}_k \cdot \bar{\sigma}^j \cdot \mathbf{s}_k \\ \alpha_{kl} = \mu \left( \delta_{kl} - \frac{1}{2}(\mathbf{n}_k \otimes \mathbf{s}_k + \mathbf{s}_k \otimes \mathbf{n}_k) : \mathbf{S}_{kl} : (\mathbf{n}_l \otimes \mathbf{s}_k + \mathbf{s}_k \otimes \mathbf{n}_l) \right) + \frac{\zeta_x}{2} \delta_{kl} \end{array} \right. \qquad (30) $$

In (30) $\delta_{kl}$ is the Kronecker symbol. We note from (28) that the coefficients $\alpha_{kl}$ satisfy the relation:

$$ |\Omega_k|\alpha_{kl} = |\Omega_l|\alpha_{lk} \qquad (31) $$

### 2.2.1. Dissipated energy

In this section, the dissipated energy $\mathcal{D}[\Delta\gamma_k^*]$ is determined. Consider the cumulative plastic strain at the j-th loading step:

$$ \gamma_k^{\text{cum},*}(\Delta\gamma_k^*) = |\Delta\gamma_k^*| + \gamma_k^{\text{cum},j-1} \qquad (32) $$

where $\gamma_k^{\text{cum},j-1}$ is the cumulative plastic strain in the k-th grain at the previous loading step. It is assumed that only one slip system is activated in each grain k, which is well validated for small applied stress as shown in the following. As shown in figure 2 a three phases hardening curve composed of a non-linear hardening phase followed by a plateau and then a linear softening is considered in this paper as suggested in [43]. It should be noted that the softening phase corresponds to rapid damage of the grain, which is completed within very few loading cycles. Therefore a linear softening phase is sufficient. The hardening phase and plateau (i.e., phase I+II) are obtained by considering a critical shear stress of the form

$$ \tau_k^C = \tau_k^0 + \mu b \sqrt{\alpha \rho_k} \qquad (33) $$

where $\tau_k^0$ (MPa) is the critical shear stress before plastic slip, b (m) is the Burger vector norm, $\alpha$ is a dimensionless hardening coefficient, and $\rho_k$ (m⁻²) is the dislocation density in the k-th grain (defined as the total dislocation length per unit volume). According to the Hall-Petch effect the initial critical shear stress $\tau_k^0$ depends on grain size as

$$ \tau_k^0 = \tau^{\min} \left( 1 + \frac{\xi}{\sqrt{D_k}} \right) \qquad (34) $$

where $D_k$ is the equivalent grain diameter, $\tau^{\min}$ (MPa) is the minimal critical shear stress (i.e., for infinite grain size), and $\xi$ is a dimensionless material coefficient. For the evolution of the dislocation density we adopt the law proposed in [44] and given by

$$ \dot{\rho}_k = \frac{\dot{\gamma}_k^{\text{cum}}}{b} \left( \frac{\sqrt{\rho_k}}{K} - 2y_c \rho_k \right) \qquad (35) $$

where K is a dimensionless coefficient, and $y_c$ (m) is an annihilation distance between dislocations. We assume that the initial plastic slip and initial dislocation density are negligible. Thus, initially $1 - 2K y_c \sqrt{\rho_k} > 0$ when $\rho_k = \rho_{\text{ini}} \approx 0$. Then $\rho_k > 0$ increases and is bounded by $\rho_{\text{max}} =$
$1/(2K y_c)^2$ as $\dot{\rho}_k = 0$ when $\rho = \rho_{max}$ according to (35). Therefore $1 - 2K y_c \sqrt{\rho_k} > 0$, and a direct integration of (35) gives:

$$ \gamma_k^{\text{cum}} = -\frac{b}{y_c} \ln(1 - 2K y_c \sqrt{\rho_k}) \quad (36) $$

By inverting (36) and injecting in (33) one obtains the following hardening law (i.e., phase I+II):

$$ \forall \gamma_k^{\text{cum}} \le \gamma_k^{\text{soft}}, \quad \tau_k^C(\gamma_k^{\text{cum}}) = \tau_k^0 + \mu \zeta (1 - \exp(-x \gamma_k^{\text{cum}})) \quad (37) $$

where $x$ and $\zeta$ are the following dimensionless coefficients:

$$ x = \frac{y_c}{b} \quad \text{and} \quad \zeta = \frac{\sqrt{\alpha}}{2Kx} \quad (38) $$

The softening phase (i.e., phase III) is assumed to be linear such as the critical shear stress reads as

$$ \forall \gamma_k^{\text{cum}} \ge \gamma_k^{\text{soft}}, \quad \tau_k^C(\gamma_k^{\text{cum}}) = \tau_k^{\text{max}} (1 - \chi (\gamma_k^{\text{cum}} - \gamma_k^{\text{soft}})) \quad (39) $$

where $\tau_k^{\text{max}} = \tau_k^C(\gamma_k^{\text{soft}})$ is obtained from (37), and $\chi$ is a dimensionless softening coefficient. It should be noted that during the softening phase the critical shear stress decreases and eventually vanishes when $\gamma_k^{\text{cum}} = \gamma_k^{\text{max}}$ where

$$ \gamma_k^{\text{max}} = \frac{1}{\chi} + \gamma_k^{\text{soft}} \quad (40) $$

Since $\tau_k^C$ continuously evolves during the hardening and softening, it is necessary to calculate the dissipated energy as an integral over each loading cycle. The infinitesimal dissipated energy per unit volume for an infinitesimal plastic slip $d\gamma$ reads $\tau_k^C d\gamma$. Therefore the total dissipated energy reads as

$$ \mathcal{D}[\Delta\gamma_1^*, \dots, \Delta\gamma_N^*] = \sum_{k=1}^{N} \int_{\Omega_k} \left[ \int_{0}^{|\Delta\gamma_k^*|} \tau_k^C d\gamma \right] d\Omega \ge 0 \quad (41) $$

For the hardening phase (i.e., phases I+II), recalling that the plastic slip is assumed uniform in each grain and plugging (32) as well as (37) into (41) one obtains:

$$ \mathcal{D}[\Delta\gamma_1^*, \dots, \Delta\gamma_N^*] = \sum_{k=1}^{N} |\Omega_k| \left[ (\tau_k^0 + \mu\zeta) |\Delta\gamma_k^*| - \frac{\mu\zeta}{x} \exp(-x \gamma_k^{\text{cum},j-1}) (1 - \exp(-x |\Delta\gamma_k^*|)) \right] \quad (42) $$

For the softening phase (i.e., phase III), plugging (32) and (39) into (41) yields

$$ \mathcal{D}[\Delta\gamma_1^*, \dots, \Delta\gamma_N^*] = \sum_{k=1}^{N} |\Omega_k| \tau_k^{\text{max}} \left( 1 + \chi \gamma_k^{\text{soft}} - \chi \gamma_k^{\text{cum},j-1} - \frac{\chi}{2} |\Delta\gamma_k^*| \right) |\Delta\gamma_k^*| \quad (43) $$
<figure><img src="image_3.png" /><figcaption>Figure 2: Three phase hardening/softening curve.</figcaption></figure>

## 2.3. Energy minimization and fully coupled solution

The stationary condition of the minimization problem (4) requires to compute the derivative of the potential and dissipated energies with respect to the plastic slips. From (29) one obtains the following derivative of the total potential energy:

$$ \frac{\partial \Phi}{\partial \Delta \gamma_k^*} = |\Omega_k| \left[ -\bar{\tau}_k^j + \mu \sum_{l=1}^{N} \alpha_{kl} (\Delta \gamma_l^* + \gamma_l^{p,j-1}) \right] \quad (44) $$

Since the absolute value $|\Delta\gamma_k^*|$ appears in the dissipated energy, one cannot consider the derivative in $\Delta\gamma_k^* = 0$. For the hardening phase (i.e., phases I+II), the derivative of the dissipated energy reads (42) as, for $\Delta\gamma_k^* \neq 0$,

$$ \frac{\partial D}{\partial \Delta \gamma_k^*} = |\Omega_k| \left[ (\tau_k^0 + \mu \zeta) - \mu \zeta \exp(-x (|\Delta \gamma_k^*| + \gamma_k^{\text{cum},j-1})) \right] \frac{\Delta \gamma_k^*}{|\Delta \gamma_k^*|} \quad (45) $$

Therefore the stationary conditions lead to the following fully coupled problem (i.e., involving all the grains in the polycrystal), for any $k$ (with $1 \le k \le N$):

$$ -\frac{\bar{\tau}_k^j}{\mu} + \sum_{l=1}^{N} \alpha_{kl} (\Delta \gamma_l^* + \gamma_l^{p,j-1}) + \left( \frac{\tau_k^0}{\mu} + \zeta - \zeta \exp(-x (|\Delta \gamma_k^p| + \gamma_k^{\text{cum},j-1})) \right) \frac{\Delta \gamma_k^p}{|\Delta \gamma_k^p|} = 0 \quad (46) $$

For the linear softening phase (i.e., phase III), from (43) one obtains:

$$ \frac{\partial D}{\partial \Delta \gamma_k^*} = |\Omega_k| \tau_k^{\text{max}} (1 + \chi \gamma_k^{\text{soft}} - \chi \gamma_k^{\text{cum},j-1} - \chi |\Delta \gamma_k^*|) \frac{\Delta \gamma_k^*}{|\Delta \gamma_k^*|} \quad (47) $$
Hence the stationary conditions lead to the following fully coupled problem (i.e., involving all the grains in the polycrystal), for any $k$ (with $1 \le k \le N$):

$$ -\frac{\bar{\tau}_k^j}{\mu} + \sum_{l=1}^{N} \alpha_{kl} (\Delta\gamma_l^* + \gamma_l^{p,j-1}) + \frac{\tau_k^{\max}}{\mu} (1 + \chi \gamma_k^{\text{soft}} - \chi \gamma_k^{\text{cum},j-1} - \chi |\Delta\gamma_k^*|) \frac{\Delta\gamma_k^*}{|\Delta\gamma_k^*|} = 0 \quad (48) $$

Using (46) and (48), the stationary condition enables to determine the plastic slip increment during the $j$-th loading step, which reads for all $k$ ($1 \le k \le N$):

$$ \Delta\gamma_k^p = \begin{cases} \Delta\gamma_k^+ & \text{if } \Delta\gamma_k^+ > 0 \\ \Delta\gamma_k^- & \text{if } \Delta\gamma_k^- < 0 \\ 0 & \text{if } \Delta\gamma_k^+ \le 0 \text{ and } \Delta\gamma_k^- \ge 0 \end{cases} \quad (49) $$

where for the hardening phase (i.e., phase I+II):

$$ \left\{ \begin{aligned} \Delta\gamma_k^+ &= -\gamma_k^{p,j-1} + \frac{\bar{\tau}_k^j - (\tau_k^0 + \mu\zeta)}{\mu\alpha_{kk}} - \eta_k \\ &\quad + \frac{1}{x} W \left[ \frac{\zeta x}{\alpha_{kk}} \exp\left(-x\left(\frac{\bar{\tau}_k^j - (\tau_k^0 + \mu\zeta)}{\mu\alpha_{kk}} + (\gamma_k^{\text{cum},j-1} - \gamma_k^{p,j-1}) - \eta_k\right)\right) \right] \\ \Delta\gamma_k^- &= -\gamma_k^{p,j-1} + \frac{\bar{\tau}_k^j + (\tau_k^0 + \mu\zeta)}{\mu\alpha_{kk}} - \eta_k \\ &\quad - \frac{1}{x} W \left[ \frac{\zeta x}{\alpha_{kk}} \exp\left(x\left(\frac{\bar{\tau}_k^j + (\tau_k^0 + \mu\zeta)}{\mu\alpha_{kk}} - (\gamma_k^{\text{cum},j-1} + \gamma_k^{p,j-1}) - \eta_k\right)\right) \right] \end{aligned} \right. \quad (50) $$

and for the softening phase (i.e., phase III):

$$ \left\{ \begin{aligned} \Delta\gamma_k^+ &= \left(1 - \frac{\tau_k^{\max}}{\mu\alpha_{kk}}\chi\right)^{-1} \left(-\gamma_k^{p,j-1} + \frac{\bar{\tau}_k^j - \tau_k^{\max}(1 + \chi(\gamma_k^{\text{soft}} - \gamma_k^{\text{cum},j-1}))}{\mu\alpha_{kk}} - \eta_k\right) \\ \Delta\gamma_k^- &= \left(1 - \frac{\tau_k^{\max}}{\mu\alpha_{kk}}\chi\right)^{-1} \left(-\gamma_k^{p,j-1} + \frac{\bar{\tau}_k^j + \tau_k^{\max}(1 + \chi(\gamma_k^{\text{soft}} - \gamma_k^{\text{cum},j-1}))}{\mu\alpha_{kk}} - \eta_k\right) \end{aligned} \right. \quad (51) $$

In (50) $\eta_k$ is defined by

$$ \eta_k = \sum_{\substack{l=1 \\ l \ne k}}^{N} \frac{\alpha_{kl}}{\alpha_{kk}} (\Delta\gamma_l^p + \gamma_l^{p,j-1}) \quad (52) $$

and $W$ is the Lambert function defined as the inverse function of $f: z \mapsto z \exp(z)$, which leads to the fundamental property $W(z \exp(z)) = z$.

It should be noted that the solution (50) and (51) to the coupled problem (46) and (48) would necessitate non-linear algorithms as for a given grain $k$ the solution non-linearly depends on all the other grains solution due to coefficient $\eta_k$ defined in (52). In the following, weak interaction between grains is considered to solve explicitly (49) and reach short computation time.
## 2.4. Weak grain interaction

In this section, we present how the proposed solution is linearized under weak interaction assumption. First consider the fully decoupled problem, where grains are independent, which corresponds to $\alpha_{kl} = 0$ if $k \neq l$, leading to $\eta_k = 0$ for any $k$ such as $1 \le k \le N$. In that case (50) and (51) can be directly used as for each grain $k$ the solution does not depend on the solution of the other grains. This solution without interaction is denoted by $\Delta^0\gamma_k^p$, and is obtained by definition when $\eta_k = 0$ in (50) and (51).

Consider now that the interaction terms $\alpha_{kl}$ (with $k \neq l$) are small compared to the terms $\alpha_{kk}$ (i.e., $\alpha_{kl}/\alpha_{kk} \ll 1$, which leads to $\eta_k \ll 1$). This assumption enables us to linearize (50) and (51) using Taylor expansion, which reads for any $k$ such as $1 \le k \le N$:

$$ \Delta\gamma_k^{\pm} \approx \Delta^0\gamma_k^{\pm} + a_k \eta_k \quad (53) $$

where for the hardening phase (i.e., phase I+II):

$$ a_k = \begin{cases} \frac{1}{1 + W \left[ \frac{\zeta x}{\alpha_{kk}} \exp\left(-x\left(\frac{\bar{\tau}_k^j - (\tau_k^0 + \mu\zeta)}{\mu\alpha_{kk}} + (\gamma_k^{\text{cum},j-1} - \gamma_k^{p,j-1})\right)\right) \right]} & \text{if } \Delta^0\gamma_k^+ > 0 \\ \frac{1}{1 + W \left[ \frac{\zeta x}{\alpha_{kk}} \exp\left(x\left(\frac{\bar{\tau}_k^j + (\tau_k^0 + \mu\zeta)}{\mu\alpha_{kk}} - (\gamma_k^{\text{cum},j-1} + \gamma_k^{p,j-1})\right)\right) \right]} & \text{if } \Delta^0\gamma_k^- < 0 \end{cases} \quad (54) $$

and for the softening phase (i.e., phase III):

$$ a_k = - \left(1 - \frac{\tau_k^{\max}}{\mu \alpha_{kk}} \chi\right)^{-1} \quad (55) $$

Therefore, the linearly coupled system reads in matrix form:

$$ \begin{pmatrix} 1 & -a_1 \frac{\alpha_{12}}{\alpha_{11}} & \cdots & -a_1 \frac{\alpha_{1N}}{\alpha_{11}} \\ -a_2 \frac{\alpha_{21}}{\alpha_{22}} & 1 & \cdots & -a_2 \frac{\alpha_{2N}}{\alpha_{22}} \\ \vdots & \vdots & \ddots & \vdots \\ -a_N \frac{\alpha_{N1}}{\alpha_{NN}} & \cdots & -a_N \frac{\alpha_{NN-1}}{\alpha_{NN}} & 1 \end{pmatrix} \begin{pmatrix} \Delta\gamma_1^p \\ \vdots \\ \Delta\gamma_N^p \end{pmatrix} \approx \begin{pmatrix} \Delta^0\gamma_1^p \\ \vdots \\ \Delta^0\gamma_N^p \end{pmatrix} + \begin{pmatrix} a_1 \sum_{l=2}^N \frac{\alpha_{1l}}{\alpha_{11}} \gamma_l^{p,j-1} \\ \vdots \\ a_N \sum_{l=1}^{N-1} \frac{\alpha_{Nl}}{\alpha_{NN}} \gamma_l^{p,j-1} \end{pmatrix} \quad (56) $$
Therefore, using the assumption $\alpha_{kl}/\alpha_{kk} \ll 1$ for $k \neq l$, one obtains:

$$ \left( \begin{array}{c} \Delta \gamma_1^p \\ \vdots \\ \Delta \gamma_N^p \end{array} \right) \approx \left( \begin{array}{cccc} 1 & a_1 \frac{\alpha_{12}}{\alpha_{11}} & \cdots & a_1 \frac{\alpha_{1N}}{\alpha_{11}} \\ a_2 \frac{\alpha_{21}}{\alpha_{22}} & 1 & \cdots & a_2 \frac{\alpha_{2N}}{\alpha_{22}} \\ \vdots & \ddots & \vdots & \vdots \\ a_N \frac{\alpha_{N1}}{\alpha_{NN}} & \cdots & a_N \frac{\alpha_{NN-1}}{\alpha_{NN}} & 1 \end{array} \right) \left[ \begin{array}{c} \Delta^0 \gamma_1^p \\ \vdots \\ \Delta^0 \gamma_N^p \end{array} \right] + \left( \begin{array}{c} a_1 \sum_{l=2}^N \frac{\alpha_{1l}}{\alpha_{11}} \gamma_l^{p,j-1} \\ \vdots \\ a_N \sum_{l=1}^{N-1} \frac{\alpha_{Nl}}{\alpha_{NN}} \gamma_l^{p,j-1} \end{array} \right] \quad (57) $$

In the latter system, all grains are coupled i.e. the effect of all grains on one particular grain is taken into account and conversely. However, coefficients $\alpha_{kl}$ rapidly decrease when the distance between grains $k$ and $l$ increases. Therefore, in practice we assume that $\alpha_{kl} \approx 0$ when the distance between grains $k$ and $l$ is higher than a threshold called interaction radius and denoted by $r$. In addition, usually only one grain $k$ undergoes most of the plastic slip among neighboring grains whose contributions are an order of magnitude lower. That critical grain maximizes the macroscopic resolved shear stress $\bar{\tau}_k^m$ given by (2). If the $k$-th grain is undergoing most of the plastic strain among its neighbors then $a_l \ll a_k$ and $(\Delta^0 \gamma_l^p + \gamma_l^{p,j-1}) \ll (\Delta^0 \gamma_k^p + \gamma_k^{p,j-1})$ for $k \neq l$ for any $l \in \mathcal{V}_k$, where $\mathcal{V}_k$ is the set of neighboring grain indexes (i.e., grains other than $k$ at a distance lower than $r$ from grain $k$). Hence the plastic slip increment in grain $k$ and its neighbors $\forall l \in \mathcal{V}_k$:

$$ \left\{ \begin{array}{l} \Delta \gamma_k^p = \Delta^0 \gamma_k^p + a_k \sum_{l \in \mathcal{V}_k} \frac{\alpha_{kl}}{\alpha_{kk}} (\Delta^0 \gamma_l^p + \gamma_l^{p,j-1}) \\ \\ \Delta \gamma_l^p = \Delta^0 \gamma_l^p + a_l \frac{\alpha_{lk}}{\alpha_{ll}} (\Delta^0 \gamma_k^p + \gamma_k^{p,j-1}) \end{array} \right. \quad (58) $$

Based on (58) the cumulative plastic strain $\gamma_k^{\text{cum},j}$ can be easily computed in all grains.

## 2.5. Integral computation and ellipsoid approximation

The plastic strain increment is determined by using (58) together with (50) and (51). However, $\alpha_{kl}$ coefficients involved in (58) are defined in (30) as integrals over each grain domain $\Omega_k$ of tensors implicitly defined in (19). The explicit determination of $S_{kk}$ involved in $\alpha_{kk}$ is analytical for specific shapes such as ellipsoids, for which $S_{kk}$ reduce to the classical Eshelby tensor [40]. However, the computation of $S_{kl}$ for $k \neq l$ involved in coupling coefficients $\alpha_{kl}$ cannot be reduced to classical formulas. As detailed in (21), $S_{kl}$ (with $k \neq l$) tensors necessitate by definition to integrate $S_l$, which is associated to the solution $\varepsilon_l(\mathbf{x})$ for which only the grain $l$ is subjected to plastic deformation, over the domain $\Omega_k$ lying fully outside $\Omega_l$. Therefore, the solution $\varepsilon_l(\mathbf{x})$ should be explicitly derived outside $\Omega_l$. For specific shapes such as ellipsoids, the coefficients $\alpha_{kl}$ are linked to fourth order symmetric tensor denoted by $H_k^{\text{ext}}(\mathbf{x})$ and called the exterior Eshelby tensor of the grain $k$ computed at any point $\mathbf{x}$ outside the grain $k$ (i.e., $\mathbf{x} \notin \Omega_k$). The resulting strain field $\sigma_k(\mathbf{x})$ outside grain $k$ associated to the plastic strain $\varepsilon_k^p$ (i.e., which takes place in grain
k) reads as

$$ \forall \boldsymbol{x} \notin \Omega_k, \varepsilon_k(\boldsymbol{x}) = (\boldsymbol{C}^{-1} : \boldsymbol{H}_k^{\text{ext}}(\boldsymbol{x}) + \mathbb{I}) : \varepsilon_k^p \quad (59) $$

Hence by definition (19) we obtain

$$ \forall \boldsymbol{x} \notin \Omega_k, \boldsymbol{S}_k(\boldsymbol{x}) = \boldsymbol{C}^{-1} : \boldsymbol{H}_k^{\text{ext}}(\boldsymbol{x}) + \mathbb{I} \quad (60) $$

Definition (21) yields

$$ \boldsymbol{S}_{lk} = \frac{1}{|\Omega_l|} \left[ \boldsymbol{C}^{-1} : \left( \int_{\Omega_l} \boldsymbol{H}_k^{\text{ext}}(\boldsymbol{x}) \mathrm{d}\Omega \right) + \mathbb{I} \right] \quad (61) $$

For ellipsoids an explicit analytical form is given for $\boldsymbol{H}_k^{\text{ext}}(\boldsymbol{x})$ in [45]. Thus, to take advantage of classical analytical formula, grain domains $\Omega_k$ are approximated by ellipsoids denoted by $\tilde{\Omega}_k$ (with $1 \le k \le N$). It should be noted that ellipsoids should be fully embedded in the grains (i.e., $\tilde{\Omega}_k \subset \Omega_k$) to make sure that $\tilde{\Omega}_k$ are disjoint from each other. Otherwise, the computation of $\tilde{\boldsymbol{H}}_k^{\text{ext}}(\boldsymbol{x})$ is singular when $\boldsymbol{x} \in \tilde{\Omega}_l \cap \tilde{\Omega}_k$ (with $k \ne l$), which makes the integral over $\tilde{\Omega}_l$ diverge; hence failing to compute coupling coefficients $\alpha_{kl}$. Therefore the mesostructure is approximated by best fitting ellipsoids under the constrain that ellipsoids should be fully embedded in the grains, which is achieved by using the MATLAB [46] tool box [38] based on convex optimization [47].

## 2.6. Asymptotic plastic strain

In the case of elastic shakedown, the plastic strain stabilizes after a certain number of cycles. Instead of computing a large number of cycles according to the solution derived in section 2.4, an asymptotic development enables to directly determine the stabilized plastic strain in grain $k$ (denoted by $\gamma_k^\infty$) as well as the cumulative plastic strain (denoted by $\gamma_k^{\text{cum},\infty}$). Such an asymptotic formula is derived under the assumption that there is no grain interaction, hence using the fully decoupled solution $\Delta^0\gamma_k^p$ given by (50) with $\eta_k = 0$. This assumption is consistent with elastic shakedown corresponding to very low applied stress $\bar{\sigma}$ and hence favoring isolated grains undergoing plastic slips. It should be noted that the proposed asymptotic formula for the elastic shakedown regime is of course derived for the hardening phase (i.e., phase I+II) as softening would correspond to rapid accumulation of plastic slip.

The asymptotic formula is derived for oscillating applied stress, where $\bar{\tau}_k^0$ denoted the average applied stress and $\Delta\bar{\tau}_k$ is half of the stress amplitude. Thus the stress cycle maximum reads as $\bar{\tau}_k^j = \bar{\tau}_k^0 + \Delta\bar{\tau}_k$ leading to $\Delta^0\gamma_k^p \ge 0$ while the minimum reads as $\bar{\tau}_k^j = \bar{\tau}_k^0 - \Delta\bar{\tau}_k$ leading to $\Delta^0\gamma_k^p \le 0$. Elastic shakedown is reached as soon as pure elastic loading cycles take place, which correspond to $\Delta^0\gamma_k^p = 0$ for both lines in (50), i.e.

$$ \left\{ \begin{array}{l} 0 = -\gamma_k^{p,\infty} + \frac{\bar{\tau}_k^0 + \Delta\bar{\tau}_k - (\tau_k^0 + \mu\zeta)}{\mu\alpha_{kk}} \\ \qquad + \frac{1}{x} W \left[ \frac{\zeta x}{\alpha_{kk}} \exp \left( -x \left( \frac{\bar{\tau}_k^0 + \Delta\bar{\tau}_k - (\tau_k^0 + \mu\zeta)}{\mu\alpha_{kk}} + (\gamma_k^{\text{cum},\infty} - \gamma_k^{p,\infty}) \right) \right) \right] \\[1.5em] 0 = -\gamma_k^{p,\infty} + \frac{\bar{\tau}_k^0 - \Delta\bar{\tau}_k + (\tau_k^0 + \mu\zeta)}{\mu\alpha_{kk}} \\ \qquad - \frac{1}{x} W \left[ \frac{\zeta x}{\alpha_{kk}} \exp \left( x \left( \frac{\bar{\tau}_k^0 - \Delta\bar{\tau}_k + (\tau_k^0 + \mu\zeta)}{\mu\alpha_{kk}} - (\gamma_k^{\text{cum},\infty} + \gamma_k^{p,\infty}) \right) \right) \right] \end{array} \right. \quad (62) $$
Hence:

$$ egin{cases} 
\exp(-x\gamma_k^{\text{cum},\infty}) = \frac{\alpha_{kk}}{\zeta} \left( \frac{-\bar{\tau}_k^0 - \Delta\bar{\tau}_k + \tau_k^0 + \mu\zeta}{\mu\alpha_{kk}} + \gamma_k^{p,\infty} \right) \\ 
\exp(-x\gamma_k^{\text{cum},\infty}) = \frac{\alpha_{kk}}{\zeta} \left( \frac{\bar{\tau}_k^0 - \Delta\bar{\tau}_k + \tau_k^0 + \mu\zeta}{\mu\alpha_{kk}} - \gamma_k^{p,\infty} \right) 
\end{cases} 
\quad (63) $$

Eventually one obtains the stabilized plastic slip $\gamma_k^{p,\infty}$ and cumulative plastic slip $\gamma_k^{\text{cum},\infty}$ as

$$ egin{cases} 
\gamma_k^{p,\infty} = \frac{\bar{\tau}_k^0}{\mu \alpha_{kk}} \\ 
\gamma_k^{\text{cum},\infty} = -\frac{1}{x} \ln \left( \frac{\alpha_{kk}}{\zeta} \left( \frac{-\Delta\bar{\tau}_k + \tau_k^0 + \mu\zeta}{\mu\alpha_{kk}} \right) \right) 
\end{cases} 
\quad (64) $$

There are some restrictions on the loading parameter $\Delta\bar{\tau}_k^0$ for Eq. (64) to apply. In more details, requiring that the argument of the logarithm is positive and that the returned value for $\gamma_k^{\text{cum},\infty}$ is positive implies that

$$ \tau_k^0 \le \Delta\bar{\tau}_k \le \tau_k^0 + \mu\zeta \quad (65) $$

Anticipating on the numerical results presented in Section 3.1, the case where $\Delta\bar{\tau}_k < \tau_k^0$ corresponds to a purely elastic evolution as represented in figure 8(a). The case where $\Delta\bar{\tau}_k > \tau_k^0 + \mu\zeta$ (corresponding to a large loading amplitude) is represented in figure 8(c). In that case, the plastic evolution converges towards a pseudo-steady state cycle as the critical shear stress $\tau_k^C$ approaches $\tau_k^{\text{max}}$ in phase II (see figure 2). The maximum value of the plastic slip $\gamma_{k,\text{max}}^{p,\infty}$ and its minimum value $\gamma_{k,\text{min}}^{p,\infty}$ over that steady state cycle can be approximated by replacing $\tau_k^C$ with $\tau_k^{\text{max}} \approx \tau_k^0 + \mu\zeta$, in which case we obtain from (50) that

$$ \gamma_{k,\text{max}}^{p,\infty} = \frac{\bar{\tau}_k^0 + \Delta\bar{\tau}_k - \tau_k^0 - \mu\zeta}{\mu\alpha_{kk}}, \quad \gamma_{k,\text{min}}^{p,\infty} = \frac{\bar{\tau}_k^0 - \Delta\bar{\tau}_k + \tau_k^0 + \mu\zeta}{\mu\alpha_{kk}} \quad (66) $$

## 2.7. Fatigue lifetime prediction

In the following we consider the grain $k$ with the highest stress and use the asymptotic response to apply relevant fatigue criteria. It should be mentioned that kinematic hardening has not been included in following numerical results. Several cases need to be considered.

**Pure elasticity with no fatigue.** If the resolved shear stress $\sigma_k^{ns} = \sigma : (n_k \otimes s_k)$ in grain $k$ is below the initial critical shear stress $\tau_k^0$ then the grain remains of course purely elastic and there is no fatigue.

**Elastic shakedown and infinite endurance.** It should be noted that a continuous single crystal behavior is considered (see figure 2), which makes the separation between phase I and II somehow arbitrary. In the following $\tau_k^{\text{pla}}$ is a notation to represent the threshold identified in the numerical experiment between elastic and plastic shakedown and is therefore an outcome of the polycrystal response to cyclic loading, not a model parameter.
If $\tau_k^0 \le \sigma_k^{ns} \le \tau_k^{pla}$ (i.e., phase I in figure 2) plastic deformations take place, and the fatigue lifetime may be infinite, which is usually determined using e.g. the Dang-Van criterion [5]. Thus infinite endurance is reached if:

$$ \sigma_k^{eq} \le \sigma^b \qquad (67) $$

where $\sigma^b$ (MPa) is a positive stress threshold and:

$$ \sigma_k^{eq} = \sup (\sigma_k^{ns} + a \sigma_k^{nn}) \qquad (68) $$

where $a$ is a positive dimensionless coefficient, and $\sigma_k^{ns} = \sigma : (n_k \otimes s_k)$ is the resolved shear stress and $\sigma_k^{nn} = \sigma : (n_k \otimes n_k)$ is the normal stress while $\sigma$ is the stress tensor computed according to (9) and $n_k$ and $s_k$ are the normal vector and slip direction of the activated slip mechanism in grain $k$. Condition (67) is related to damage indicators proposed by [48] and can be interpreted as a condition for crack initiation in the plastic slip plane ($n_k, s_k$). It should be noted that the equivalent stress $\sigma_k^{eq}$ is easily computed using the stabilized plastic slip $\gamma^{p,\infty}$ defined in (64). In more detail, the stabilized stress in the critical grain $k$ is given by $\sigma = \bar{\sigma} + C : (S_{kk} - I) : \varepsilon^{p,\infty}$ where $\varepsilon^{p,\infty} = \frac{1}{2}\gamma^{p,\infty}(n_k \otimes s_k + s_k \otimes n_k)$ is the plastic strain in the shakedown state. Using (64) we obtain

$$ \sigma_k^{ns} = \bar{\tau}_k - \bar{\tau}_k^0, \quad \sigma_k^{nn} = \bar{\sigma}_{nn} - \frac{\beta_{kk}}{\alpha_{kk}} \bar{\tau}_k^0 \qquad (69) $$

where $\bar{\sigma}_{nn} = n_k : \bar{\sigma} : n_k$ and $\beta_{kk} = (n_k \otimes n_k) : C : (S_{kk} - I) : (n_k \otimes s_k + s_k \otimes n_k)$. Eq. (69) shows that the equivalent stress $\sigma_k^{eq}$ in (68) differs from its macroscopic counterpart $\bar{\sigma}_k^{ns} + a \bar{\sigma}_k^{nn}$.

**Elastic shakedown and Basquin law.** When (67) is not satisfied, i.e. when $\sigma_k^{eq} > \sigma^b$, fatigue failure is expected to occur after a finite number of cycles denoted by $N_f$. In this case if elastic shakedown is reached, a common approach to estimate $N_f$ using a Basquin law [49]:

$$ N_f = A \times \left( \frac{\sigma_k^{eq}}{\sigma^b} \right)^{-B}, \qquad (70) $$

where $A$ and $B$ are material constants and where $\sigma_k^{eq}$ is computed based on (64).

**Plastic shakedown and Manson-Coffin law.** When $\tau_k^{pla} \le \sigma_k^{ns} \le \tau_k^{max}$ the plateau of the hardening curve (i.e., phase II in figure 2) is reached and a repetitive elasto-plastic loop is reached in grain $k$ if the cumulative plastic strain is such as $\gamma_k^{pla} \le \gamma_k^{cum} \le \gamma_k^{soft}$. Then a plastic shakedown regime is reached and the plastic deformation amplitude of the repetitive elasto-plastic loop denoted by $\Delta\varepsilon_k^{loop}$ is computed from direct cycling using (58). In this case the Manson-Coffin's law is used to compute the number of cycle to reach fatigue failure, which reads:

$$ N_f = A' \times (\Delta\varepsilon_k^{loop})^{B'}, \qquad (71) $$

where $A'$ and $B'$ are material coefficients. In the simplest case where the interactions between grains are neglected, the amplitude $\Delta\varepsilon_k^{loop}$ of the plastic strain loop can be obtained from the asymptotic solution (66) as $\Delta\varepsilon_k^{loop} = (\Delta\bar{\tau}_k - \tau_k^0)/\mu\alpha_{kk}$.
**Softening and number of cycles to failure.** Furthermore, if $\gamma_k^{\text{cum}} \ge \gamma_k^{\text{soft}}$ the softening phase (i.e., phase III in figure 2) is reached, and the number of cycles until the plastic slip reaches $\gamma_k^{\text{max}}$ defined in (40) is considered as the number of cycles to failure $N_f$.

A logic chart that summarizes the proposed procedure to analyze fatigue failure is presented in figure 3.

<figure><img src="image_4.png" /><figcaption>Figure 3: Logic chart of fatigue failure procedure.</figcaption></figure>

# 3. Results

## 3.1. Model validation

The proposed approach leading to the recursive formula (58) to compute the increment of plastic slip in all grains is validated by comparison with CPFEM. This numerical computation is carried out in *FreeFEM++* [42] under a 2D assumption and a single loading cycle to spare computation time. It should be noted that for simpler implementation in *FreeFEM++*, $\zeta$ is set to 0 in the non-linear hardening curve (37) leading to a constant critical shear stress in all grains denoted by $\tau^0 = 100$ MPa. In addition it should be noted that for a given applied stress, hardening would decrease plastic strains. Without hardening the effect of grain interaction is therefore maximized (because of larger plastic strains). Since the model relies on weak grain interaction assumption, the
proposed validation without hardening is in fact more severe as hardening may hide discrepancies between the CPFEM computation and the proposed approach due to more limited grain interaction.

A 500 grains tessellation generated with the free software Neper [50] is considered. A single pure shear loading step is applied in the form $\bar{\sigma} = \bar{\tau} (e_x \otimes e_y + e_y \otimes e_x)$ where $(e_x, e_y)$ denotes the plane directions and $\bar{\tau} = 101$ MPa. It should be noted that any external load could have been used instead of pure shear. Shear-type loading were initially chosen for single grain analysis in order to facilitate the comparison between the loading direction and the activated shear slip-plane as a function of crystal orientation. All grains are FCC and oriented so that their <111> plane belongs to the computation plane $(e_x, e_y)$. Crystal orientation then reduces to a single angle denoted by $\theta (^\circ)$ ranging from 0 to 60° because of crystal hexagonal periodicity in the <111> plane. Hence, there are six slip systems (i.e., [101], [101], [011], [011], [110], [110]). The resolved shear stress $\sigma_k^{ns}$ depends on crystal orientation. Considering the 6 slip planes, the resolved shear stress is maximal for $\theta = 0^\circ$ and minimal for $\theta = 15^\circ$ and the same pattern is repeated as shown in figure 4. Thus for the sake of comparison, crystal orientations are randomly assigned between 0 and $15^\circ$ following a normal distribution of mean $\mu_\theta = 7.5^\circ$ and standard deviation $\sigma_\theta = 2^\circ$. The resulting oriented tessellation is presented in figure 6(a). To avoid edge effects in the FreeFEM++ computation, the polycrystal is embedded in a purely elastic domain so that boundary conditions $\bar{\sigma}.n$ are applied sufficiently far from the grains. The tessellation is meshed using 80,000 uniform triangular elements, and the surrounding elastic domain is meshed using non-regular triangular elements as shown in figure 5. The displacement field is discretized as a piecewise linear function (linear per element). The computation time of the FreeFEM++ computation is around 6 hours on a personal computer. It should be noted that 3-node elements are usually associated with poor performance in comparison to 4-node elements (or higher order elements). However, 3-node elements are easier to use considering the implementation of crystalline plasticity in FreeFEM++. Indeed the strain field is uniform in each element and therefore easily accessible in FreeFEM++ whereas for elements with 4 or more nodes, it would have been necessary to use strains at Gauss points, which are unfortunately not directly accessible in FreeFEM++. Nevertheless, using 3-node elements enabled us to use a high mesh density, which mitigates inaccuracy issues.

Resulting plastic slips are presented with the same color scale in figure 6.c for the CPFEM computation, and in figure 6.d for the present work. A good agreement is observed. In both computations the same 51 grains undergo plastic deformation with the same activated slip system. In addition, to assess more quantitatively the quality of the proposed approach a dissipation distance is introduced for each grain $k$ undergoing plastic deformation:

$$ \text{err}_k = \frac{|\Omega_k| |\gamma_{k,FEM}^p| - |\hat{\Omega}_k| |\gamma_k^p|}{\max_{1 \le k \le 51} [\Omega_k |\gamma_{k,FEM}^p|]} \quad (72) $$

where $|\Omega_k|$ and $|\hat{\Omega}_k|$ are the real and ellipsoid approximated volumes of grain $k$, whose sums over all grains are both normalized to 1. In addition, $\gamma_{k,FEM}^p$ and $\gamma_k^p$ are the average plastic slip in grain $k$ for the CPFEM computation and the proposed approach respectively. It should be noted that the relative dissipation distance (72) is computed with respect to the maximal grain dissipation over all the grains. This choice is relevant as only the grains undergoing the highest plastic deformation
play a role in fatigue lifetime prediction. The dissipation distance (72) is presented as a percentage of the maximal grain dissipation in figure 6(b). The average over the 51 grains undergoing plastic deformation is around 4.5%. This comparison therefore partly validates the quality of the present model. In particular, the relative dissipation distance is small in the grains with the highest plastic grains, which are the most critical with respect to fatigue. A closer inspection at the results in Fig. 6 shows that large relative dissipation distance corresponds either to clusters of neighboring grains with similar plastic strains (such as the group of 4 grains near the middle of the left side), or to isolated grains whose shape is poorly approximated by an ellipsoid (such as the isolated plastifying grain near the upper right corner).

<figure><img src="image_5.png" /><figcaption>Figure 4: Resolved stress $\sigma_k^{ns}$ as a function of crystal orientation $\theta$.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Figure 5: Mesh of the CPFEM computation.</figcaption></figure>

## 3.2. *Strain stress curves and sensitivity analysis*

It should be noted that the purpose of the paper is to demonstrate computational capabilities of the proposed approach regardless of a specific material. Material parameters used in the following
<figure><img src="image_7.png" /><figcaption>Figure 6: a) Oriented tessellation (500 grains), b) Relative dissipation distance with respect to the maximal grain dissipation, c) Average plastic slip per grain for the CPFEM computation, d) Average plastic slip per grain for the present work.</figcaption></figure>

are not identified from experimental measurements but are chosen more or less arbitrarily. In addition, even though there is significant number of model parameters, all are common material parameters for which experimental identification procedures exist. For instance the critical shear stress $\tau_k^0$ and corresponding hardening curves may be identified on polycrystalline structures under tension/compression tests using digital image correlation in a scanning electron microscope.

To characterize typical response of the proposed solution derived in section 2, four different cyclic loading paths are computed for a single isolated grain $k$ to demonstrate that the solution (58) is able to capture the four elasto-plastic regimes in fatigue, namely pure elasticity, elastic shakedown, plastic shakedown and ratcheting (even though model assumptions make the model more appropriate for low applied stress and hence for high cycle fatigue). In this example the shear modulus is set to $\mu = 80769$ MPa, the initial critical shear stress and hardening coefficients involved in (37) read $\tau_k^0 = 100$ MPa, $\zeta = 100/\mu$ and $x = 100$, and the dimensionless coefficient defined in (30) reads $\alpha_{kk} = 0.6$. From the hardening rule (37) the critical shear stress increases with deformation from $\tau_k^0$ to $\tau_k^{\max} \approx \tau_k^0 + \mu\zeta$. The four different cyclic loading are presented in figure 7. Of course when the resolved shear stress remains below $\tau_k^0$ the grain remains fully elastic for all cycles as shown in figure 8(a). When the resolved shear stress reaches its maximum amplitude between $\tau_k^0$ and $\tau_k^{\max}$ an elastic shakedown regime is obtained after a certain number of cycles as shown in figure 7(b). Indeed because of hardening, once the critical shear stress reaches
the maximum applied resolved shear stress $\bar{\tau}_k^j$, subsequent loading cycles remain elastic. As shown in figure 8(c) and (d), cyclic loading whose maximum amplitude is higher than the maximal critical shear stress $\tau_k^{\max}$ results in plastic shakedown regime or even in ratcheting if the loading amplitude increases along with cycles. Indeed, at each cycle the applied resolved shear stress $\bar{\tau}_k^j$ overcomes the critical shear stress, which results in plastic slips.

<figure><img src="image_8.png" /><figcaption>Figure 7: Four different loading paths.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Figure 8: Stress strain curves for a) elastic regime, b) elastic shakedown, c) plastic shakedown, and d) ratcheting.</figcaption></figure>

In addition, a sensitivity analysis is carried out to identify the parameters that have a significant influence on the plastic slip $\gamma_k^{p,j}$ and the cumulative plastic slip $\gamma_k^{cum,j}$ of a given loading cycle $j$.
and grain $k$, and consequently on the lifetime. Results (50) and (51) show that the microstructural parameters likely to influence lifetime are : crystalline orientation, which plays a role through $\bar{\tau}_k^j$, the shape of the ellipse describing the grain, which appears in the expression of the coefficient $\alpha_{kk}$, and grain size through the Hall-Petch relationship (34). In the following only the effect of crystal orientation is demonstrated for the sake of simplicity.

The fatigue lifetime of a given RVE is usually determined by the most critical grains that have the shortest lifetime. The proposed sensitivity analysis though is performed for a fixed single critical grain $k$ for which crystal orientation evolves over the entire range $0^\circ$ to $60^\circ$. It is thus assumed that the most critical grain can have any crystal orientation with the same probability. This corresponds to situations with high stress gradient for which RVEs contain only few grains, and hence the probability density function of crystal orientation of the critical grain is uniform.

The effect of crystal orientation is demonstrated for a fixed critical grain shape corresponding to an ellipse for which the minor to the major axis ratio is set to 0.5, and the angle between the major axis and the horizontal axis $e_x$ is denoted by $\Psi (^\circ)$ and set to $\Psi = 0^\circ$. Hardening parameters involved in (34) are listed in table 1, and parameters involved in fatigue criteria (67), (68), (70) and (71) are listed in table 2. The macroscopic applied stress given in (1) where $f(j)$ oscillates between +1 and -1 with a zero average value, and $\bar{\sigma} = \bar{\tau} (e_x \otimes e_y + e_y \otimes e_x)$. Three different conditions are tested with **condition 1**: $\bar{\tau} = 79.21$ MPa, **condition 2**: $\bar{\tau} = 100.53$ MPa, and **condition 3**: $\bar{\tau} = 121.84$ MPa.

The equivalent stress $\sigma_k^{eq}$ and the fatigue lifetime $N_f$ are presented for **condition 1** as a function of crystal orientation $\theta$ respectively in figures 9 and 10 (where it should be noted that each dot correspond to a different computation according to a given crystal orientation of the critical grain $k$). Loading conditions are such as elastic shakedown is reached for the entire range $0^\circ - 60^\circ$ so that the fatigue lifetime is computed according to the Basquin law (70). According to (68) since the equivalent stress $\sigma_k^{eq}$ not only includes the resolved shear stress $\sigma_k^{ns}$ but also the resolved normal stress $\sigma_k^{nn}$ the evolution presented in figures 9 and 10 is different from the resolved shear stress evolution alone presented in figure 4. Significant variations of the fatigue lifetime with a maximum at $N_f \approx 12 \times 10^8$ cycles are observed around $\theta = 17^\circ$, while being rather uniform with $N_f \approx 10^8$ cycles for the rest of the range. Thus for the tested cyclic load the fatigue lifetime is very often limited to $N_f \approx 10^8$ cycles excepted when the material is strongly textured with a narrow distribution of crystal orientations around $17^\circ$.

The fatigue lifetime $N_f$ is presented for **condition 2** as a function of crystal orientation $\theta$ in figure 11. The applied stress is higher than for condition 1, and is such as elastic shakedown (i.e $N_f$ is computed according to the Basquin law (70)) is reached only for a part of the range of crystal orientations, while plastic shakedown (i.e $N_f$ is computed according to the Manson-Coffin law (71)) is reached for crystal orientations for which the resolved shear stress $\sigma_k^{ns}$ is higher than the threshold $\tau_k^{pla}$ without reaching the softening zone (i.e., phase III in figure 2). Steep variations are observed at the transition between elastic and plastic shakedown in figure 11 depending on the chosen coefficients for each fatigue law. It is clear from figure 11 that for the given cyclic load, the fatigue lifetime can be significantly increased if the material is textured so that crystal orientations around $0^\circ$ and $30^\circ$ are avoided.

The fatigue lifetime $N_f$ is presented for **condition 3** as a function of crystal orientation $\theta$ in figure 12. The applied stress is higher than for condition 2 and therefore assumptions underlying
the model are not necessarily verified. Results may therefore be unreliable. However, even though condition 3 is beyond the scope of the model, it is still presented to complete the analysis and demonstrate model consistency. The applied stress is sufficient to reach ratcheting regime over the entire range 0 – 60° so that the fatigue lifetime is computed according to softening (i.e., phase III in figure 2). Thus the fatigue lifetime is obtained as the number of cycles so that the cumulative plastic shear stress $\gamma_k^{cum}$ reaches $\gamma_k^{max}$. The fatigue lifetime is therefore maximal for crystal orientations corresponding to $\theta = 15°$ and $\theta = 45°$ for which the resolved shear stress is minimal (see. figure 4).

<table><caption>Table 1: Hardening parameters involved in (34).</caption><thead><tr><th>μ<br>(MPa)</th><th>ζ × 10<sup>-4</sup><br>(-)</th><th>ξ<br>(-)</th><th>τ<sup>min</sup><br>(MPa)</th><th>x<br>(-)</th><th>D<sub>k</sub><br>(mm)</th></tr></thead><tbody><tr><td>80769</td><td>4.4</td><td>0.05</td><td>59.2</td><td>1</td><td>0.3536</td></tr></tbody></table>

<table><caption>Table 2: Fatigue criteria parameters.</caption><thead><tr><th>a<br>(-)</th><th>σ<sup>b</sup><br>(MPa)</th><th>A<br>(-)</th><th>B<br>(-)</th><th>A'<br>(-)</th><th>B'<br>(-)</th></tr></thead><tbody><tr><td>0.7</td><td>30</td><td>1.558 × 10<sup>20</sup></td><td>3.6</td><td>0.65</td><td>2.13</td></tr></tbody></table>

<figure><img src="image_10.png" /><figcaption>Figure 9: Equivalent stress $\sigma_k^{eq}$ as a function of crystal orientation $\theta^°$ for condition 1.</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>Figure 10: Fatigue lifetime $N_f$ as a function of crystal orientation $\theta^\circ$ for condition 1.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>Figure 11: Fatigue lifetime $N_f$ as a function of crystal orientation $\theta^\circ$ for condition 2.</figcaption></figure>
<figure><img src="image_13.png" /><figcaption>Figure 12: Fatigue lifetime $N_f$ as a function of crystal orientation $\theta^\circ$ for condition 3.</figcaption></figure>

The proposed approach is also used to derive probabilistic S-N curves when crystal orientation is the only variability of the critical grain. Hence, for each applied stress the fatigue lifetime $N_f$ is computed for a large number of crystal orientations of the critical grain. In figure 13 the resulting S-N curve is presented as a box plot. For each applied stress, the blue box spreads from the 1st to the 3rd quartile of the distribution of fatigue lifetime $N_f$ and the 2nd quartile (i.e., median) is marked in red. A clear change of regime is observed at $\bar{\tau} = 107.632$ MPa for which results are widely dispersed due to the effect of crystal orientation leading to elastic or plastic shakedown. For higher applied stress ratcheting and plastic shakedown take place but model assumptions may limit the reliability of results. On the contrary, for lower applied stress, model assumptions are valid and then very high cycle fatigue is computed accurately. This behavior is qualitatively the same as experimental results reported on steels, e.g., [51] shows similar S-N curves with a rapid decrease corresponding to ratcheting before $10^5$ cycles, no data point between $10^5$ and $10^7$ cycles corresponding to plastic shakedown, and a decrease corresponding to high cycle fatigue after $10^7$ cycles; [52] also shows the same behavior for ratcheting and plastic shakedown.
<figure><img src="image_14.png" /><figcaption>Figure 13: Computed probabilistic S-N curve when crystal orientation is the only variability of the critical grain.</figcaption></figure>

## 4. Conclusion

In this paper an energetic three-dimensional exact solution of crystal plasticity in polycrystals subjected to cyclic loading has been derived for arbitrary shaped grains. The model mainly relies on the assumption that the plastic strain is uniform in each grain and that elastic stiffness can be considered isotropic. In addition, isotropic hardening is considered as well as a linear softening phase when the cumulative plastic slip reaches a certain threshold. The plastic slip in each grain is obtained in the form of a recursive formula coupling all grains and involving integrals over the arbitrary shaped grain domains.

To reach short computation time, weak grain interaction has been assumed, hence decoupling the recursive formula, and the arbitrary shaped grains have been approximated with ellipsoids to facilitate integral computation. The resulting solution enables us to rapidly compute plastic slips and cumulative plastic slips in all grain for very high number of cycles (i.e., up to $10^9$). The proposed approach has been validated in 2D against a classical crystal plasticity finite element model, and satisfying agreement has been observed. In addition, an asymptotic formula has been established for elastic shakedown by neglecting grain interactions, which enables to directly determine the stabilized plastic slip and cumulative plastic slip without computing their evolution through loading cycles.

This model has been combined with classical fatigue criteria to compute the fatigue lifetime associated to any polycrystal and any cyclic load. The effect of crystal orientation of the most loaded grain has been investigated and a probabilistic S-N curve has been derived to demonstrate the model applicability and potential.

**Funding** This work was supported by the Pierre Ailleret EDF fund.

## References

[1] S. Pearson, Initiation of fatigue cracks in commercial aluminium alloys and the subsequent propagation of very short cracks, Engineering Fracture Mechanics 7 (1975) 235–247.
[2] E. Ferrie, J.-Y. Buffiere, W. Ludwig, 3d characterisation of the nucleation of a short fatigue crack at a pore in a cast al alloy using high resolution synchrotron microtomography, International Journal of Fatigue 27 (2005) 1215-1220.

[3] M. Herbig, A. King, P. Reischig, H. Proudhon, E. M. Lauridsen, J. Marrow, J.-Y. Buffière, W. Ludwig, 3-d growth of a short fatigue crack within a polycrystalline microstructure studied using combined diffraction and phase-contrast x-ray tomography, Acta Materialia 59 (2011) 590-601.

[4] H. Maitournam, C. Krebs, A. Galtier, A multiscale fatigue life model for complex cyclic multiaxial loading, International Journal of Fatigue 33 (2011) 232-240.

[5] K. Dang Van, Introduction to fatigue analysis in mechanical design by the multiscale approach, in: High-cycle metal fatigue, Springer, 1999, pp. 57-88.

[6] A. Wöhler, Über die festigkeitsversuche mit eisen und stahl, Ernst & Korn, 1870.

[7] Eurocode 3: design of steel structures: part 1-5: plated structural elements, BSI, London, 2010. Incorporating corrigendum April 2009.

[8] J.-B. Cartiaux, A. Ehrlacher, F. Legoll, A. Libal, J. Reygner, Probabilistic formulation of miner's rule and application to structural fatigue, Probabilistic Engineering Mechanics (2023).

[9] M. Peigney, C. Stolz, An optimal control approach to the analysis of inelastic structures under cyclic loading, Journal of the Mechanics and Physics of Solids 51 (2003) 575-605.

[10] S. Sakout, D. Weisz-Patrault, A. Ehrlacher, Energetic upscaling strategy for grain growth. I: Fast mesoscopic model based on dissipation, Acta Materialia 196 (2020) 261-279.

[11] D. Weisz-Patrault, S. Sakout, A. Ehrlacher, Energetic upscaling strategy for grain growth. II: Probabilistic macroscopic model identified by Bayesian techniques, Acta Materialia (2021).

[12] K. Dang Van, Sur la résistance à la fatigue des métaux, Sciences et Techniques de l'Armement 47 (1973) 641.

[13] K. Dang Van, Y. Papadopoulos, Multiaxial fatigue failure criterion: a new approach, Fatigue'87. 2 (1987) 997-1008.

[14] K. Dang Van, Macro-micro approach in high-cycle multiaxial fatigue, in: Advances in multiaxial fatigue, ASTM International, 1993.

[15] I. Papadopoulos, A high-cycle fatigue criterion applied in biaxial and triaxial out-of-phase stress conditions, Fatigue & Fracture of Engineering Materials & Structures 18 (1995) 79-91.

[16] K. D. Van, I. V. Paradopoulos, High-cycle metal fatigue: from theory to applications, volume 392, Springer, 2014.

[17] E. Charkaluk, A. Constantinescu, H. Maitournam, K. Dang Van, Revisiting the Dang Van criterion, Procedia Engineering 1 (2009) 143-146.

[18] V. Monchiet, E. Charkaluk, D. Kondo, Plasticity-damage based micromechanical modelling in high cycle fatigue, Comptes Rendus Mécanique 334 (2006) 129-136.

[19] J. R. Rice, D. M. Tracey, On the ductile enlargement of voids in triaxial stress fields, Journal of the Mechanics and Physics of Solids 17 (1969) 201-217.

[20] G. Bertolino, A. Constantinescu, M. Ferjani, P. Treiber, A multiscale approach of fatigue and shakedown for notched structures, Theoretical and Applied Fracture Mechanics 48 (2007) 140-151.

[21] C. Robert, N. Saintier, T. Palin-Luc, F. Morel, Micro-mechanical modelling of high cycle fatigue behaviour of metals under multiaxial loads, Mechanics of Materials 55 (2012) 112-129.

[22] R. Guerchais, G. Scalet, A. Constantinescu, F. Auricchio, Micromechanical modeling for the probabilistic failure prediction of stents in high-cycle fatigue, International Journal of Fatigue 87 (2016) 405-417.

[23] R. Guerchais, F. Morel, N. Saintier, Microstructure-dependent predictions of the effect of defect size and shape on the high-cycle fatigue strength, Open repository Art et Métiers (2016).

[24] M. Peigney, C. Stolz, Optimal control approach of elastoviscoplastic structures under cyclic loading, Comptes Rendus de l'Académie des Sciences-Series IIB-Mechanics 329 (2001) 643-648.

[25] H. Maitournam, B. Pommier, J.-J. Thomas, Détermination de la réponse asymptotique d'une structure anélastique sous chargement thermomécanique cyclique, Comptes Rendus Mécanique 330 (2002) 703-708.

[26] K. V. Spiliopoulos, K. D. Panagiotou, A direct method to predict cyclic steady states of elastoplastic structures, Computer Methods in Applied Mechanics and Engineering 223 (2012) 186-198.

[27] T. Nguyen-Tajan, B. Pommier, H. Maitournam, M. Houari, L. Verger, Z. Du, M. Snyman, Determination of the
stabilized response of a structure undergoing cyclic thermal-mechanical loads by a direct cyclic method, in: 16th ABAQUS Users Conference, Munich.

[28] B. Fedelich, A. Ehrlacher, An analysis of stability of equilibrium and of quasi-static transformations on the basis of the dissipation function, European journal of mechanics. A. Solids 16 (1997) 833–855.

[29] A. Mielke, Energetic formulation of multiplicative elasto-plasticity using dissipation distances, Continuum Mechanics and Thermodynamics 15 (2003) 351–382.

[30] C. Miehe, Strain-driven homogenization of inelastic microstructures and composites based on an incremental variational formulation, International Journal for numerical methods in engineering 55 (2002) 1285–1322.

[31] C. Miehe, M. Lambrecht, E. Gürses, Analysis of material instabilities in inelastic solids by incremental energy minimization and relaxation methods: evolving deformation microstructures in finite plasticity, Journal of the Mechanics and Physics of Solids 52 (2004) 2725–2769.

[32] A time-integration scheme for thermomechanical evolutions of shape-memory alloys, Comptes Rendus Mécanique 334 (2006) 266–271.

[33] M. Peigney, J.-P. Seguin, An incremental variational approach to coupled thermo-mechanical problems in anelastic solids. application to shape-memory alloys, International Journal of Solids and Structures 50 (2013) 4043–4054.

[34] G. Scalet, M. Peigney, A robust and efficient radial return algorithm based on incremental energy minimization for the 3d souza-auricchio model for shape memory alloys, European Journal of Mechanics-A/Solids 61 (2017) 364–382.

[35] J. Bluthé, D. Weisz-Patrault, A. Ehrlacher, Energetic approach for a sliding inclusion accounting for plastic dissipation at the interface, application to phase nucleation, International Journal of Solids and Structures 121 (2017) 163–173.

[36] M. Peigney, G. Scalet, F. Auricchio, A time integration algorithm for a 3d constitutive model for smas including permanent inelasticity and degradation effects, International Journal for Numerical Methods in Engineering 115 (2018) 1053–1082.

[37] D. Weisz-Patrault, M. Gantier, A. Ehrlacher, Mixed analytic/energetic approach for a sliding orthotropic hollow cylinder. Application to coil sagging, International Journal of Solids and Structures 165 (2019) 75–92.

[38] A. Mutapcic, Maximum volume inscribed ellipsoid in a polyhedron (2006).

[39] J. D. Eshelby, The determination of the elastic field of an ellipsoidal inclusion, and related problems, in: Proceedings of the Royal Society of London A: Mathematical, Physical and Engineering Sciences, volume 241, The Royal Society, pp. 376–396.

[40] T. Mura, Micromechanics of defects in solids, Springer Science & Business Media, 1987.

[41] F. Hofmann, G. Bertolino, A. Constantinescu, M. Ferjani, Numerical exploration of the Dang Van high cycle fatigue criterion: application to gradient effects, Journal of Mechanics of Materials and Structures 4 (2009) 293–308.

[42] F. Hecht, New development in freefem++, J. Numer. Math. 20 (2012) 251–265.

[43] I. Papadopoulos, A high-cycle fatigue criterion applied in biaxial and triaxial out-of-phase stress conditions, Fatigue & Fracture of Engineering Materials & Structures 18 (1995) 79–91.

[44] L. Tabourot, Vers une vision unifiée de la plasticité cristalline. Mécanique, Ph.D. thesis, Université de Savoie, 2001.

[45] X. Jin, Z. Wang, Q. Zhou, L. M. Keer, Q. Wang, On the solution of an elliptical inhomogeneity in plane elasticity by the equivalent inclusion method, Journal of Elasticity 114 (2014) 1–18.

[46] T. M. Inc., Matlab version: 9.13.0 (r2022b), 2022.

[47] S. Boyd, L. Vandenberghe Schmid, Convex Optimization, Cambridge University Press, 2004.

[48] Z. Mróz, A. Seweryn, Damage description with related crack initiation and propagation conditions, Journal de Physique IV Proceedings 06 (1996) 529–538.

[49] O. Basquin, The exponential law of endurance tests, in: Proc Am Soc Test Mater, volume 10, pp. 625–630.

[50] R. Quey, P. Dawson, F. Barbe, Large-scale 3d random polycrystals for the finite element method: Generation, meshing and remeshing, Computer Methods in Applied Mechanics and Engineering 200 (2019) 1729–1745.

[51] B. M. Schönbauer, M. Fitzka, U. Karr, H. Mayer, Variable amplitude very high cycle fatigue of 17-4ph steel with a stepwise sn curve, International Journal of Fatigue 142 (2021) 105963.
[52] A. Mohabeddine, J. Correia, P. A. Montenegro, A. De Jesus, J. M. Castro, F. Berto, Probabilistic sn curves for cfrp retrofitted steel details, International Journal of Fatigue 148 (2021) 106205.