# Comptes Rendus
## Mécanique

Bastien Agard, Landry Giraud, Françoise Fauvin, Jean-Christophe Roux, Pierre Monnet and Eric Feulvarch

Fast computation of critical planes for fatigue life analysis of metals

Volume 350 (2022), p. 495-506

Published online: 3 November 2022

https://doi.org/10.5802/crmeca.139

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

This article is licensed under the
CREATIVE COMMONS ATTRIBUTION 4.0 INTERNATIONAL LICENSE.
http://creativecommons.org/licenses/by/4.0/

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

Les Comptes Rendus. Mécanique sont membres du
Centre Mersenne pour l'édition scientifique ouverte
www.centre-mersenne.org
e-ISSN : 1873-7234
Short paper / *Note*

# Fast computation of critical planes for fatigue life analysis of metals

Bastien Agard<sup>a</sup>, Landry Giraud<sup>b</sup>, Françoise Fauvin<sup>c</sup>, Jean-Christophe Roux<sup>c</sup>, Pierre Monnet<sup>d</sup> and Eric Feulvarch<sup>*</sup>,<sup>c</sup>

<sup>a</sup> ESI Group, Le Récamier, 70 Rue Robert, 69458 Lyon Cedex 06, France

<sup>b</sup> TRA-C industrie, ZAC les Olmes, 69490 Vindry-sur-Turdine, France

<sup>c</sup> Univ Lyon, Ecole Centrale de Lyon, LTDS UMR 5513 CNRS, 58 rue Jean Parot, 42023 Saint-Etienne cedex 02, France

<sup>d</sup> Haulotte Group, Rue Emile Zola, 42420 Lorette, France

*E-mails:* bastien.acm@gmail.com (B. Agard), l.giraud@tra-c.com (L. Giraud), francoise.fauvin@enise.fr (F. Fauvin), jean-christophe.roux@enise.fr (J.-C. Roux), pmonnet@haulotte.com (P. Monnet), eric.feulvarch@enise.fr (E. Feulvarch)

**Abstract.** The main objective of this paper is to propose an efficient computational strategy for critical plane type approaches dedicated to the fatigue life analysis of metals. The formulation developed consists in determining the critical planes in the space of the strain or stress tensors describing the mechanical history. To the authors' knowledge, such an approach has never been reported before in the literature. Thus, the new numerical method developed avoids scanning all the possible planes as standard approaches. As shown in the examples proposed, this particularity allows decreasing the CPU time: from 20 days by scanning planes to 10 min for the FEM analysis of a welded structure with the criterion of Fatemi and Socie.

**Keywords.** Shear amplitude, Critical planes, Fatigue analysis, Tresca norm, Finite element.

Manuscript received 3 July 2022, revised 13 September 2022, accepted 4 October 2022.

## 1. Introduction

For the multiaxial fatigue life analysis of solid metals, the critical plane type approaches are widely used. One can cite the Fatemi-Socie criterion for low-cycle fatigue [1]. It is based on seeking the greatest value of the shear strain amplitude by scanning all the possible planes at a given material point. In practice, 180° planes need to be tested by considering a 1° angle increment. The ones giving the maximum value of the shear amplitude are retained. In the case of high-cycle fatigue, Bernasconi and Papadopoulos [2] proposed to compute the shear amplitude by means of the radius of the minimal enclosing circle containing all the shear history on each plane. Weber et al. [3] proposed to optimize the angle increment for reducing the

<footnote>* Corresponding author.</footnote>
calculation duration. Unfortunately, experience shows that such an approach cannot be used for analyzing all the integration points of a FEM model because it is too much time consuming.

An alternative approach consists in seeking for an equivalent shear amplitude in the space of the stress or strain tensors as proposed by Dang Van et al. [4]. It consists in computing the shear amplitude by means of the radius of the minimal enclosing hypersphere containing all the tensors describing the cycle loading. Fauvin et al. [5] recently implemented an efficient computational technique for this purpose and proved that this equivalent shear amplitude estimates the one obtained by scanning all the possible planes. Even if this approach avoids scanning all the planes, it is unable to give the critical planes needed for the criterion of Fatemi-Socie applied in this paper.

Rather than seeking for the minimal enclosing hypersphere, Lemaitre and Chaboche [6] proposed to define the shear amplitude by means of the von Mises norm between the two most distant tensors. In the case of the strain history, we propose here to estimate the shear amplitude by means of the longest distance between the tensors in the sense of Tresca. The interest lies in the fact that the critical planes can be deduced from the Mohr's circles without scanning all the possible planes. To the authors' knowledge, such an approach has never been reported before in the literature. For computational efficiency, an algorithm in optimal time $O(nk)$ with $k \ll n$ is also proposed in this paper for searching at the two most distant tensors.

The paper is organized as follows. Section 2 is devoted to the general formulation based on the Tresca norm in the case of strains. The associated numerical algorithm is also presented. The examples detailed in Section 3 show the potential and the efficiency of the method proposed in terms of results and computation times. Even if the mathematical formalism is different from the one of Bernasconi and Papadopoulos [2], comparisons show that similar results are obtained with a highly reduced computation time in the case of multiaxial loadings.

## 2. An efficient computation of critical planes

### 2.1. Theory

The critical planes are the planes on which the shear strain amplitude denoted $\Delta\gamma/2$ is maximal. In order to seek them, a first approach consists in considering a plane defined with the normal vector $\mathbf{n}$ as shown in Figure 1. The shear strain at the instant $i$ of the loading cycle is denoted $\gamma_n^{(i)}$. It can be calculated from the deviatoric part of the strain tensor $\boldsymbol{\epsilon}^{(i)}$ as follows:

$$ \gamma_n^{(i)} = \boldsymbol{\epsilon}_D^{(i)} \cdot \mathbf{n} - ((\boldsymbol{\epsilon}_D^{(i)} \cdot \mathbf{n}) \cdot \mathbf{n})\mathbf{n}. \quad (1) $$

The shear strain amplitude is then evaluated by means of the diameter $\Delta\gamma_n$ of the set of points belonging to the plane. In this paper, it is based on the euclidian distance between the two most farthest tips of the vectors $\gamma_n^{(i)}$:

$$ \Delta\gamma_n = \max_{(i,j)} \|\gamma_n^{(i)} - \gamma_n^{(j)}\|_2. \quad (2) $$

The critical planes are maximizing this diameter, and more precisely they satisfy the following double maximization:

$$ \Delta\gamma = \max_n \Delta\gamma_n = \max_n \left( \max_{(i,j)} \|\Delta\gamma_n^{(i,j)}\|_2 \right), \quad (3) $$

where $\Delta\gamma_n^{(i,j)} = \gamma_n^{(i)} - \gamma_n^{(j)}$.

In order to determine them, one must exhaustively scan all planes around each point of interest. By choosing an increment of 1°, 180° planes must be considered. This double maximization is therefore very expensive in terms of storage area and computational time. We propose in this
<figure><img src="image_3.png" /><figcaption>Figure 1. Diameter Δγn on the plane with normal n.</figcaption></figure>

work to use a permutation of this double maximization in order to solve the problem of finding the critical planes in the space of the deviatoric tensors ε_D^{(k)}. Thus, we propose to write (3) as follows:

$$ \Delta\gamma = \max_{(i,j)} \left( \max_n \| \Delta\gamma_n^{(ij)} \|_2 \right). \quad (4) $$

The shear strain difference vector between the shear states (i) and (j) can be determined using the deviatoric part of the deformation tensor according to the formula:

$$ \begin{gathered} \Delta\gamma_n^{(ij)} = \Delta\epsilon_D^{(ij)} \cdot n - ((\Delta\epsilon_D^{(ij)} \cdot n) \cdot n)n, \quad (5) \\ \text{where } \Delta\epsilon_D^{(ij)} = \epsilon_D^{(i)} - \epsilon_D^{(j)}. \end{gathered} $$

As shown in Figure 2, the Mohr's circles associated with the tensor Δε_D^{(ij)} shows that max_n ||Δγ_n^{(ij)}|| is equal to the radius of the largest Mohr's circle. Thus, we can write:

$$ \max_n \| \Delta\gamma_n^{(ij)} \|_2 = \frac{1}{2} |\Delta\epsilon_1^{(ij)} - \Delta\epsilon_{III}^{(ij)}| = \frac{1}{2} \| \Delta\epsilon^{(ij)} \|_{\text{Tresca}}, \quad (6) $$

where Δε_1^{(ij)}, Δε_II^{(ij)} and Δε_III^{(ij)} are the three eigenvalues of the tensor Δε^{(ij)} such as Δε_1^{(ij)} > Δε_II^{(ij)} > Δε_III^{(ij)}.

The maximum shear strain amplitude Δγ/2 can then be expressed by means of the Tresca norm as follows:

$$ \frac{\Delta\gamma}{2} = \frac{1}{4} \max_{(i,j)} \| \Delta\epsilon^{(ij)} \|_{\text{Tresca}}. \quad (7) $$

In contrast to the standard scanning approaches, one can note that the maximum shear strain amplitude is estimated in the deviatoric space by means of the Tresca norm without paying special attention to the orientation of planes. Because of the angle of π/4 on the Mohr's circles in Figure 2, the two critical planes associated with Δγ/2 are obtained very quickly. Indeed, they are oriented by the two bisectors of the eigenvectors I and III associated with Δε_1^{(ij)} and Δε_III^{(ij)} of the solution of (7). If two eigenvalues are equal, it is necessary to consider four bisectors whereas the shear strain amplitude is zero in the case of three identical eigenvalues.

From the simulation point of view, the use of the Tresca norm has the advantage of avoiding to scan all the possible planes. The main difficulty lies in solving the double maximization given by (7).
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Mohr's circles of the strain amplitude tensor $\Delta\varepsilon$ in 3D.</figcaption></figure>

## 2.2. Computational algorithm

By considering the strain tensors as points in dimension five in the deviatoric space, the double maximization given by (7) comes down to determining the diameter of a set of points by means of the Tresca norm. The trivial and “greedy” algorithm which compares the distances between all the pairs of the $n$ points of the set has a quadratic time complexity $O(n^2)$ in the worst case. In low dimensions (2 or 3), the calculation of the diameter of a set of points is based on optimal algorithmic methods (see e.g. [7–9]). Some of these algorithms could be extended to higher dimensions, but the generalization is quite difficult. However, Preparata and Shamos [10] showed that it is possible to determine the diameter of a set of $n$ points in dimension $m$ in optimal time $O(n \log n)$. In this paper, we propose to apply a computational method working for points in any finite dimension $m$. The algorithmic complexity is $O(nk)$ with $k \ll n$. Thus, the method is very fast in practice.

### 2.2.1. Geometrical and mathematical properties

Let us consider a set $\mathcal{E}$ composed of $n$ points in a normed vector space $\mathcal{P}$. The distance $d$ is then its induced metric $d(P, Q) = \|Q - P\|$.

One can define, in the metric space $\mathcal{P}$:

**Definition 1.** The diameter of $\mathcal{P}$ is the maximum distance between any two points of $\mathcal{P}$:

$$ \operatorname{diam}(\mathcal{P}) = \max_{\mathbf{P}, \mathbf{Q} \in \mathcal{P}} d(\mathbf{P}, \mathbf{Q}). \quad \square $$

**Property 1.** The endpoints of a diameter belong to the convex hull<fnref n="1" /> of set $\mathcal{P}$.

Moreover, if $[\mathbf{P}, \mathbf{Q}]$ is a diameter of $\mathcal{P}$, then $\mathbf{P} \in \{\text{points } S \in \mathcal{P} \text{ such that } d(S, \mathbf{Q}) \text{ is maximum}\}$ and $\mathbf{Q} \in \{\text{points } S \in \mathcal{P} \text{ such that } d(S, \mathbf{P}) \text{ is maximum}\}$. $\square$

If $\mathcal{B}_{[\mathbf{P}, \mathbf{Q}]}$ is the ball with diameter $[\mathbf{P}, \mathbf{Q}]$ we deduce from Property 1 the following result:

**Property 2.** If two points $(\mathbf{P}, \mathbf{Q}) \in \mathcal{P}$ are such that $\mathcal{P} - \{S \in \mathcal{P}_{[\mathbf{P}, \mathbf{Q}]}\} = \emptyset$, then $[\mathbf{P}, \mathbf{Q}]$ is a diameter of $\mathcal{P}$. $\square$

These results suggest the main idea of the algorithm: compute the diameter of a set of points by considering and finding the points belonging to the convex hull, i.e., finding pairs of points

<footnote n="1">¹The convex hull of a set of points $\mathcal{P}$ is the smallest convex set that contains all the points of $\mathcal{P}$.</footnote>
which are the farthest from each other. Note that in our case, the norm will be chosen as the Tresca norm in the deviatoric space [5].

### 2.2.2. Computing the diameter of a set of points—algorithm

In the sequel, $\boldsymbol{\epsilon}^{(i)}$ denotes the $i$th tensor considered as a point of set $\mathcal{E}$; $d_T\boldsymbol{\epsilon}^{(ij)}$ is defined as the distance (based on the Tresca norm) between two points (tensors) $\boldsymbol{\epsilon}^{(i)}$ and $\boldsymbol{\epsilon}^{(j)}$; $d_T\boldsymbol{\epsilon}^{(ij)} = d(\boldsymbol{\epsilon}^{(i)}, \boldsymbol{\epsilon}^{(j)}) = \|\Delta\boldsymbol{\epsilon}^{(ij)}\|_{\text{Tresca}}; \mathcal{B}^{\circ}_{[\text{diameter}]} = \mathcal{B}^{\circ}[\boldsymbol{\epsilon}^{(i)}, \boldsymbol{\epsilon}^{(j)}]$ is the open ball with diameter $[\boldsymbol{\epsilon}^{(i)}, \boldsymbol{\epsilon}^{(j)}]$. The computational method for the diameter computation is described by Algorithm 1.

<caption>**Algorithm 1 – Diameter computation**</caption>

**Data:** Set of points (tensors) $\mathcal{E}$
**Result:** Diameter = (list of) $[\boldsymbol{\epsilon}^{(i)}, \boldsymbol{\epsilon}^{(j)}]$:
// (list of) two endpoints: diameter(s) of set $\mathcal{E}$.
**Initialization:** $\boldsymbol{\epsilon}^{(i)}$ ← (random) point of $\mathcal{E}$;
$d_T\boldsymbol{\epsilon}^{(ij)} ← 0$; diameter ← $\emptyset$; stop ← false
**repeat**
> **Find** ($\boldsymbol{\epsilon}^{(k)} \in \mathcal{E}_{\text{current}}$) **such as** ($d_T\boldsymbol{\epsilon}^{(ik)}$ is maximum)
> &nbsp;&nbsp;&nbsp;&nbsp;**if** $d_T\boldsymbol{\epsilon}^{(ik)} > d_T\boldsymbol{\epsilon}^{(ij)}$ **then**
> &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;**diameter** ← $[\boldsymbol{\epsilon}^{(i)}, \boldsymbol{\epsilon}^{(k)}]$ // two endpoints
> &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$d_T\boldsymbol{\epsilon}^{(ij)} ← d_T\boldsymbol{\epsilon}^{(ik)}$
> &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$\mathcal{E}_{\text{current}} ← \mathcal{E}_{\text{current}} - \{\boldsymbol{\epsilon} \in \mathcal{B}^{\circ}_{[\text{diameter}]}\}$ // remove from $\mathcal{E}_{\text{current}}$
> &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;// all points $\boldsymbol{\epsilon}$ in the open ball with diameter $[\boldsymbol{\epsilon}^{(i)}, \boldsymbol{\epsilon}^{(k)}]$
> &nbsp;&nbsp;&nbsp;&nbsp;**if** $\mathcal{E}_{\text{current}} - \{\boldsymbol{\epsilon}^{(i)}, \boldsymbol{\epsilon}^{(k)}\} \neq \emptyset$ **then**
> &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$$ \boldsymbol{\epsilon}^{(m)} \leftarrow \frac{1}{2}(\boldsymbol{\epsilon}^{(i)} + \boldsymbol{\epsilon}^{(k)}) $$
> &nbsp;&nbsp;&nbsp;&nbsp;**Find** ($\boldsymbol{\epsilon}^{(i)} \in \mathcal{E}_{\text{current}}$) **such as** ($d_T\boldsymbol{\epsilon}^{(im)}$ is maximum)
> &nbsp;&nbsp;&nbsp;&nbsp;// another iteration with a new point $\boldsymbol{\epsilon}^{(i)}$ on the convex hull
> &nbsp;&nbsp;&nbsp;&nbsp;**end**
> **else**
> &nbsp;&nbsp;&nbsp;&nbsp;stop ← true
> **end**
> **until** $(\mathcal{E}_{\text{current}} - \{\boldsymbol{\epsilon}^{(i)}, \boldsymbol{\epsilon}^{(k)}\} = \emptyset)$ or (stop = true)
> $\mathcal{E}_{\text{out}} ← \mathcal{E} - \{\boldsymbol{\epsilon} \in \mathcal{B}^{\circ}_{[\text{diameter}]}\}$
> **if** $\mathcal{E}_{\text{out}} - [\text{diameter}] = \emptyset$ **then**
> Diameter ← diameter // unique diameter found
> **else**
> &nbsp;&nbsp;&nbsp;&nbsp;**foreach** point $\boldsymbol{\epsilon}^{(k)} \in \mathcal{E}_{\text{out}}$ **do**
> &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;**foreach** point $\boldsymbol{\epsilon}^{(l)} \in \mathcal{E}$ **do**
> &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;**if** $d_T\boldsymbol{\epsilon}^{(kl)} > d_T\boldsymbol{\epsilon}^{(ij)}$ **then**
> &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$$ d_T\boldsymbol{\epsilon}^{(ij)} ← d_T\boldsymbol{\epsilon}^{(kl)} $$
> &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;list_of_diameters ← $[\boldsymbol{\epsilon}^{(k)}, \boldsymbol{\epsilon}^{(l)}]$
> &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;**else**
> &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;**if** $d_T\boldsymbol{\epsilon}^{(kl)} = d_T\boldsymbol{\epsilon}^{(ij)}$ **then**
> &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$$ list_of_diameters ← list_of_diameters \oplus [\boldsymbol{\epsilon}^{(k)}, \boldsymbol{\epsilon}^{(l)}] $$
> &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;**end**
> &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;**end**
> &nbsp;&nbsp;&nbsp;&nbsp;**end**
> **end**
> Diameter ← list_of_diameters
> **end**
<figure><img src="image_5.png" /><figcaption>**Figure 3.** Evolution of the principal directions during one loading cycle.</figcaption></figure>

<table><caption>**Table 1.** Components of the normal vectors defining the critical planes</caption><thead><tr><th></th><th></th><th>n<sub>1</sub></th><th>n<sub>2</sub></th><th>n<sub>3</sub></th></tr></thead><tbody><tr><td>Bernasconi and Papadopoulos [2]</td><td><b>P</b><sub>1</sub></td><td>-0.6645</td><td>-0.2418</td><td>-0.7071</td></tr><tr><td>New approach</td><td><b>P</b><sub>2</sub></td><td>-0.6645</td><td>-0.2418</td><td>0.7071</td></tr><tr><td></td><td><b>P</b><sub>1</sub></td><td>-0.6662</td><td>-0.2371</td><td>-0.7071</td></tr><tr><td></td><td><b>P</b><sub>2</sub></td><td>-0.6662</td><td>-0.2371</td><td>0.7071</td></tr></tbody></table>

### 3. Applications

#### 3.1. Analysis of a nonproportional multiaxial strain state

We propose here to apply the algorithm detailed above in the case of cyclic loading with varying principal strain directions. The first objective is to compare the critical planes and the corresponding maximal shear strain amplitude $\Delta\gamma/2$ given by (7) and by scanning all the possible planes as proposed by Bernasconi and Papadopoulos [2]. The second objective is to compare the computational time to show the interest in the technique developed.

Let us consider a material point subjected to cyclic loading inducing the following strain tensor:

$$ \epsilon(t) = \begin{bmatrix} \epsilon_{11}(t) & \epsilon_{12}(t) & 0 \\ \epsilon_{12}(t) & \epsilon_{22}(t) & 0 \\ 0 & 0 & \epsilon_{33}(t) \end{bmatrix}, \qquad (8) $$

where $\epsilon_{11}(t) = 0.005 \cos(3\pi t)$, $\epsilon_{22}(t) = -0.001 \cos(2\pi t)$, $\epsilon_{33}(t) = 0.0075 \cos(4\pi t)$ and $\epsilon_{12}(t) = -0.002 \cos(5\pi t)$. During the loading cycle, the principal directions vary strongly in the plane perpendicular to the direction numbered 3 as shown in Figure 3. Therefore, the multiaxial strain state is highly nonproportional.

Table 1 summarizes the critical planes. It is to be noted that the approach of Bernasconi and Papadopoulos [2] gives two perpendicular critical planes as predicted by the formulation developed. Moreover, the two planes obtained by means of each of these methods are very close because the maximum difference is lower than 0.3° in terms of directions.

Table 2 gives the shear strain amplitudes. The value computed by seeking for the critical planes as proposed by Bernasconi and Papadopoulos [2] is equal to the one obtained by means of the new approach proposed.
<figure><img src="image_6.png" /><figcaption>**Figure 4.** Comparison of the computational time.</figcaption></figure>

<table><caption>**Table 2.** Shear strain amplitude</caption><thead><tr><th></th><th>Δγ / 2 (%)</th></tr></thead><tbody><tr><td>Bernasconi and Papadopoulos [2]</td><td>1.246</td></tr><tr><td>New approach</td><td>1.246</td></tr></tbody></table>

Figure 4 shows the comparison of the computational times with a constant angle increment of 1° for the approach of Bernasconi and Papadopoulos [2]. All the numerical experiments are based on a regular distribution of the temporal points along the whole loading cycle. One can note that the computational times corresponding to the new approach do not evolve significantly and remain notably below those of the classical approach of Bernasconi and Papadopoulos [2].

## 3.2. Analysis of strain gauge measurements

We propose here to apply the computational method developed for the analysis of an industrial structure subjected to a loading cycle and equipped with a rosette composed of three gauges with an angle of 45°. The acquisition frequency is equal to 100 Hz for a cycle duration of 725 s.

Because of the plane stress configuration on the surface of the structure perpendicular to the direction 3, the strain tensor is of the following form:

$$ \epsilon(t) = \begin{bmatrix} \epsilon_{11}(t) & \epsilon_{12}(t) & 0 \\ \epsilon_{12}(t) & \epsilon_{22}(t) & 0 \\ 0 & 0 & \epsilon_{33}(t) \end{bmatrix}. \qquad (9) $$

As for the previous example, the principal directions vary significantly as shown in Figure 5. This means that the multiaxial strain state is highly nonproportional.

The computed critical planes are summarized in Table 3. Once again, the approach of Bernasconi and Papadopoulos [2] gives two critical planes which are perpendicular (see Figure 6). The planes obtained by means of the methods are again very close because the maximum difference is here lower than 0.05° in terms of directions.

Table 4 gives the shear strain amplitudes. One can see that the value computed by seeking the critical planes [2] is once again equal to the one obtained by means of the new approach proposed.
<figure><img src="image_7.png" /><figcaption>**Figure 5.** Evolution of the principal directions.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**Figure 6.** Critical planes.</figcaption></figure>

<table><caption>**Table 3.** Components of the normal vectors defining the critical planes</caption><thead><tr><th></th><th></th><th>n<sub>1</sub></th><th>n<sub>2</sub></th><th>n<sub>3</sub></th></tr></thead><tbody><tr><td>Bernasconi and Papadopoulos [2]</td><td><b>P</b><sub>1</sub></td><td>0.5997</td><td>-0.3747</td><td>0.7071</td></tr><tr><td rowspan="2">New approach</td><td><b>P</b><sub>2</sub></td><td>0.5997</td><td>-0.3747</td><td>-0.7071</td></tr><tr><td><b>P</b><sub>1</sub></td><td>0.6</td><td>-0.3742</td><td>0.7071</td></tr><tr><td></td><td><b>P</b><sub>2</sub></td><td>0.6</td><td>-0.3742</td><td>-0.7071</td></tr></tbody></table>
<figure><img src="image_9.png" /><figcaption>**Figure 7.** Geometry (left) and failure at the edge of the weld bead (right).</figcaption></figure>

<table><caption>**Table 4.** Shear strain amplitude</caption><thead><tr><th></th><th>Δγ / 2 (%)</th></tr></thead><tbody><tr><td>Bernasconi and Papadopoulos [2]</td><td>0.1134</td></tr><tr><td>New approach</td><td>0.1134</td></tr></tbody></table>

The computational time is less than 10 s for the new approach whereas it is greater than 4 h for the technique of Bernasconi and Papadopoulos [2] with the software Matlab® [11].

### 3.3. FEM analysis of a welded assembly

This section focuses on the analysis of the welded assembly shown in Figure 7. A cyclic overload- ing has been applied on the whole structure to test the fatigue resistance in the case of a mechan- ical loading greater than the maximum authorized one to investigate low-cycle fatigue. Experi- mental tests showed that failure always appears at the edge of the weld bead as shown in Figure 7 in the heat affected zone just after about 30,000 loading cycles.

Following the computational strategy developed by Giraud [12] and improved by Agard [13], the simulation applied in this paper is composed of three stages:

(1) calculation of the metallurgical and mechanical state after welding;

(2) computation of the stabilized cyclic state;

(3) estimation of the lifetime by means of the Fatemie-Socie model.

The finite element computation of the stabilized elastoplastic cyclic state of a welded assembly requires considering the residual effects of welding in the case of a S355 steel: the metallurgical composition field through the heat affected zone and the associated residual stresses. For this purpose, the thermo-metallurgical and mechanical simulation of welding has been carried out by means of the software Sysweld® [14]. The mesh is made of 194,312 quadratic tetrahedral finite elements containing four integration points. The welding parameters and the associated heat source are detailed in [15]. Figure 8 shows the evolution of the temperature computed during welding and the final distribution of the von Mises residual stress.

It is to be noted that the von Mises residual stress achieves more than 400 MPa in the corner near the edge of the weld bead. As far as metallurgical transformations are concerned, Figure 9 shows that the bainite and martensite phases are both produced in the heat affected zone, in particular near the failure area.
<figure><img src="image_10.png" /><figcaption>**Figure 8.** Distribution of the temperature (°C) during welding (left) and the von Mises residual stress (MPa) after welding (right).</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>**Figure 9.** Proportion of phases after welding: martensite (left) and bainite (right).</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>**Figure 10.** Distribution of the von Mises residual stress (MPa) at the end of a stabilized loading cycle (left) and the corresponding estimation of the lifetime $N_f$ (right).</figcaption></figure>

Once the thermo-metallurgical state after welding is obtained, the computation is restarted by applying the cyclic loading with Giraud's cyclic multi-phase elastoplastic model [16] until a stabilized cyclic state is achieved. As plotted in Figure 10, the residual stresses are highly reduced near the edge of the weld bead because of the plastic shakedown.
As far as the estimation of the lifetime, denoted $N_f$, is concerned, it can be obtained by means of the Fatemi–Socie critical plane model [1] for low carbon steels:

$$ \frac{\Delta\gamma}{2} \left( 1 + k \frac{\sigma_n^{\max}}{\sigma_y^{\text{ref}}} \right) = (1 + \nu) \frac{\sigma_f'}{E} (2N_f)^b + \frac{k}{2} (1 + \nu) \frac{\sigma_f'^2}{E\sigma_y^{\text{ref}}} (2N_f)^{2b} \\ + (1 + \nu_p) \varepsilon_f' (2N_f)^c + \frac{k}{2} (1 + \nu_p) \frac{\varepsilon_f' \sigma_f'}{\sigma_y^{\text{ref}}} (2N_f)^{b+c}, \quad (10) $$

where $k = 0.6$, $\sigma_f' = 1027$ MPa, $\varepsilon_f' = 0.322$, $\nu_p = 0.5$, $\nu = 0.3$, $E = 210,000$ MPa, $b = -1027$ and $c = -0.487$. $\sigma_y^{\text{ref}}$ is taken equal to the yield stress of the theoretical initial homogeneous S355 steel: 355 MPa. The formula (10) combines the maximum normal stress $\sigma_n^{\max}$ on the critical plane associated with the maximum shear strain amplitude $\Delta\gamma/2$.

For the numerical treatment of the multiaxial critical plane criterion of Fatemi and Socie, the technique developed in Section 2 is applied for the treatment of the mechanical history of each integration point with Matlab® [11]. It consists of finding the critical plane having the maximum normal stress $\sigma_n^{\max}$ for obtaining $N_f$ by solving (10) by means of the Newton–Raphson method.

The number of cycles $N_f$ related to each finite element in Figure 10 corresponds to the lowest lifetime of the four integration points. The computational time is about 10 min with the strategy developed in Section 2 whereas it is more than 20 days for the approach of Bernasconi and Papadopoulos [2]. It is to be noted that the results are identical.

## 4. Conclusions

The aim of this paper was to propose an efficient computational strategy for searching at critical planes and determining the associated shear amplitude.

For this purpose, a numerical method has been proposed in Section 2:

* A formulation based on the Tresca norm was developed for the determination of the critical planes without scanning all the possible planes.

* An algorithm in optimal time $O(nk)$ with $k \ll n$ was detailed for numerical efficiency.

In Section 3, the applications show the capability of the developed approach in analyzing strain gauge measurements or FEM results. For example, the new numerical method allows decreasing the CPU time significantly: from 20 days by scanning planes to 10 min for the FEM analysis of a welded structure with the criterion of Fatemi and Socie.

## Conflicts of interest

Authors have no conflict of interest to declare.

## References

[1] A. Fatemi, D. Socie, "A critical plane approach to multiaxial fatigue damage including out-of-phase loading", *Fatigue Fract. Eng. Mater. Struct.* **1** (1988), p. 149-165.

[2] A. Bernasconi, I. V. Papadopoulos, "Efficiency of algorithms for shear stress amplitude calculation in critical plane class fatigue criteria", *Comput. Mater. Sci.* **34** (2005), p. 355-368.

[3] B. Weber, B. Kenmeugne, J. C. Clement, J. L. Robert, "Improvements of multiaxial fatigue criteria computation for a strong reduction of calculation duration", *Comput. Mater. Sci.* **15** (1999), p. 381-399.

[4] K. Dang Van, B. Griveau, O. Message, "On a new multiaxial fatigue limit criterion: theory and applications", in *Biaxial and Multiaxial Fatigue*, EGF 3 (M. W. Brown, K. J. Miller, eds.), Mechanical Engineering Publications, London, 1989, p. 479-496.
[5] F. Fauvin, J. C. Roux, P. Monnet, E. Feulvarch, "Fast estimation of the shear stress amplitude for fatigue life analysis of metals", *Eur. J. Mech. - A/Solids* **8** (2020), article no. 103928.

[6] J. Lemaitre, J. L. Chaboche, *Mechanics of Solid Materials*, Cambridge University Press, Cambridge, UK, 1990.

[7] E. Ramos, "Intersection of unit-balls and diameter of a point set in R3", *Comput. Geom. Theory Appl.* **8** (1997), p. 57-65.

[8] K. L. Clarkson, P. W. Shor, "Applications of random sampling in computational geometry", *Discrete Comput. Geom.* **4** (1989), p. 387-421.

[9] S. Har-Peled, "A practical approach for computing the diameter of a point-set", in *17th ACM Symposium on Computational Geometry (SOCG'2001), Tufts Univ., Medford (MA), USA*, 2001, p. 177-186.

[10] F. P. Preparata, M. I. Shamos, *Computational Geometry*, Monographs in Computer Science, Springer, New York, 1985.

[11] Matlab®, *Mathematical Computing Software*, The MathWorks Inc., USA, 2022.

[12] L. Giraud, "Apports au dimensionnement en fatigue à faible nombre de cycles des assemblages soudés en S355", PhD Thesis, Univ. Lyon, 2020 (in French).

[13] B. Agard, "Détermination d'une stratégie de dimensionnement en fatigue à faible nombre de cycles adaptée au contexte industriel", PhD Thesis, Univ. Lyon, 2021 (in French).

[14] Software Sysweld®, *Reference Manual*, ESI Group, Lyon, France, 2022.

[15] E. Cottier, P. Anglade, A. Brosse, E. Feulvarch, "Fast 3D simulation of a single-pass steel girth weld", *Mech. Ind.* **17** (2016), article no. 401.

[16] L. Giraud, C. Pouvreau, F. Josse, W. Berckmans, F. Lefebvre, C. Carrillo, E. Feulvarch, "A multi-phase linear kinematic elastoplastic model for the HAZ of welded S355J2 steel under low-cycle fatigue", *C. R. Mécanique* **348** (2020), no. 3, p. 175-190.