# Shakedown within polycrystals: a direct numerical assessment

Domenico Magisano, Eric Charkaluk, Géry de Saxcé, Toufik Kanit

► To cite this version:

Domenico Magisano, Eric Charkaluk, Géry de Saxcé, Toufik Kanit. Shakedown within polycrystals: a direct numerical assessment. Advances in Direct Methods for Materials and Structures, 2018, 10.1007/978-3-319-59810-9_3. hal-01772204

HAL Id: hal-01772204

https://hal.science/hal-01772204

Submitted on 24 Sep 2022

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
# Shakedown Within Polycrystals:
A Direct Numerical Assessment

**D. Magisano**

Università della Calabria, Rende, Italy
e-mail: domenico.magisano@unical.it

**E. Charkaluk**

École Polytechnique, Palaiseau, France
e-mail: charkaluk@lms.polytechnique.fr

**G. de Saxcé** **T. Kanit**

Université Lille 1, Villeneuve-d'Ascq, France
e-mail: gery.desaxce@univ-lille1.fr
e-mail: toufik.kanit@univ-lille1.fr

**Abstract** It is well known that in high cycle fatigue (HCF), macroscopically, structures undergo elastic shakedown and the stress level commonly determines the lifetime. In this domain, the fatigue phenomena is due to local plasticity at the grain scale. Therefore, some multiscale HCF multiaxial fatigue criteria were proposed, among them the well-known Dang Van criterion. This criterion supposes that in a polycrystal, some misoriented grains can undergo plastic shakedown which conducts to crack initiation. The objective of this work is to validate this assumption by conducting numerical simulations on polycrystalline aggregates. As it is necessary to estimate the stabilized state in each grain of the polycrystal, classical incremental simulations are not the best way as it will be highly time-consuming because of the size of the aggregate. In the recent years, Pommier proposed a method called Direct Cyclic Algorithm to obtain the stabilized response of a structure under cyclic periodic loading, which it is shown to be more efficient compared to an incremental analysis in such situation. However, errors can be obtained in certain case with respect to the incremental solution. In this work, a Crystal Plasticity FEM model, based on dislocation densities, was used. As a first step, an aggregate of 20 grains of AISI 316L stainless steel under strain controlled cyclic loading was studied. Precise comparisons were conducted with incremental analysis and the results show that DCA seems to be an efficient solution in order to estimate the shakedown state of polycrystalline aggregates.
# 1 Introduction

The fatigue of metals under cyclic loadings is the consequence of crack initiation and growth until the complete failure of the concerned structure. The design of metal components requires physically based crack initiation criteria compatible with structural computation, performed for example by Finite Element Method. In this aim many theoretical developments were done in the last decades and this conducted to different fatigue criteria. High cycle fatigue (HCF) criteria are generally stress tensor based while low cycle fatigue (LCF) ones are often defined from strain variables [22]. However, if elementary mechanisms have been extensively studied, the complete understanding of transition from cyclic plasticity to strain localization and crack initiation is still an open problem.

A common approach in LCF and HCF and an improvement of such criteria can be researched by taking into account the crystalline microstructure and not only the macroscopic behaviour. In fact, metallic materials are made of an aggregate of grains more or less well-oriented, with respect to the loading axis. Under mechanical loading, this leads to heterogeneous deformation at the microstructure scale [10, 19]. In polycrystals with random crystallographic orientations, inhomogeneous stress and plastic strain fields are established because of orientation, grain shape and size, and geometrical effects. In LCF, bulk plastic deformation takes place, but this is by no means homogeneous. It is expected, therefore, that localized regions of preferential slip develop, leading to similarly localized regions of crack initiation. In HCF, the material undergoes elastic deformation at the macroscopic level, and it is only in small areas that plastic deformation occurs. For the case of a polycrystalline metal subjected to HCF with uniform macroscopic stress, the localized regions deforming plastically are generated again by the inhomogeneity of the stress resulting from crystallographic orientation, grain shape, etc. Even in HCF, as very localized plastic straining develops, hardening can occur, leading to a redistribution of stress and a further localization of plasticity. In other words, this is shakedown occurring, and it is the basis of the Dang Van criterion for crack initiation [3]. Dang Van argues that local plastic flow is essential for crack initiation and that, if shakedown occurs at the grain level, the material would have an infinite fatigue life.

Since many years, experimental tests and numerical simulations on polycrystalline aggregates are made in order to understand and model such mechanisms. Finite element models, incorporating crystal plasticity [8], were developed in order to study the evolution of mechanical quantities in the microstructure, the variables which rule the heterogeneous behaviour, for example the role of grain orientation with respect to the loading axis and misorientation with its neighbours in causing load shedding and stress localizations, and for trying to find a link between macroscopic and microscopic behaviour. Numerical simulations can be conducted on real microstructures, built by using an EBSD technique in order to obtain informations on grain boundaries and orientations, or on fictive microstructures, based on a random distribution of crystallographic orientations and grain sizes [16]. The simulations are then made by using a FE code, after the introduction
of the crystal plasticity models and the definition of boundary conditions and load-
ings. The classical incremental method (Newton-Raphson) is the standard technique
to solve this kind of nonlinear analysis. Due to the size of the aggregates and the
nonlinear character of the differential equations, such computations are very time-
consuming. Under cyclic loadings these simulations are quite not possible due to
the necessary high number of cycles before obtaining the stabilized response of the
aggregates. A method called Direct Cycle Analysis (DCA) has been recently imple-
mented in the commercial software ABAQUS [15] and it is today commonly used
to obtain the stabilized cycle for LCF design.

The aim of this work is to test the performances and the accuracy of the DCA
in the evaluation of the steady state of policrystalline aggregates under cyclic load-
ing. As a first step, an aggregate of 20 grains of AISI 316L stainless steel under
strain controlled cyclic loading is studied. A Crystal Plasticity FEM model based
on dislocation densities, already implemented in a User-subroutine for ABAQUS/
Standard [20], is considered. The fictive microstructure and the finite element mesh
are generated by NEPER, a software based on Voronoi tessellation [16]. This work
starts with an introduction to the behaviour of metals under cyclic loads. The crys-
tal plasticity model is then described together with how to obtain the geometrical
model of a microstructure. The DCA is briefly recalled and it is tested in the case of
crystal plasticity. Precise comparisons are conducted with an incremental analysis in
order to validate the DCA as an efficient and accurate tool for the evaluation of the
shakedown state of polycrystalline aggregates.

2 Steady State and Fatigue Criteria

Among different methods of fatigue design of metal structures, a decoupled approach
can be used. It plans to perform a thermomechanical analysis without taking into
account the damage and then appropriate fatigue criteria allow to link the mechan-
ical parameters of the stabilized state to the lifetime of the structure. A structure
whose material obeys Drucker's postulate, will reach, after some cycles of loading,
a steady cycle in which the stresses and the strain rates gradually stabilize and remain
unaltered on passing to the next cycle. There are three different categories of steady
stress cycles:

* Ratcheting (or incremental collapse), where non-vanishing plastic strain rates $\dot{\varepsilon}^p$ and a non-vanishing plastic strain ratchet $\Delta\varepsilon^p \neq 0$ occurs at various parts of V. This happens for sufficiently high load amplitudes and is a dangerous long-term response as the plastic strains grow bigger from cycle to cycle leading to large displacements so that the structure becomes unserviceable or the ultimate strain of the material is reached.

* Alternating plasticity (or plastic shakedown), in which non-vanishing plastic strain rates $\dot{\varepsilon}^p$ exist but there is no increment of plastic strains, i.e. $\Delta\varepsilon^p = 0$. This also
occurs for high load amplitudes and this type of long-term response leads to low cycle fatigue, which reduces the working life of a structure or of a component.

* Shakedown (or elastic shakedown), where plastic strain rate $\dot{\epsilon}^p$ vanish over the whole body. In this case further plastic straining stops when the steady cycle is reached and the structure subsequently responds, to further cycling, purely elastically. This long-term response occurs for relative low load amplitudes and is a favourable situation provided the plastic deformation that has been produced at the transient phase is sufficiently small. In this case the structure or a structural component can have infinite fatigue life, i.e. there is no failure even after many cycles, or finite life due to high cycle fatigue.

In the case of elastic shakedown the structure, suddenly or after some cycles, has an elastic response. Even if the stress is less than yield stress failure can occurs because of HCF. The HCF design is based on the stress control and the stress is computed by considering the metal as homogeneous, without taking into account stress concentration due to the material heterogeneity. Some metals, if the stress level is low, show an infinite life, i.e. fatigue does not occur even if the number of cycles is very high. When the structure undergoes plastic shakedown, fatigue occurs after a lower number of cycle and the LCF design is usually strain controlled. Both criteria are purely phenomenological, relying directly on the interpretation of experimental results at the macroscopic scale. An example of fatigue criterion that takes into account the difference between the macroscopic and the mesoscopic mechanical fields is the Dang Van criterion for infinite life.

## 2.1 Dang Van Criterion

One of the fatigue models including grain level phenomena in a macroscopic fatigue criterion is the Dang Van-Papadopulos criterion based on the following assumption:

* In the LCF regime, physical observations at both macroscopic and mesoscopic scale show extensive plastic strains. Moreover homogenisation theory shows that strains and stresses at the two scales tend to be closer to each other with increasing plastic strain. This can be translated into saying that the higher the applied load, the more similar mesoscopic and macroscopic scales will behave.

* In the HCF regime, two fatigue domains corresponding to finite and infinite lifetime can be considered. Physical observations at the macroscopic scale show that structures are macroscopically in an elastic shakedown state. At the mesoscopic scale of the grains, it is now commonly accepted that elastic shakedown occurs only in the case of infinite lifetime. If lifetime is finite, some grains will be oriented such that they can not reach an elastic shakedown state, but will experience a plastic shakedown or ratcheting state leading to failure after a finite number of cycles. The stress concentration due to this mesoscopic failure marks the initiation of a macroscopic crack associated with failure on the macroscopic scale.
Focusing on the case of HCF, one can imagine a case where only one misoriented grain is subject to plastic slip. Then a simple homogenisation scheme of a plastic inclusion in an elastic matrix can be used to derive closed-form relations between mesoscopic and macroscopic fields. Examples of possible homogenisation assumptions are:

* Lin–Taylor strain equality: $\epsilon = E$, that is the hypothesis of the initial Dang Van–Papadopoulos fatigue criterion;

* Sachs stress equality $\sigma = \Sigma$;

where $\epsilon$ and $\sigma$ are strain and stress at grain level, $E$ and $\Sigma$ are strain and stress at macroscopic level.

If, in all the cases, the same elastic behaviour at the mesoscopic and the macroscopic scale is assumed, the relation between mesoscopic and macroscopic fields can be written in the general form:

$$ \sigma = \Sigma - C^* \epsilon^p = \Sigma + \rho^* \quad (1) $$

where $\rho^*$ should be interpreted as mesoscopic residual stress field. The particular case of each homogenisation model is obtained depending on the form of $C^*$:

* Lin–Taylor’s model: $C^* = C$

* Sachs’s model $C^* = 0$;

As we will see better in the next chapter the plastic deformation is the consequence of atomic slip in preferential directions on particular atomic planes (slip systems). However under the assumption that the grain orientations statistically cover all directions [14], and so even the slip systems, according to Dan Vang fatigue criterion infinite life occurs if

$$ \max_t(\tau^T(t) + a\sigma^H(t)) < b \quad (2) $$

where $\tau^T$ is the Tresca norm of mesoscopic shear, $\sigma^H = 1/3 \text{tr}(\sigma)$ is the hydrostatic mesoscopic stress, $a$ and $b$ are material parameters.

The homogenisation assumptions seen previously cannot represent at the best the reality. Numerical simulations on polycrystalline aggregates can be the way to improve the link between mesoscopic and macroscopic quantities.

# 3 The Crystal Plasticity Model

## 3.1 Single Crystal Plasticity

In crystal plasticity theory, plastic deformation is modelled by using the slip system activity concept. It is a physically based plasticity theory that represents the deformation of a metal at the microscale. It is important to note that, in this theory,
the continuum framework is kept even if at such scale the physics is not continuous anymore. Dislocations are assumed to move across the crystal lattice along specific crystallographic planes and directions. When the material is subjected to loading, the local resolved shear stress occurs on a slip plane and along a slip direction and its magnitude controls the movement, the creation and the annihilation of dislocations. Consequently the material is locally loaded on specific directions while the volume remains constant. Moreover the crystal lattice could deform elastically while the resolved shear stress has not reached its critical level. Two laws are necessary to describe the single crystal plasticity:

* the flow rule, that describes, for intragranular variable, the slip initiation for each system;

* the hardening rule, that describes the hardening which occurs on each system after successive loadings.

The crystal plasticity model described in the following is that already implemented as a User-subroutine for ABAQUS/Standard [20].

### 3.1.1 Flow Rule

The Resolved Shear Stress vector $\tau^s$ is the projection of the local stress tensor expressed in the global reference system on every possible slip system $s$. The Schmid law takes the following general form [17]:

$$ \tau^s = \sigma : D^s \quad (3) $$

$$ D^s = \frac{1}{2}(b^s \otimes n^s + n^s \otimes b^s) \quad (4) $$

where $\sigma$ is the local stress tensor, $n^s$ is the slip plane normal, $b^s$ is the slip direction and $D^s$ is the Schmid matrix, that is the symmetric part of the Schmid tensor $L^s = b^s \otimes n^s$ and we denote here $W^s$ its antisymmetric part. One can easily find that, in the case of a single crystal submitted to a uniaxial loading, the Schmid matrix $D^s$ is reduced to a simple composition of crystallographic orientations $cos\psi cos\lambda$, where $\lambda$ is the slip direction angle related to the load direction and $\psi$ is the slip plane normal angle. According to the Schmid law, a single crystal leaves the elasticity and shear occurs when the Resolved Shear Stress on a slip system $s$ reaches its threshold value, the Critical Resolved Shear Stress (CRSS), noted $\tau_c^s$. It could basically be written as

$$ \left\{ \begin{array}{l} \dot{\gamma}^s = 0, \tau^s < \tau_c^s \\ \dot{\gamma}^s > 0, \tau^s \ge \tau_c^s \end{array} \right. \quad (5) $$

with $\dot{\gamma}^s$ the slip rate associated to system $s$. The flow rule introduces a rate dependent formulation. Classically one use a viscoplastic formulation more convenient for numerical simulation [9], that consists in a power law:
$$ \dot{\gamma}^s = \dot{\gamma}_0 \left| \frac{\tau^s}{\tau_c^s} \right|^n \text{sign}(\tau^s) \quad (6) $$

where $\dot{\gamma}_0$ is a reference value of the slip rate and $n$ is the strain rate sensibility parameter. The rate dependent approach in fact avoids to deal with discontinuous conditions because all slip systems are allowed to be active but only some systems have non-negligible slips (when $n \gg 1$, $\dot{\gamma}^s \cong 0$ unless $\tau^s \cong \tau_c^s$).

### 3.1.2 Hardening Rule

When sliding occurs on a specific system, it interacts with different obstacles: additive elements, precipitates, dislocations “forest”, etc. The increasing of dislocation density and these interactions lead to a local material resistance which results in critical shear stress increasing. In some cases an annihilation process could occurs leading to a decrease of dislocation density storage rate. One experimentally observes that there is a critical distance between dislocation systems with opposite sign leading to their annihilation. This process could conduct to a saturation process of hardening and also softening process. A physical based hardening rule has to take into account these basic mechanisms of dislocations generation and annihilation. It can take the following general expression proposed by [11]:

$$ \tau_c^s = \sum_u H^{su} |\dot{\gamma}^s| \quad (7) $$

where the $H^{su}$ terms are the components of the hardening matrix. This matrix takes into account the slip on the considered system also on the others. Diagonal terms $H^{ss}$, called self-hardening, account for the hardening on slip system $s$ due to its own slip activity. Remaining terms, called latent-hardening, account for the hardening on slip system $s$ due to the slip activity on the whole set of other systems $u$. The form of the hardening matrix depends on the considered physical mechanisms.

To account for more physics Eq. 7 should be based on dislocation density. Basically, the total dislocation density could be divided in two parts: statistically stored dislocations $\rho_s^s$, which trap each others in a random way and/or are required for compatible deformation of various part of the crystal, and geometrically stored dislocations $\rho_g^s$, which are required when gradient of plastic shear exists. For simplicity's sake in this model it's implicitly assumed that only $\rho_s^s$ influences the work hardening. No geometrically stored dislocations are therefore considered. The form of the slip resistance classically used was introduced in [5]:

$$ \tau_c^s = \tau_0 + \mu b \sqrt{\sum_u d^{su} \rho^u} \quad (8) $$

where $\tau_0$ is an initial value of Critical Resolved Shear Stress, $\mu$ is the shear modulus, $b$ the magnitude of the Burgers vector, $d^{su}$ the interaction coefficient between
the slip systems *s* and *u* and $\rho^u$ represents the local density of the statistically stored dislocations on the slip system *u*. In the case of FCC structure the hardening interaction matrix is a 12 × 12 matrix composed by only 4 independent coefficients due to crystal symmetries. The identification of these coefficients is still complicated and is classically done by Discrete Dislocation Dynamics simulations. In practice one uses here a simple form of the matrix depending on only two parameters: $d_p$ for the self-hardening ($s = u$) and $d_f$ for the latent-hardening ($s \neq u$).

The hardening rule needs to be completed by a flow rule for dislocation density on each slip system *s*, function of slip rate $\dot{\gamma}^s$. The description of dislocation generation was firstly figured out by Frank and Read [6]. It takes generally the following form:

$$ \dot{\rho}^s = \frac{1}{b} \left( \frac{\sum_u (a^{su} \rho^u)}{L} - 2y_c \rho^s \right) |\dot{\gamma}^s| \quad (9) $$

in which the term $\frac{\sum_u (a^{su} \rho^u)}{Lb}$ accounts for dislocation creation and the term $2y_c \rho^s$ accounts for dislocation annihilation.

The material parameter $2y_c$ associated to annihilation is related to the critical annihilation distance, which is here taken constant. The parameter *L* is the mean free path of the mobile dislocations in system *s*, also assumed to be constant in this study. Finally $a^{su}$ accounts for the geometrical interaction between dislocation densities. Assuming that collinear, coplanar and orthogonal crystal interactions have the same impact on mean free paths allows to reduce the number of $a_i$ unknown to only four independent coefficients as shown in Table 1. Note that the Eq. 9 implies that there exists a saturation of dislocation production which depends on the deformation path, and more precisely on the dislocation quantity accumulated on the whole set of slip systems. Note also that dislocation production rate $\dot{\rho}^s$ depends on the absolute value

<table><caption>**Table 1** Dislocation density interaction matrix of FCC material</caption><thead><tr><th></th><th>A<sub>2</sub></th><th>A<sub>3</sub></th><th>A<sub>6</sub></th><th>B<sub>2</sub></th><th>B<sub>4</sub></th><th>B<sub>5</sub></th><th>C<sub>1</sub></th><th>C<sub>3</sub></th><th>C<sub>5</sub></th><th>D<sub>1</sub></th><th>D<sub>4</sub></th><th>D<sub>6</sub></th></tr></thead><tbody><tr><th>A<sub>2</sub></th><td>a<sub>1</sub></td><td>a<sub>3</sub></td><td>a<sub>3</sub></td><td>a<sub>3</sub></td><td>a<sub>5</sub></td><td>a<sub>5</sub></td><td>a<sub>3</sub></td><td>a<sub>5</sub></td><td>a<sub>6</sub></td><td>a<sub>3</sub></td><td>a<sub>6</sub></td><td>a<sub>5</sub></td></tr><tr><th>A<sub>3</sub></th><td></td><td>a<sub>1</sub></td><td>a<sub>3</sub></td><td>a<sub>5</sub></td><td>a<sub>3</sub></td><td>a<sub>6</sub></td><td>a<sub>5</sub></td><td>a<sub>3</sub></td><td>a<sub>5</sub></td><td>a<sub>6</sub></td><td>a<sub>3</sub></td><td>a<sub>5</sub></td></tr><tr><th>A<sub>6</sub></th><td></td><td></td><td>a<sub>1</sub></td><td>a<sub>5</sub></td><td>a<sub>6</sub></td><td>a<sub>3</sub></td><td>a<sub>6</sub></td><td>a<sub>5</sub></td><td>a<sub>3</sub></td><td>a<sub>5</sub></td><td>a<sub>5</sub></td><td>a<sub>3</sub></td></tr><tr><th>B<sub>2</sub></th><td></td><td></td><td></td><td>a<sub>1</sub></td><td>a<sub>3</sub></td><td>a<sub>3</sub></td><td>a<sub>3</sub></td><td>a<sub>6</sub></td><td>a<sub>5</sub></td><td>a<sub>3</sub></td><td>a<sub>5</sub></td><td>a<sub>6</sub></td></tr><tr><th>B<sub>4</sub></th><td></td><td></td><td></td><td></td><td>a<sub>1</sub></td><td>a<sub>3</sub></td><td>a<sub>6</sub></td><td>a<sub>3</sub></td><td>a<sub>5</sub></td><td>a<sub>5</sub></td><td>a<sub>3</sub></td><td>a<sub>5</sub></td></tr><tr><th>B<sub>5</sub></th><td></td><td></td><td></td><td></td><td></td><td>a<sub>1</sub></td><td>a<sub>5</sub></td><td>a<sub>5</sub></td><td>a<sub>3</sub></td><td>a<sub>6</sub></td><td>a<sub>5</sub></td><td>a<sub>3</sub></td></tr><tr><th>C<sub>1</sub></th><td></td><td></td><td></td><td></td><td></td><td></td><td>a<sub>1</sub></td><td>a<sub>3</sub></td><td>a<sub>3</sub></td><td>a<sub>3</sub></td><td>a<sub>5</sub></td><td>a<sub>5</sub></td></tr><tr><th>C<sub>3</sub></th><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td>a<sub>1</sub></td><td>a<sub>3</sub></td><td>a<sub>5</sub></td><td>a<sub>3</sub></td><td>a<sub>6</sub></td></tr><tr><th>C<sub>5</sub></th><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td>a<sub>1</sub></td><td>a<sub>5</sub></td><td>a<sub>6</sub></td><td>a<sub>3</sub></td></tr><tr><th>D<sub>1</sub></th><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td>a<sub>1</sub></td><td>a<sub>3</sub></td><td>a<sub>3</sub></td></tr><tr><th>D<sub>4</sub></th><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td>a<sub>1</sub></td><td>a<sub>3</sub></td></tr><tr><th>D<sub>6</sub></th><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td>a<sub>1</sub></td></tr></tbody></table>
of slip rate $\dot{\gamma}^s$ and this hardening is therefore isotropic. Nevertheless, as each slip direction has its own dislocation and yield evolution, the yield surface in spatial domain evolves from an initial sphere of radius $\tau_0$ to a more complex form (Fig. 1).

## 3.2 Microstructure Geometry

Numerical simulation can be computed on two different kind of polycrystalline aggregates:

* real microstructure, built by E.S.B.D. analysis which provides information on grain boundaries and crystallographic orientation

* fictive microstructure, based on a random distribution of crystallographic orientations and grain sizes.

### 3.2.1 Fictive Microstructures

Several authors have proposed analytical and numerical methods to construct random polycrystalline morphologies. These methods can be rooted in the basic principles of phase transformations (leading in some cases to Voronoi tessellations), physically-driven simulations of annealing or recrystallization, or algorithms that attempt to directly reproduce statistical data coming from experimental characterizations. Among randomly generated morphologies, Voronoi tessellations have the advantages of being defined analytically and having straight triple lines and flat grain

<figure><img src="image_2.png" /><figcaption>Fig. 1 Hardening evolution with dislocation density</figcaption></figure>
boundaries. Still, Voronoi tessellations show important variabilities in grain size and shape which are representative of real polycrystalline morphologies [16].

Mathematically a Voronoi tessellation of a n-D space is a collection of n-D entities that fills the space with no overlaps and no gaps. These entities are polyhedra and are formally defined as zones of influence of a particular set of points, corresponding to their centres. Being given a spatial domain $D \in R^n$, a set of points $G_i(x_i)$ within D and a norm $d(\cdot, \cdot)$, a Voronoi polyhedron $C_i$ is associated to every point $G_i$ as follows:

$$C_i = \{ P(x) \in D | d(P, G_i) \le d(P, G_j) \quad \forall j \ne i \} \quad (10)$$

There are different available software which provide a Voronoi tessellation of 2-D and 3-D domains, like for example NEPER [16], which provides the polycrystalline morphology, described by sets of points, lines, surfaces and volumes, and the free meshing of the morphology. It uses an Euclidean distance like norm and the set of points $G_i$ is taken randomly distributed. By construction, a Voronoi polyhedron is convex; hence in 3-D the intersection of two Voronoi polyhedra is a plane, called “tessellation face”, the intersection of three Voronoi polyhedra is a straight line, called “tessellation edge”, and the intersection of four Voronoi polyhedra is a point, called “tessellation vertex”. From a physical point of view, the generation of Voronoi tessellations corresponds to a process of solidification or recrystallization where all grains nucleate at the same time and grow isotropically at the same rate. Voronoi tessellations qualitatively reproduce some important properties of real polycrystalline morphologies, like the distribution of grain size and the number of first neighbours.

## 4 Review of the DCA

The classical approach in FE codes to obtain the stabilized response of an elastic-plastic structure subjected to cyclic loading is to apply the periodic loading cycles repetitively to the unstressed structure until a stabilized state is obtained. At each instant in time it typically involves using Newton's method to solve the nonlinear equilibrium equations. To avoid the considerable numerical expense associated with such a transient analysis, a DCA has been suggested in [15] and has been recently implemented in the commercial software ABAQUS. It is based on the following assumption:

* quasi-static analysis, i.e. without taking into account dynamic effect

* geometrically linear behavior

* use of Fourier series to describe the stabilized state

* direct research of the stabilized cyclic response of the structure iteratively using the Modified Newton's method with the elastic stiffness matrix
## 4.1 Description of the Algorithm

In a quasi-static analysis with cyclic periodic loads the nonlinear equilibrium equations can be written as:

$$R(t) = F(t) - I(t) = 0 \qquad (11)$$

where $F(t)$ is the discretized form of a cyclic load that has the characteristic $F(t + T) = F(t)$ at all times $t$ during a load cycle with period $T$, $I(t)$ represents the internal force vector generated by the stress, and $R(t)$ is the residual vector. We are looking for a displacement function that describes the stabilized response of the structure, i.e. a displacement function $u(t)$ that at all times $t$ during a load cycle with period $T$ has the characteristic $u(t) = u(t + T)$. Since the periodicity of the solution, the best choice for this purpose seems to be a truncated Fourier series:

$$u(t) = u_0 + \sum_{k=1}^{N} [u_k^s \sin k\omega t + u_k^c \cos k\omega t] \qquad (12)$$

where $n$ stands for the number of terms in the Fourier series; $\omega = 2\pi/T$ is the angular frequency; and $u_0, u_k^s$ and $u_k^c$ are unknown displacement coefficients.

We also expand the residual vector in a truncated Fourier series in the same form as the displacement solution:

$$R(t) = R_0 + \sum_{k=1}^{N} [R_k^s \sin k\omega t + R_k^c \cos k\omega t] \qquad (13)$$

where each residual vector coefficient $R_0, R_k^s$ and $R_k^c$ in the Fourier series corresponds to a displacement coefficient. At each instant in time in the cycle the residual vector $R(t)$ is obtained by using standard element-by-element calculations and provides the Fourier coefficients:

$$\begin{aligned} R_0 &= \frac{2}{T} \int_0^T R(t) dt \\ R_k^s &= \frac{2}{T} \int_0^T R(t) \sin(k\omega t) dt \\ R_k^c &= \frac{2}{T} \int_0^T R(t) \cos(k\omega t) dt \end{aligned} \qquad (14)$$

For the computing of the residual Fourier terms we need to evaluate the residual vector at some time points obtained by a discretization of time period in time increments. ABAQUS uses a trapezoidal rule, which assumes a linear variation of the residual over a time increment, to integrate the residual coefficients. For accurate integration the number of time points must be larger than the number of Fourier coefficients (which is equal to $2n + 1$, where $n$ represents the number of Fourier terms).
The unknown displacement coefficients can now be computed iteratively by using the Modified Newton method. Starting from their approximations obtained after an iteration (i), the coefficient after an iteration (i + 1) are:

$$
\begin{align}
u_0^{(i+1)} &= u_0^{(i)} + c_0^{(i+1)} \nonumber \\
u_k^{(i+1)} &= u_k^{(i)} + c_k^{(i+1)} \tag{15} \\
u_k^{c(i+1)} &= u_k^{c(i)} + c_k^{c(i+1)} \nonumber
\end{align}
$$

where the corrections to the coefficients of the displacement solution are found by
solving the following linear systems:

$$
\begin{align}
K c_0^{(i+1)} &= R_0^{(i)} \notag \\
K c_k^{(i+1)} &= R_k^{(i)} \tag{16} \\
K c_k^{c(i+1)} &= R_k^{c(i)} \notag
\end{align}
$$

with *K* the initial elastic stiffness matrix. The updated displacement coefficients are used in the next iteration to obtain displacements at each instant in time. This process is repeated until convergence is obtained. Each pass through the complete load cycle can therefore be thought of as a single iteration of the solution to the nonlinear problem [1]. Since the use of the Modified Newton iterative scheme, the DCA cannot be utilized in the presence of strong nonlinearities, such as geometric nonlinearities and contact.

During the iterative solution process, the DCA imposes periodic conditions by using the state obtained at the end of the previous iteration as the starting state for the current iteration, i.e. $s_{i=0}^{i+1} = s_{i=T}^i$, where *s* is a solution variable such as plastic strain.

**5 Numerical Simulations on a Polycrystal**

**5.1 Geometric Model and Material**

In this section we analyse a polycrystal composed by 20 grains. The specimen is a
square cuboid in which two dimensions are ten times the other one. The microstruc-
ture, the crystal orientations and the mesh are generated by NEPER. The geometrical
model and the boundary conditions are shown in Fig. 2. The model is composed by
2094 tetrahedral elements C3D4 and 636 nodes. The material used in the simulation
is the AISI 316L austenitic stainless steel, a material used for many structural appli-
cations. It is a polycrystalline aggregates with FCC structure and its crystal plasticity
parameters are available in literature and are reported below. The definition of elastic
fourth rank stiffness tensor C for a FCC crystal need only 3 parameters, due to the
<figure><img src="image_3.png" /><figcaption>**Fig. 2** Geometrical model and boundary condition of the polycrystal</figcaption></figure>

cubic symmetry, which for 316 L are shown in Table 2, where $a = 2C_{44}/(C_{11} - C_{22})$ is the coefficient of anisotropy ($a = 1$ for isotropic material). The flow rule consists in a power law defined by 2 parameters reported in Table 3: the reference value of the slip rate $\dot{\gamma}_0$ and the strain rate sensitivity parameter. Eleven parameters are required to define the hardening behaviour of an FCC crystal and are summarized in Tables 4 and 5.

<table><caption>**Table 2** Elastic parameters of CPFEM for 316L [13]</caption><thead><tr><th>C<sub>11</sub> (MPa)</th><th>C<sub>12</sub> (MPa)</th><th>C<sub>44</sub> (MPa)</th><th>a</th></tr></thead><tbody><tr><td>248000</td><td>142000</td><td>71000</td><td>1.34</td></tr></tbody></table>

<table><caption>**Table 3** Flow parameters of CPFEM for 316L [2, 21]</caption><thead><tr><th>γ̇<sub>0</sub> (s<sup>-1</sup>)</th><th>n</th></tr></thead><tbody><tr><td>3.4 × 10<sup>-6</sup></td><td>10</td></tr></tbody></table>
<table><caption>**Table 4** Hardening parameters of CPFEM for 316L [4, 7]</caption><thead><tr><th>tau<sub>0</sub> (MPa)</th><th>b (mm)</th><th>y<sub>c</sub></th><th>L (mm<sup>-2</sup>)</th><th>ρ<sub>0</sub></th><th>d<sub>p</sub></th><th>d<sub>f</sub></th></tr></thead><tbody><tr><td>35</td><td>3.2 × 10<sup>-7</sup></td><td>3.2 × 10<sup>-6</sup></td><td>33</td><td>10<sup>-3</sup></td><td>0.06</td><td>0.004</td></tr></tbody></table>

<table><caption>**Table 5** Dislocation density interaction parameters of CPFEM for 316L [12]</caption><thead><tr><th>a<sub>1</sub></th><th>a<sub>3</sub></th><th>a<sub>5</sub></th><th>a<sub>6</sub></th></tr></thead><tbody><tr><td>0.12</td><td>0.07</td><td>0.14</td><td>0.12</td></tr></tbody></table>

## 5.2 Monotonic Tensile Test

A strain controlled monotonic tensile test with a constant strain rate of 0.005 s⁻¹ is realized. The simulation is first performed without taking into account the geometric nonlinearities, i.e. the rotation of crystal axis during the deformation.

Figure 3 shows the heterogeneous mechanical average behaviour of the different grains in loading direction. At 2.5% of average strain the heterogeneity on strain is about 1% and about 144 MPa on stress. We can see on the average strain vs stress curve that, at 2.5% of average strain, the stress is 263 MPa and is consistent with experimental and other numerical results [18]. The behaviour obtained is qualitatively acceptable although the number of grains is very low and the microstructure

<figure><img src="image_4.png" /><figcaption>**Fig. 3** Monotonic test, stress-strain curve in load direction for each grain</figcaption></figure>
is fictive. The influence of the mesh on the results is not tested. The same test is per-
formed with both models, geometrically linear and nonlinear, up to 5% of average
strain in order to quantify the error of neglecting the rotation of crystal axis during
the deformation. In Fig. 4 we can see that, even at 5% of average strain, the error
on the strain is certainly negligible, while the error on the average strain and on the
maximum stress becomes to be of few percent. However in fatigue design the values
of strain reached are usually much lower that 5% and so the neglecting of geomet-
ric nonlinearities are surely acceptable. After this observation we can safely use the
Direct Cyclic Algorithm on polycrystalline aggregates. In fact it cannot be used for
geometrically nonlinear analysis since it uses the modified Newton-Raphson method
to solve the equilibrium equations.

5.3 Cyclic Test

A strain controlled cyclic symmetric tensile-compression test up to 0.25% with a
constant strain rate of 0.0025 s⁻¹ is now performed on the same aggregate and the
results are discussed. Figures 5 and 6 show the average behaviour and the behaviour

<figure><img src="image_5.png" /><figcaption>Fig. 4 Error of neglecting geometric nonlinearities on maximum stress, maximum strain and average stress in load direction</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>**Fig. 5** Stress-strain curve in load direction for the misoriented grain: direct cyclic and incremental (cycle 1, 20, 100, 200, 600 and 1000) solution</figcaption></figure>

of the misoriented grain. We can see that the steady cycle is reached by an incremental analysis about 600 cycles. The DCA requires with 19 terms, the maximum number of Fourier terms, 682 iterations to converge. We can also see that the results obtained by DCA and Incremental analysis are roughly equivalent by looking at the two stress-strain curves.

Figure 7 shows that the strain field obtained with the two method has practically the same form and the same values. Looking at Fig. 8 we can see that even with only 5 Fourier terms and practically the same number of iterations we can obtain a solution very close to the incremental one. In fact, if we exclude the decomposition of the stiffness matrix and the evaluation of the residual vector in each time increment, which are the same in both analysis, we can link the remaining computational time to the number of linear systems required to achieve the convergence:

* 19 F.T. ⇒ 682 iterations × (2 × 19+1) F.C. = 26598 L.S.

* 5 F.T. ⇒ 682 iterations × (2 × 5+1) F.C. = 7480 L.S.

In this case the number of linear systems with 5 Fourier terms is 3.5 times lower than the number required with 19. The error is quantified in Tables 6 and 7 for the mis-oriented grain and for the average behaviour for different number of terms. Looking to these results we can make some consideration:
<figure><img src="image_7.png" /><figcaption>**Fig. 6** Stress-strain average curve in load direction: direct cyclic and incremental (cycle 1, 20, 100, 200, 600 and 1000) solution</figcaption></figure>

* the DCA solution seems to converge to the solution with more Fourier terms;

* at convergence, the DCA solution is not coincident with the incremental one and then a lower number of Fourier terms can provide better results.

In Fig. 9 it is shown the behaviour of two different grains with respect to the average behaviour. The load transfer between groups of grains of different orientation occurs due to plastic deformation governed by crystal slip. In this respect the grain 7 and 8 in this polycrystalline aggregates are representative of the extremes of behaviour, grain 7 being the strongest orientation, and grain 8 being the softest and the most compliant. It is interesting to note that grain 7 undergoes practically elastic shakedown, while the response of grain 8 exhibits a significant hysteresis area. This is indicative of inelastic deformation by crystal slip, which leads to energy dissipation. It is interesting to interpret this result within the framework of energy-related criteria of fatigue crack initiation. For example, Dang Van postulates that infinite life may be expected, but only if after a certain initial period grain level shakedown occurs, the cyclic deformation becomes purely elastic, and no energy dissipation takes place. The result shown in Fig. 9 indicates that the situation is different between grains even within the same macroscopic sampling volume. Some grains, e.g. grain 8, satisfy the initiation criterion, while others, e.g. grain 7, must be expected to have infinite fatigue life (Table 8).
<figure><img src="image_8.png" /><figcaption>(a) Incremental</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>(b) Direct Cyclic</figcaption></figure>

<caption>**Fig. 7** Strain field in load direction at 0.25% of average axial strain at steady cycle</caption>
<figure><img src="image_10.png" /><figcaption>**Fig. 8** Comparison of the stress-strain curve in load direction of the misoriented grain and the average behaviour obtained by DCA with 19 and 5 Fourier terms</figcaption></figure>

<table><caption>**Table 6** Direct cyclic error for the misoriented grain</caption><thead><tr><th>Analysis</th><th>ΔE</th><th>err (%)</th><th>ΔS (MPa)</th><th>err (%)</th><th>Energy (MPa)</th><th>err (%)</th></tr></thead><tbody><tr><td>D.C. 5 F.T.</td><td>0.0060</td><td>-4.8</td><td>599</td><td>0.3</td><td>1.710</td><td>1.5</td></tr><tr><td>D.C. 9 F.T.</td><td>0.0062</td><td>-1.7</td><td>602</td><td>0.9</td><td>1.753</td><td>4.0</td></tr><tr><td>D.C. 14 F.T.</td><td>0.0063</td><td>-0.2</td><td>605</td><td>1.4</td><td>1.754</td><td>4.1</td></tr><tr><td>D.C. 19 F.T.</td><td>0.0064</td><td>1.3</td><td>597</td><td>0.0</td><td>1.734</td><td>2.9</td></tr><tr><td>Incremental</td><td>0.0063</td><td></td><td>597</td><td></td><td>1.685</td><td></td></tr></tbody></table>

<table><caption>**Table 7** Direct cyclic error for the average behaviour</caption><thead><tr><th>Analysis</th><th>ΔE</th><th>err (%)</th><th>ΔS (MPa)</th><th>err (%)</th><th>Energy (MPa)</th><th>err (%)</th></tr></thead><tbody><tr><td>D.C. 5 F.T.</td><td>0.0046</td><td>-6.2</td><td>625</td><td>-2.4</td><td>0.743</td><td>-2.0</td></tr><tr><td>D.C. 9 F.T.</td><td>0.0048</td><td>-3.2</td><td>639</td><td>-0.2</td><td>0.761</td><td>0.4</td></tr><tr><td>D.C. 14 F.T.</td><td>0.0049</td><td>-1.6</td><td>640</td><td>0.0</td><td>0.761</td><td>0.4</td></tr><tr><td>D.C. 19 F.T.</td><td>0.0050</td><td>0.0</td><td>642</td><td>0.2</td><td>0.753</td><td>-0.6</td></tr><tr><td>Incremental</td><td>0.0050</td><td></td><td>640</td><td></td><td>0.758</td><td></td></tr></tbody></table>
<figure><img src="image_11.png" /><figcaption>**Fig. 9** Comparison of the stress-strain curve in load direction of different grains with respect to the average one</figcaption></figure>

<table><caption>**Table 8** Comparison of the behaviour of different grains with respect to the average one</caption><thead><tr><th>Variable</th><th>Average</th><th>Grain 7</th><th>Grain 8</th></tr></thead><tbody><tr><td>ΔE</td><td>0.0050</td><td>0.0043</td><td>0.0064</td></tr><tr><td>ΔS</td><td>642</td><td>784</td><td>597</td></tr><tr><td>Energy</td><td>0.75</td><td>0.19</td><td>1.73</td></tr></tbody></table>

# 6 Conclusions

Fatigue criteria based on stress and strain control are commonly used to link the steady response to the number of cycles that provides failure. In HCF the structure undergoes elastic shakedown and the life time is commonly determined by the macroscopic stress level. In LCF the structure undergoes alternating plasticity and the lifetime is estimated with the difference of strain of the steady cycles. However, a metal is a polycrystalline aggregate made of many grains with different orientations of the crystals. A heterogeneous response is expected due to the anisotropy of each crystal with regard to the elastic and, much more, the plastic behaviour. An improvement of fatigue criteria could be made by considering the behaviour of metals at the microstructure scale. To understand and to model such heterogeneity Crystal
Plasticity theory was used in the last years to build Crystal Plasticity FEM models. In this work a model based on dislocation density and large deformation was considered. The numerical model was already implemented as a UMAT for ABAQUS standard. We have seen that this model allows to achieve a first approximation of the heterogeneity of an aggregate of 20 grains of AISI 316L stainless steel under strain controlled monotonic and cyclic loading. The aggregate was built with the software NEPER that provides a random distribution of grain size and orientations and at first the analysis was performed incrementally with a full Newton method. This approach is quite expensive in the case of cyclic loading, since the application of many loading cycles may be required before the stabilized response is obtained, about 600 cycles for our aggregate.

Recently, Pommier proposed a method called Direct Cyclic Algorithm to obtain the stabilized response of a structure under cyclic periodically loading, which it is shown to be more efficient of an incremental analysis, much more increasing the size of the problem and the number of cycles required to obtain the steady cycle, but with a small error with respect to the incremental solution. The same cyclic analysis of the aggregate was performed by DCA, with the same time increment of the incremental analysis and the maximum number of Fourier terms allowed, that provides the steady cycle in 682 iterations. The convenience of this kind of analysis in terms of time is clear because the number of linear systems to be solved for a single cycle with the incremental method is surely higher that the number of linear system required for a single iteration of DCA. Moreover the full Newton method needs to assemble and decompose the stiffness matrix many times during a cycle, operation that dominates the entire computational effort, while the DCA uses always the elastic matrix. Each pass through the complete load cycle with DCA can therefore be thought of as few iterations of incremental method. In the simulation performed in this work, the DCA is about four time faster than the Incremental method. The DCA solution is close to the incremental solution but it is not exactly coincident. The error with the maximum number of Fourier terms on the strain difference, stress difference and energy of steady cycle is, however, a few percent for the softest grain and negligible for the average behaviour. The DCA is not able to deal with geometric non-linearities, but the error of neglecting the rotation of crystal axis during the deformation has tested and it was shown to be acceptable. We have also seen that a further reduction of the computational cost can be achieved by using a lower number of Fourier terms with a slightly larger error.

This result can be surely used in a future work where the model will take into account the real morphology, number and orientations of grains, using EBSD analysis of specimens used in experimental investigations. Since the size of this aggregates, the DCA is maybe the only possible way to compute numerical simulations, allowing a real comparison between mechanical fields at grain scale and to study the interaction between the local constitutive law of the single crystal and the macroscopic behaviour of the specimen, for example with respect to crack initiation and fatigue life. These informations could lead to the definition of new fatigue criteria.
**Acknowledgements** This work is the result of the internship of Domenico Magisano at the Laboratory of Mechanics of Lille, funded by an Erasmus Placement Grant of the Lifelong Learning Programme and by a Master Grant provided by the Lille 1 University.

## References

1. ABAQUS 6.10 (2010) Theory and user's manual. Dassault systemes

2. Cailletaud G, Doquet V, Pineau A (1991) Cyclic multiaxial behaviour of an austenitic stainless steel. In: Fatigue under biaxial and multiaxial loading 12

3. Dang Van K, Papadopoulos IV (1999) High-cycle metal fatigue: theory to applications. CISM courses and lectures 392. Springer, Vienna

4. Evard P, El Bartali A, Aubin V, Rey C, Degallaix S, Kondo D (2010) Influence of boundary conditions on bi-phased polycrystal microstructure calculation. Int J Solid Struct 47(16):1979–1986

5. Franciosi P (1984) Etude théorique et expérimentale du comportement elastoplastique des monocrystaux métalliques se déformant par glissement. PhD thesis

6. Frank FC, Read WT (1950) Multiplication processes for slow moving dislocation. Phys Rev 79(4):722–723

7. Guilhem Y (2011) Numerical investigation of local mechanical fields in 316L steel polycrystalline aggregates under fatigue loading. PhD thesis

8. Huang Y (1991) A user-material subroutine incorporating single crystal plasticity in the ABAQUS finite element program. Division of applied sciences, Harvard University, Cambridge, MA, pp 02138

9. Hutchinson JW (1976) Bounds and self-consistent estimates for creep of polycrystalline materials. Proc R Soc Lond A 348

10. Korsunsky AM, James KE, Daymond MR (2004) Intergranular stresses in polycrystalline fatigue: diffraction measurement and self-consistent modelling. Eng Fract Mech 71:805–812

11. Mandel J (1965) Generalization de la theorie de la plasticité de W. T. Koiter. Int J Solids Struct 1(3):273–295

12. Monnet G (2009) A crystalline plasticity law for austenitic stainless steel

13. Mu P (2011) Study of crack initiation in low-cycle fatigue of an austenitic stainless steel. PhD thesis

14. Papadopoulos IV (1994) A new criterion of fatigue strength for out-of phase bending and torsion of hard metals. Int J Fatigue 16(6):377–384

15. Pommier B (2003) Détermination de la réponse asymptotique d'une structure anélastique sous chargement thermomécanique cyclique. PhD thesis

16. Quey R, Dawson PR, Barbe F (2011) Large-scale 3D random polycrystals for the finite element method: generation, meshing and remeshing. Comput Methods Appl Mech Eng 200:1729–1745

17. Schmid E, Boas W (1950) Plasticity of crystals with special reference to metals. F. A. Hughes (London)

18. Seghir R, Charkaluk E, Dufrénoy P, Bodelot L (2010) Thermomechanical coupling in crystalline plasticity under fatigue loading. Proc Eng 2:1155–1164

19. Seghir R, Bodelot L, Charkaluk E, Dufrénoy P (2011) Numerical and experimental estimation of thermochemical fields heterogeneity at grain scale of 316L stainless steel. Comput Mater Sci 53:464–473

20. Seghir R (2012) Experimental and numerical investigation of thermomechanical couplings and energy balance in metallic polycrystals. PhD thesis, Ecole Centrale de Lille

21. Simonovski I, Cizelj L (2009) The influence of grain structure size on microstructurally short crack. J Eng Gas Turbines Power 113(4)

22. Socie DF, et Marquis GB (2000) Multiaxial fatigue. SAE Inc Warrendale