# Large plastic deformation of voids in crystals

Jalal Smiri, Joseph Paux, Oğuz Umut Salman, Ioan R Ionescu

► To cite this version:

Jalal Smiri, Joseph Paux, Oğuz Umut Salman, Ioan R Ionescu. Large plastic deformation of voids in crystals. International Journal of Solids and Structures, In press, 324, pp.113657. 10.1016/j.ijsolstr.2025.113657. hal-05390211

HAL Id: hal-05390211

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

Submitted on 2 Dec 2025

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
Large plastic deformation of voids in crystals

Jalal Smiri*†, Joseph Paux*, Oğuz Umut Salman* and Ioan R. Ionescu*‡

version : 2025-04-07

**Abstract**

The mechanisms of void growth and coalescence are key contributors to the ductile failure of crystalline materials. At the grain scale, single crystal plastic anisotropy induces large strain localization leading to complex shape evolutions. In this study, an Arbitrary Lagrangian-Eulerian (ALE) framework for 2D crystal plasticity combined with dynamic remeshing is used to study the 2D shape evolution of cylindrical voids in single crystals. The large deformation and shape evolution of the voids under two types of loading are considered: (i) radial and (ii) uni-axial loadings. In both cases, the voids undergo complex shape evolutions induced by the interactions between slip bands, lattice rotations and large strain phenomena.

In case (i), the onset of the deformation revealed the formation of a complex fractal network of slip bands around the voids. Then, large deformations unearth an unexpected evolution of the slip bands network associated with significant lattice rotations, leading to a final hexagonal shape for the void. In case (ii), we obtain shear bands with very large accumulated plastic strain (> 200%) compared to the macroscopic engineering strains (< 15%). A high dependence between crystalline orientations, slip band localization and therefore shape evolution was observed, concluding in a high dependency between crystalline orientation and void shape elongation, which is of prime importance regarding coalescence of the voids, thus to the formation of macro-cracks.

**keywords** : Void growth, Crystal plasticity, Slip bands, Single crystal, large deformation, Eulerian approach

<footnote>*LSPM, University Sorbonne-Paris-Nord, Villetaneuse, France</footnote>

<footnote>†Universite Marie et Louis Pasteur, CNRS, Institut FEMTO-ST, F-25000 Besançon, France</footnote>

<footnote>‡IMAR, Romanian Academy, Bucharest, Romania</footnote>
# 1 Introduction

The mechanisms of void growth and coalescence are key contributors to the ductile failure of crystalline materials, occurring through plastic flow around pre-existing voids or nucleated cavities at second-phase particles. These cavities can emerge on two distinct scales within polycrystalline material: at the grain scale (i.e., intragranular and intergranular porosities) and at the polycrystalline scale. The progression of void growth is highly dependent on the nature of the surrounding plastic deformation, which can be modeled by (i) crystal plasticity at the grain scale and (ii) isotropic or texture-induced anisotropic plasticity at the polycrystal scale. Thus, micromechanical modeling of void growth and coalescence must account for these varying local behaviors in order to adequately capture the influence of porosity on the overall macroscopic response. However, it is generally infeasible to derive comprehensive analytical models for the plasticity of ductile porous solids due to the complexity involved.

The aim of this paper is to carry out an analysis of large deformation voids growth in single crystals (HCP, FCC) using Eulerian crystal plasticity finite element simulations. This analysis is driven by previous work on the analytical homogenization of porous single crystals with internal cavities, as demonstrated in several models that address void growth [17, 24, 26, 30, 33, 34, 36, 45, 46]. In particular, it extends the work of Paux et al. (2022), who contribute to the understanding of ductile failure in crystalline materials by addressing the anisotropic plastic behavior of single crystals using an innovative approach that includes a new class of piecewise constant velocity fields derived from numerical FFT analysis restricted to small strain, which effectively improves the estimation of the macroscopic yield stress, particularly in models with periodic void distributions in single crystals.

While the model proposed in [33] provides a robust framework to estimate yield criteria, it lacks the ability to predict the evolution of void shapes during deformation, a critical factor in understanding the long-term behavior of materials under stress. To address this shortcoming, the present paper makes use of an Eulerian numerical model that not only captures the onset of yielding but also accurately predicts the evolution of void shapes. We aim to extend the contributions of [33] by providing a more comprehensive tool for the analysis of ductile failure in crystalline materials. The model is applied to the well-known experimental study by [7], which observed that voids in an HCP crystal evolve into a polyhedral cavity, a prediction also made in the same study. In that work, void growth was successfully modeled using a finite element approach that incorporated elastic-plastic behavior and kinematic hardening, with growth rates associated with pris-
matic slip activation. In this context, the proposed Eulerian rigid (visco) plastic proposed model offers an alternative or complementary approach for analyzing void evolution.

Conventionally, large deformation simulations are performed in a Lagrangian framework, where the mesh follows the deforming material. While effective for moderate strains, as the deformation progresses, mesh distortion introduces significant errors and ultimately leads to simulation failure. To address these challenges, this paper adopts an Eulerian approach, which is particularly effective for simulating extreme deformations. In the Eulerian formulation, a fixed grid allows materials to flow through without distorting the mesh, making it ideal for analyzing complex material behavior during severe plastic deformation.

The elastic-(visco-)plastic model needs a permanent interplay between both Eulerian (more appropriate for plasticity) and Lagrangian (more appropriate for elasticity) descriptions. However, in applications involving large deformations of metals, the elastic component of the deformation is small with respect to the inelastic one, and can be neglected by using a rigid-(visco)plastic approach (see for instance [22, 23, 25]). This simplified model takes important theoretical and numerical advantages by using only one (Eulerian) configuration. This work starts from the Eulerian framework in modeling large deformations in crystal plasticity introduced by [5] and uses the high-resolution numerical methods developed in [6].

Let us outline the contents of the paper. The aim of the first section is to recall from [5] the 3D rigid-viscoplastic crystal model and its 2D simplification. In the next two sections, the proposed model is evaluated under two distinct plane-strain configurations involving intergranular porous crystals : (1) Growth evolution under radial loading (hydrostatic tension) applied on a hexagonal close-packed (HCP) crystal and (2) Growth evolution under uniaxial loading (tensile test) applied to a face-centered cubic (FCC) crystal. These cases are chosen to illustrate the model’s capability in different crystallographic structures and loading conditions, thereby demonstrating its broader applicability. Finally, after the conclusion, an appendix recalls the numerical algorithm form [6] and details the re-meshing procedure used to obtain the results.

## 2 Eulerian rigid-plastic approach of crystal plasticity

Plastic deformation in crystalline materials is highly complex and involves multiple length scales. These spatial heterogeneities range from the atomistic scale (dislocation cores and grain boundary structures) to mesoscale dislocation patterns and grain microstructures,
extending to the macroscopic scale of the specimen. At the macroscopic level, plasticity manifests itself as a smooth flow described by a continuous stress-strain response.

Several theories have been developed to describe and understand plastic deformation at the macroscopic level. Continuum Crystal Plasticity (CP) theory is the most widely used approach for modeling crystal plasticity. This theory incorporates lattice-based kinematics into the classical continuum framework. Its original mathematical formulation was introduced by Hill [19] and Hill and Rice [20], with initial applications by Asaro and Rice [1, 3] and Pierce, Asaro, and Needleman [35]. Since then, many authors have further developed CP theory (see [39] for a comprehensive overview). The classical Continuum crystal theory of plasticity assumes that crystalline materials undergo irreversible flow when applied stresses exceed certain thresholds. This theory models the stress-strain response with continuous curves, implying that the heterogeneities are averaged out. This approach has been successful in reproducing key plasticity phenomena such as yielding, hardening, and shakedown.

CP theory spans macroscopic length scales, from micrometers (μm) to millimeters (mm), due to its coarser representation of plastic deformation compared to other methods such as molecular dynamics or discrete dislocation dynamics approaches. In recent CP formulations, deformations are considered finite, and the continuum description clearly distinguishes between reference and deformed configurations.

## 2.1 Multiplicative decomposition of the deformation gradient

Consider a single crystal at time $t = 0$, free of any surface tractions and body forces and let choose this configuration, say $K_0 \subset \mathbb{R}^d$ (here $d = 2, 3$ is the space dimension) as reference configuration of the crystal. Let $K = K(t) \subset \mathbb{R}^d$ denote the current configuration. The incorporation of lattice features is achieved through a multiplicative decomposition of the total deformation gradient $F$ into elastic and plastic components:

$$F = F^e P. \qquad (2.1)$$

This decomposition implies a two-stage deformation process. First, $P$ transforms the initial reference state $K_0$ to an intermediate state $\tilde{K}$, characterized by plastic deformation only with no change in volume. $P$ is called the (visco)plastic deformation with respect to the reference configuration of a material neighborhood of the material point $X$ at time $t$. Then, $F^e$ brings the body to the final configuration $K$ through elastic deformation and rigid lattice rotation, i.e. $F^e = RU^e$ where $R$ denotes the rotation of the crystal axes
with respect to its isoclinic orientation, and $U^e$ is the elastic right stretch tensor.

Following [3, 27], $P$ is assumed to leave the underlying lattice structure undeformed and unrotated, ensuring the uniqueness of the decomposition in (2.1). The unique feature of CP theory is its construction of the plastic component $P$ by constraining dislocation kinematics. Plastic flow evolves along pre-selected slip directions via volume-preserving shears, leaving the crystal lattice undistorted and stress-free [31].

Since in applications involving large deformations and high strain rates, the elastic component of the deformation is small with respect to the inelastic one, it can be neglected and a rigid-viscoplastic approach will be adopted (such a hypothesis is generally used e.g. Hutchinson [22], Lebensohn and Tome [25], Kok et al.[23], etc.). That means that we neglect the elastic lattice strain $U^e$ by supposing that $U^e \approx I$. This leads to the following decomposition for the deformation gradient $F$ (see for example Kok et al. [23]):

$$F = RP. \qquad (2.2)$$

Such a hypothesis is valid since during forming or other industrial processes, the elastic component of deformation is negligibly small (typically $10^{-3}$) in comparison to the plastic component (typically $>10^{-1}$). It is also to be noted that once the elastic/plastic transition is over, the stress evolution in the grains is controlled by plastic relaxation (see Tomé and Lebensohn [25]).

## 2.2 Eulerian description of the lattice rotations

Crystal slip systems are labeled by integers $s = 1, \dots, N$, with $N$ denoting the number of slip systems. Each slip system $s$ is specified by the unit vectors $(b_s^0, m_s^0)$, where $b_s^0$ is the slip direction and $m_s^0$ is the normal to the slip plane in the perfect undeformed lattice. Since the viscoplastic deformation does not produce distortion or rotation of the lattice, the lattice orientation is the same in the reference and intermediate configurations $K_0$ and $\tilde{K}$ and is specified by $(b_s^0, m_s^0)$, $s = 1 \dots N$.

Let's note $b_s = b_s(t)$ and $m_s = m_s(t)$ the glide direction and glide plane normal, respectively, in the deformed configuration $K$. At $t = 0$, it reads $b_s(0) = b_s^0$ and $m_s(0) = m_s^0$. Since elastic effects are neglected,

$$b_s = Rb_s^0, \quad m_s = Rm_s^0. \qquad (2.3)$$
Note that (2.3) implies that $b_s$ and $m_s$ are unit vectors. Furthermore,

$$b_s \otimes m_s = R (b_s^0 \otimes m_s^0) R^T. \qquad (2.4)$$

We seek to express the lattice evolution equations only in terms of vector and tensor fields associated with the current configuration. Let $v = v(t, x)$ be the Eulerian velocity field, $L$ the velocity gradient, $D$ the rate of deformation, and $W$ the spin tensor,

$$L = L(v) = \nabla v, \quad D = D(v) = (\nabla v)^{\text{symm}}, \quad W = W(v) = (\nabla v)^{\text{skew}}. \qquad (2.5)$$

The viscoplastic deformation is only due to slip ; the slip contribution to the viscoplastic deformation being ([38], [48]),

$$\dot{P}P^{-1} = \sum_{s=1}^{N} \dot{\gamma}^{s} b_{s}^{0} \otimes m_{s}^{0}, \qquad (2.6)$$

where $\dot{\gamma}^s = \dot{\gamma}^s(t)$ is the viscoplastic shear rate on the slip system $s$. If we denote by

$$M_s = (b_s \otimes m_s)^{\text{symm}}, \quad Q_s = (b_s \otimes m_s)^{\text{skew}}. \qquad (2.7)$$

then, using $L = \dot{F}F^{-1} = \dot{R}R^T + R\dot{P}P^{-1}R^T$, Eqs. (2.4) and (2.6), the rate of deformation $D$ can be written as

$$D = \sum_{s=1}^{N} \dot{\gamma}^{s} M_{s}. \qquad (2.8)$$

Taking the anti-symmetric part of $L$, we obtain that the spin tensor is $W = \dot{R}R^T + \sum_{s=1}^{N} \dot{\gamma}^s Q_s$ and a differential equation for the rotation tensor $R$:

$$\dot{R} = (W - \sum_{s=1}^{N} \dot{\gamma}^{s} Q_{s}) R. \qquad (2.9)$$

The evolution equations (2.9) describe the evolution of the lattice in terms of vector and tensor fields associated with the current configuration.

## 2.3 Plastic and visco-plastic flow rules

In order to complete the model, we need to provide the constitutive equation for the slip rate $\dot{\gamma}_s$ as a function of $\tau_s$, the stress component acting on the slip plane of normal $m_s$
in the slip direction $b_s$. In the current configuration, $\tau_s$ is expressed as

$$ \tau^s = \sigma : M_s, \qquad (2.10) $$

where $\sigma = \sigma(t)$ is the Cauchy stress tensor acting in the current configuration $\mathcal{K}$ while $M_s$ is defined by (2.7). Note that $\{\tau^s\}_{s=1,N}$ are not independent; they belong to a fifth dimensional space of $\mathbb{R}^N$ corresponding to the dimension of the space of deviatoric stresses.

To determine the shear strain rates $\dot{\gamma}_\alpha$ relative to the local stress, a constitutive law is needed. Various proposals exist, ranging from phenomenological to more physically based approaches. One simple phenomenological approach assumes that $\dot{\gamma}_\alpha$ depends on the stress only through the resolved shear stress $\tau^s$. The associated plastic dissipation functional is neither strongly convex nor differentiable, so the solution could not be unique. Additional assumptions are needed in order to restrict the number of solutions.

One way to overcome this difficulty of determining the active slip systems problem is to adopt a rate-dependent approach for the constitutive response of the single crystal. A widely used rate-dependent (viscoplastic) model is the Norton type model, which relates the shear strain rate $\dot{\gamma}^s$ on a slip system $s$ to the resolved shear stress $\tau^s$ through a power-law (see Asaro and Needleman [2])

$$ \dot{\gamma}^s = \dot{\gamma}_0^s \left| \frac{\tau^s}{\tau_c^s} \right|^n \mathrm{sign}(\tau^s), \qquad (2.11) $$

where $\dot{\gamma}_0$ is a reference shear strain rate, while the exponent $n$ has a fixed value. The Schmid model is recovered for large values of $n$ ($n \to \infty$).

Another regularization of the Schmid law, which will be used here in the numerical computations, can be done by using a Perzyna-like viscoplastic law of the form:

$$ \dot{\gamma}^s = \frac{1}{\eta_s} [|\tau^s| - \tau_c^s]_+ \mathrm{sign}(\tau^s), \qquad (2.12) $$

where $\eta_s$ is the viscosity, which may depend on the slip rate, and $[x]_+ = (x + |x|)/2$ denotes the positive part of any real number $x$. The Schmid model is recovered for vanishing viscosities $\eta_s$ (i.e. $\eta_s \to 0$).

## 2.4 2-D model with 3 slip systems

As a first approach to consider the shape evolution of voids in single crystal under large strain, a 2-D model ($d = 2$) is considered. This simplified model (by comparison with
3-D models) provides the possibility to investigate closely the formation of slip bands and the evolution of the lattice rotation around the void with reasonable computation time, therefore enabling an interpretation of the link between crystal plasticity features and shape evolution of the voids.

We will introduce here a two dimensional model with $N = 3$ slip systems. In this case let us denote with $\phi$ and $-\phi$ the angles formed by slip system $r = 1$ with the other two systems $r = 2, 3$ and let $\theta$ be the angle formed by the slip system 1 with the $Ox_1$ axis (see Figure 1). The three composite in-plane slip systems $b_1, b_2, b_3$ are specified by the angles $\theta, \theta + \phi, \theta - \phi$.

<figure><img src="image_2.png" /><figcaption>Figure 1: Two dimensional model with three slip systems.</figcaption></figure>

The main simplification for the 2-D problem comes from the lattice rotation $\mathbf{R}$ which is now a rotation $\mathbf{R}(\theta, \mathbf{e}_3)$ with angle $\theta$ along $Ox_3$ axis and we have

$$Q_r = \frac{1}{2}((1,0) \otimes (0,1) - (0,1) \otimes (1,0)). \qquad (2.13)$$

for all $r = 1, 2, 3$. (2.9) simplifies into a scalar form,

$$\dot{\theta} = \frac{\partial \theta}{\partial t} + \mathbf{v} \cdot \nabla \theta = \frac{1}{2} \left( \sum_{r=1}^{3} \dot{\gamma}^{r} - \left( \frac{\partial v_1}{\partial x_2} - \frac{\partial v_2}{\partial x_1} \right) \right). \qquad (2.14)$$

As already stated in the literature, this model is physically sound in two situations : i) the in-plane deformation of a FCC crystal [12], and ii) the slip in the basal plane of a hexagonal crystal [11].
## In-plane deformation of a FCC crystals

Rice [38] showed that certain pairs of the three-dimensional systems that are potentially active need to combine in order to achieve plane-strain deformation. For an FCC crystal, with 12 potentially active slip systems, we consider $Ox_3$ axis to be parallel to [110] in the crystal basis, which means that the plane-strain plane ($Ox_1x_2$) is the plane [110] - [001]. Some geometrical constraints (see [6]) have to be satisfied such that $N = 3$ pairs of composite systems will give deformation in the plane ($Ox_1x_2$). For all $r = 1, 2, 3$ we have denoted by $b_r$ and $m_r$ the normalized projections of the corresponding three-dimensional slip directions and normal directions ($k, l$) onto the $x_1x_2$-plane. The angle $\phi$ between the slip systems 1 and 2 is given by

$$ \phi = \arctan(\sqrt{2}) \approx 54.7^{\circ}. \qquad (2.15) $$

Since there are some scalar factors between the first two components of the in-plane systems and the 3-dimensional ones given by $q_1 = \frac{1}{\sqrt{3}}$, $q_2 = q_3 = \frac{\sqrt{3}}{2}$, the 2-D composite slipping rate $\dot{\gamma}^r$ corresponds to the 3-D slipping rate $\dot{\gamma}^k$ multiplied by $2q_r$. As it follows from [6], the 2-D yield limit $\tau_0^r$ corresponds to the 3-D yield limit $\tau_0^k$ divided by $q_r$ while the 2-D viscosity $\eta^r$ corresponds to the 3-D yield limit $\eta^k$ divided by $q_r^2$.

## Slip in the basal plane of a hexagonal crystal

Alternatively, this situation corresponds to hexagonal close-packed (HCP) crystals (such as Ti, Mg, Zr, etc) under plane strain, with the plane of deformation aligned with the basal plane (0001) of the hexagonal lattice. This situation has been considered experimentally in [7]. In this situation, the plastic strain is mainly accommodated by the three prismatic slip systems, i.e. the $(10\bar{1}0)\langle 1\bar{2}10 \rangle$ slip family. Each of those slip systems can act independently from the two others, as its strain is a plane strain. One might remark that combinations of other slip systems of HCP crystals from the basal and the pyramidal slip families could also lead to plane strain, but it would require significantly more energy, and therefore never occur. Thus, those slip families are not considered in this work. The prismatic slip systems are symmetric, leading to

$$ \phi = \pi/3 = 60^{\circ}. \qquad (2.16) $$
## 2.5 Initial and boundary value problem formulation

We begin by presenting the equations governing the motion in a domain $D = D(t)$ of an incompressible rigid-viscoplastic crystal for the simplified 2D model with $N = 3$ slip systems. In an Eulerian description of a crystal visco-plasticity theory, the unknowns are: the velocity $v: [0, T] \times D \to \mathbb{R}^2$, the crystal lattice orientation, i.e. the rotation $\theta: [0, T] \times \mathbb{R}$ and the Cauchy stress $\sigma: [0, T] \times D \to \mathbb{R}^{2 \times 2}$. Let $\sigma = \sigma' + pI$, where $\sigma'$ is the stress deviator while $p: [0, T] \times D \to \mathbb{R}$ is the pressure.

The momentum balance in the Eulerian coordinates reads

$$ \rho^{\text{mass}}(\partial_t v + v \cdot \nabla v) - \text{div} \sigma' + \nabla p = \rho^{\text{mass}} f \quad \text{in } D, \qquad (2.17) $$

where the mass density $\rho^{\text{mass}} > 0$ and the body forces $f$ are supposed to be known. The incompressibility condition is

$$ \text{div} v = 0 \quad \text{in } D. \qquad (2.18) $$

The momentum balance equations are completed by the constitutive equation, which relates the stress tensor $\sigma$ and the rate of deformation tensor $D(v)$ (see (2.8)) through the evolution equations for each slip system $s$ given by the viscous regularizations (2.11) or (2.12).

The boundary $\partial D$ of the domain $D$ is decomposed into two disjoint parts, $\Gamma_v$ and $\Gamma_s$, such that the velocity is prescribed on $\Gamma_v$ and traction is prescribed on $\Gamma_s$, at any time $t$:

$$ v(t) = V(t) \quad \text{on} \quad \Gamma_v, \quad \sigma(t)n = S(t) \quad \text{on} \quad \Gamma_s, \qquad (2.19) $$

where $n$ stands for the outward unit normal on $\partial D$, $V$ is the imposed velocity and $S$ is the prescribed stress vector.

We also consider another partition of $\partial D$ into $\partial_{in}D(t)$ and $\partial_{out}D(t)$ corresponding to incoming $(v \cdot n < 0)$ and outcoming $(v \cdot n \ge 0)$ flux. The boundary conditions associated to the lattice evolution equations (2.14) reads

$$ \theta(t) = \theta^{\text{in}}(t), \quad \text{on} \quad \partial_{\text{in}}D(t), \qquad (2.20) $$

To the field equations, we add the initial conditions

$$ v(0) = v^0, \quad \theta(0) = \theta^0, \quad \text{in} \quad D. \qquad (2.21) $$

where $v^0$ is the initial velocity and $\theta^0$ gives the initial orientation of the crystal lattice.
Note that in this model, there is no need to prescribe initial conditions for the stress.
This is very convenient since the initial stress field is generally not known or it cannot be
easily measured.

## 2.6 Numerical strategy

To integrate the governing equations, a mixed finite-element and Galerkin discontinuous
strategy, developed in [5], will be used. We will give here only the principles of the
method used in the next sections while a detailed description is provided in Appendix
A. Time implicit (backward) Euler scheme is used for the field equations, which gives
a set of nonlinear equations for the velocities $v$ and lattice orientation $R$. At each
iteration in time, an iterative algorithm is used in order to solve these nonlinear equations.
Specifically, the variational formulation for the velocity field is discretized using the finite
element method, while a finite volume method with an upwind choice of the flux is
adopted for solving the hyperbolic equations that describe the evolution of the lattice
orientation. One of the advantages of the proposed numerical strategy is that it does
not require the consideration of elastic deformation since it is not based on an elastic
predictor/plastic corrector scheme (e.g. [42]).

It is to be noted that in the case of the proposed rigid-viscoplastic model (2.12),
numerical difficulties arise from the non-differentiability of the viscoplastic terms. That
means that one cannot use finite element techniques developed for Navier-Stokes fluids
(see for instance [14, 47]). To overcome these difficulties the iterative decomposition-
coordination formulation coupled with the augmented Lagrangian method of [10, 13] was
modified. The reason for this modification is that in the used crystal model, there is not
co-axiality between the stress deviator and the rate of deformation as it is the case in
the Bingham model used by [10, 13]. Note that this iterative decomposition-coordination
algorithm permits to solve at each iteration the equations for the unit vectors that define
the lattice orientation.

Furthermore, if the domain $\mathcal{D}$ occupied by the single crystal (or poly-crystals) varies
in time, then an arbitrary Eulerian-Lagrangian (ALE) description was adopted. In the
numerical simulations presented in the next sections we use the followings spatial dis-
cretization: P2 for the velocity field, P1 for the pressure, P1 discontinuous for the stress
field, slip rates, dislocation densities and lattice orientations.
<figure><img src="image_3.png" /><figcaption>Figure 2: Schematic representation (not to scale) of a porous single-crystal disc unit cell with a single circular void at its center.</figcaption></figure>

# 3 Void growth evolution under radial loading

## 3.1 Setting of the numerical experiment

In the present section, we present the numerical setting, including geometrical and material characteristics, initial and boundary conditions, followed by the results and discussion of the 2D void growth simulations under plane strain in an HCP single crystal.

Fig. 2 shows a schematic representation of the initial state of the one-void unit cell represented by a disc of radius $R_0$, containing an initially circular void at its center with a radius $r_0 = R_0/30$. This configuration results in an initial void fraction of $f_0 = \left(\frac{r_0}{R_0}\right)^2 = 0.111\%$.

On the external boundary, denoted by $\Gamma_0(t)$, we impose an expansion velocity $V_0$, corresponding to the in-compressible Eulerian velocity field $v^* = \frac{V_0 R_0}{r} e_r$, i.e.

$$ \mathbf{V}(t)(x, y) = \frac{V_0 R_0}{x^2 + y^2} (xe_1 + ye_2), \quad \text{for } (x, y) \in \Gamma_0(t), \qquad (3.1) $$

while the internal boundary, denoted by $\Gamma_1(t)$, is traction free (i.e. $S(t) = 0$). The initial conditions in velocity and in the crystal orientation where chosen to be

$$ \mathbf{v}^0(x, y) = \frac{V_0 R_0}{x^2 + y^2} (xe_1 + ye_2), \quad \theta^0(x, y) = \theta_0. \qquad (3.2) $$
The expansion velocity $V_0$ was chosen to be small such that the process is expected to be quasi-static. By integration on the Eulerian coordinate $r$ gives the external radius can be deduced to be $R^2(t) = R_0^2 + 2V_0R_0t$, which means that the engineering area deformation is given by

$$ \epsilon^{eng}(t) = \frac{R(t)^2 - R_0^2}{R_0^2} = \frac{2V_0t}{R_0}. \qquad (3.3) $$

The time interval $[0, T]$, the velocity $V_0$ and the rayon $R_0$ where chosen such that the final engineering volume/area deformation is $\epsilon_{final}^{eng} = 1.32\%$. Note that, even if the final $\epsilon_{final}^{eng}$ is small, the strains around the cavity are expected to be much larger. Based on the engineering area deformation formula above, we estimate it to be of order of $\epsilon_{final}^{eng}R_0/r_0 = 30\epsilon_{final}^{eng} \approx 40\%$. The material parameters used in this quasi-static computation is not directly specific to a particular material but reflect the order of magnitude typically associated with hexagonal close-packed (HCP) crystals. The material parameters are: density $\rho = 3200 \text{ kg/m}^3$, slip resistance $\tau_c = 20 \text{ MPa}$ for all systems with no hardening, and a viscosity $\eta$ chosen to be as small as possible (of order of $1\%$ of $\tau_c R_0 / V_0$) for numerical reasons. The initial orientation of the crystal was chosen to be $\theta_0 = 0^\circ$.

In the following subsections, we examine the evolution of slip bands around the cavity, from their initiation to large deformations under hydrostatic tension loading as predicted by the proposed numerical model. Given the distinct phases observed in the deformation fields, the analysis is divided into two stages, each discussed in a separate subsection. The first stage focuses on small deformations, analyzing the onset of shear and kink bands and comparing the findings with existing models from the literature. The second stage explores the large deformation regime, following the evolution of the deformation patterns, described in the onset phase.

## 3.2 Onset of the shear/kink band deformation

At low deformation levels, although the size and morphology of the void show minor variations, slip lines are present, allowing for their characteristics to be analyzed. In Figure 3, we present a qualitative comparison of the equivalent von Mises strain distribution in a periodic hexagonal unit cell with aligned crystal orientation $\theta = 0^\circ$, computed using the FFT method described in [33] with the deformation rate predicted by the current model. This comparison is particularly valid at low deformation levels where the two models are valid and where the localization patterns of deformation and deformation rate are similar. Despite the difference in cell micro-structure, hexagonal in [33] and disk-shaped in the current model, this comparison remains justified at equivalent low porosity levels (0.111%
and 0.24%, respectively), where edge effects are minimal and do not significantly impact the results around the cavity. All results are presented with an enlarged view around the void of the circular cell to clearly visualize the different fields. The zoomed region is shown as a rectangular area, with dimensions set to 26$r_0$ in width and 19$r_0$ in height, centered on the void, where $r_0$ is the initial radius of the void.

At the onset of strain localization, a fractal network emerges as observed in the spatial distribution of strain and strain rate (Figure 3). A (redundant) base pattern of twenty four slip lines, composed of two superimposed six-pointed stars with a 45-degree phase shift between them, recursively repeats himself with self-similarity from the exterior of the cell up to the void, forming successive layers. The base pattern is made up of two types of lines: lines parallel to slip directions, known as classical slip bands, and lines parallel to the slip plane normals, referred to as kink bands.

A slip band is defined as a region characterized by high plastic strain and low crystallographic reorientation, where a single slip system is predominantly active and the slip direction is parallel to the band. In contrast, a kink band is a region with high plastic strain and sharp crystallographic reorientation, also dominated by one slip system, but with the Burgers vector perpendicular to the band direction [9, 49].

The consistent appearance of this twenty-four-slip-line pattern in both models suggests that it is an intrinsic feature of the strain field in porous single crystals. Notably, there is a clear agreement between the two models in their description of the emergence of slip and kink bands.

<figure><img src="image_4.png" /><figcaption>Figure 3: Onset of the deformation: qualitative comparison using two different approaches. Left: equivalent strain obtained under small strain hypothesis, as reported in [33]. Right: zoom of the equivalent strain rate at a low deformation level ($\epsilon_{vol}^{eng} = 0.0033\%$) obtained with the Eulerian model.</figcaption></figure>

In this context of incipient plasticity, the three prismatic slip systems demonstrate
<figure><img src="image_5.png" /><figcaption>Figure 4: Onset of the deformation. The three slip systems of HCP crystal based on finite element simulation results : (left) $s = 1$, (middle) $s = 2$, and (right) $s = 3$. Top: slip deformation rates computed with the presented Eulerian model. Bottom: slips deformation computed with a strain gradient crystal plasticity model without hardening from [4].</figcaption></figure>

equivalent amounts of glide, as shown in Fig. 4 (Top), with all slip systems being simultaneously active. Each system is active in specific angular areas, forming, for each system, an orthogonal cross. The crosses of the different systems are equivalent considering the crystal symmetries. The system's activities maxed out at the contact with the void. A similar finding is reported by [4] (see Figure 4 bottom), where it is noted that slip contours are periodic with a period of 90°. Specifically, slip system one is predominantly active in the region $0 \le \theta \le 30°$, slip system two in $30 \le \theta \le 60°$, and slip system three in $60 \le \theta \le 90°$ of the upper right quadrant.

Having explored the initiation and early evolution of slip and kink bands at small deformation levels and seen an agreement between the present model and the model of Paux et. al. [33], we now shift our focus to the large deformation regime. In this stage, the behavior of the slip bands undergoes significant changes as deformation progresses, with noticeable effects on the microstructure and deformation patterns. In the following section, the response of the model is used to examine how these patterns evolve under increasing strain, highlighting the transformation of slip lines and the eventual disappearance of shear bands and/or kink bands.
<figure><img src="image_6.png" /><figcaption>Figure 5: Evolution of equivalent deformation strain rate and void morphology at different stages of the deformation: $\epsilon^{eng} = 0.2475\%$, $\epsilon^{eng} = 0.66\%$, $\epsilon^{eng} = 0.99\%$ and $\epsilon^{eng} = 1.32\%$.</figcaption></figure>

## 3.3 Large deformations towards a hexagonal shape

Figure 5 illustrates the evolution of the deformation rate field within the cell and the corresponding changes in void shape at large local deformation for various intervals. As deformation progresses, the void shape loses its isotropy, slip bands progressively become dominant in the very surrounding of the void, while the kink bands progressively lose activity and disappear from this vicinity. It is interesting to note that the second stage of the slip line network from the void undergoes opposite evolution and becomes mainly made of kink bands. The persistent activity of the slip bands at elevated local deformation levels drives the transformation of the initially circular void into a hexagonal shape, with its facets aligned along the slip directions. As deformation advances further, as shown from (c) to (d) in Figure 5, the deformation rate pattern, and therefore the (quasi) hexagonal void undergoes a relatively uniform expansion. The shear band activity tends to become less intense compared to the period (a) and (b), which can be explained by
<figure><img src="image_7.png" /><figcaption>Figure 6: Qualitative comparison between the current model finite element computations and the experimental observations. (Left) Void morphology at $\epsilon^{eng} = 0.66\%$. (Right) SEM image of the hexagonal cross-section of the cavity observed on the sample surface, showing the three prismatic slip systems (1, 2, 3) as identified in [7].</figcaption></figure>

the increase of the void size, as the strain rate is expected to be roughly proportional to $1/r^2$ [16, 33].

The final state corresponding to the deformation rate contour derived from Figure 5 (d) is again presented in Figure 6 (a), where we compare the final void shape from finite element simulation results (with its deformation rate) with an experimental observation of a void after plastic deformation reported by [7]. From a morphological perspective, both images exhibit a similar hexagonal shape of the void. The slip traces that are present in the experimental results suggest that the hexagonal shape results from slip bands around the void similar to those predicted by our numerical simulation.

Although both the experimental and numerical model consider an HCP crystal with three prismatic systems, it is important to recognize that the loading conditions affecting the void may not be identical. The experimental method employs poly-crystal specimens, making it challenging to determine the exact loading conditions applied in the vicinity of the void and complicates the direct comparison between the simulation and the experiment. Nonetheless, numerical simulations show a strong overall correlation with the experimental results. In particular, it was able to predict the disappearance of the kink bands around the voids under large strain, leading to the same slip activity around the void as predicted by [4] without relying on strain gradient plasticity. It is interesting to note that this result cannot be predicted by small strain simulations, where slip and kink bands are strictly equivalent. As the main difference between kink and shear band under
<figure><img src="image_8.png" /><figcaption>Figure 7: Variations in circular cell microstructure orientation θ (expressed in degrees) during hydrostatic tension loading in the basal plane of an HCP crystal, with the initial crystal orientation set to θ = 0° at : (a) ϵ<sub>eng</sub> = 0.125%, (b) ϵ<sub>eng</sub> = 0.2475%, (c) ϵ<sub>eng</sub> = 0.33%, (d) ϵ<sub>eng</sub> = 0.66%, (e) ϵ<sub>eng</sub> = 0.99% and (f) ϵ<sub>eng</sub> = 1.32%. The color scale of the contour plot is customized.</figcaption></figure>

large strain is the induced lattice rotation, the evolution of the crystalline orientation in
the domain is further investigated.

**3.4 Rotation of crystal lattice**

As noted earlier, at large deformations, the slip lines around the voids evolve into slip
bands, while the kink bands, which are present around the voids at smaller deformations,
cease to exist in the very vicinity of the void. One might remark that, at large deforma-
tion, the second star-like pattern from the void seems to follow the opposite evolution
and is mainly made of kink bands (Figure 5). Since kink bands play a key role in lattice
rotation mechanisms, this transition prompts an investigation into how crystal orientation is affected. The evolution of the lattice rotation around the void, presented in Figure 7, shows that, in regions far from the void, the crystal orientation remains unchanged (close to 0 degrees). We remark that the orientation exhibits a discontinuity at the angles $\varphi = \pi/6 + k\pi/3$ as the attractor does. However the discontinuity is more present near the void, where the deformations are much larger, than near the external boundary. This is in concordance with the stability analysis of [44]: the distance between the lattice orientation and its attractor decreases significantly for large strains (more than 50%) which are present near the void but are too small on the external boundary. Moreover, near the void, the orientation distribution formed a star-like pattern with six identical branches, matching the six-star pattern associated with kink bands, which is expected. This pattern is consistent with results at low deformation provided by [12]. Within each branch, the orientation varies from $-10^\circ$ to $10^\circ$, with one half of the branch rotating positively and the other half negatively. This bidirectional rotation seems to create a repulsive effect that contributes to the kink bands fading. Then, the majority of the accumulated plastic strain and slip activity are concentrated in the slip bands. The corner defining the angles of the hexagon forms at the points of contact between the slip bands and the void. In these corner regions, the single crystal experiences a minimal rotation, which allows the stability of the slip line pattern and therefore the stability of the hexagonal morphology.

As we move closer to the cavity, the highest rotation level of $30^\circ$ degrees appears at the center of each edge of the hexagon, represented by a small half branch motif mirrored on the opposite side. These points seem to correspond to the location where the kink bands used to touch the void. One half branch shows negative rotation, while the other exhibits positive rotation. The disappearance of kink bands at high deformation levels implies that the orientation contour plot shows no significant differences between small and large deformations. The only notable change is the emergence of a second star-like pattern in the final simulation stage, slightly farther from the first, which is consistent with the presence of kink bands in the second stage of the slip line network from the void.

# 4 Void growth evolution under uniaxial loading

## 4.1 Setting of the numerical experiment

For FCC crystal under uniaxial loading, the one-void domain is represented by a rectangular pillar of dimension $H_0 \times L_0$ with $H_0 = 2L_0$, containing an initially circular void at the center of radius $r_0 = \frac{L_0}{30}$ (see Figure 8). This configuration results in an initial void
<figure><img src="image_9.png" /><figcaption>Figure 8: Schematic diagram of a tensile test with four rigid zones. The considered initial crystalline orientations are illustrated on the left. Two diagonal slip bands forms of angle $\phi$, indicating the expected localized plastic deformation.</figcaption></figure>

fraction of $f_0 = \frac{\pi r_0^2}{H_0 L_0} = 0.1745\%$.

The pillar's top is animated with a vertical velocity $V = Ve_2$ (small enough to deal with a quasi-static loading), while the bottom part is kept fixed ($V = 0$). The lateral boundary $\Gamma_c(t)$ is stress free ($S = 0$). The FCC single crystal has initially the orientation $\theta_0$, i.e.

$$ \theta(0, x, y) = \theta_0, $$

with three choices of $\theta_0 \in \{0^\circ, \frac{\phi}{2} = \frac{54.7^\circ}{2}, 65^\circ\}$ to gain an understanding of the role played by the angle between the loading axis and the crystal lattice.

The engineering strain is given by

$$ \epsilon^{eng}(t) = \frac{H(t) - H_0}{H_0}, $$

where $H(t)$ is the current height of the pillar. The material parameters used for this simulation are: density $\rho^{mass} = 16650 \text{ kg/m}^3$, slip resistance for all systems equal to $\tau_c = 220 \text{ MPa}$ with no hardening, and a viscosity $\eta$ chosen, for numerical reasons, to be as small as possible to ensure that the system is rate-independent (around 1% of $\tau_c H_0 / V$). The upper velocity $V$ and the time interval $[0, T]$ where chosen such that the
final engineering strain is $\epsilon_{final}^{eng} = 15\%$.

In this test, the top surface of the specimen remains fixed horizontally, ensuring that both the top and bottom surfaces stay parallel and retain the same horizontal coordinates at their extremities, despite significant deformation. The analysis focuses on how shear and kink bands influence the evolution of the void, which is affected by the initial crystallographic orientation of the FCC crystal.

In this study, we analyzed three cases that represent three different initial crystallographic orientations of the:

* The first case corresponds to an orientation of $\theta_0 = 0°$, aligning with one of the three initial slip systems, namely the slip system 1.

* The second case involves an orientation of $\theta_0 = 54.7°/2$, which lies between two slip systems, one at $0°$ and the other at $\theta_0 = 54.7°$.

* The third case considers an orientation of $\theta_0 = 65°$, which intends to represent an unparticular orientation.

## 4.2 Evolution of the cavity shape and shear bands localization

The Figures 9, 10 and 11 respectively present the accumulated equivalent strain, the accumulated slip activity on each system, and the lattice rotation in the three cases at different stages of the tensile test. A first general remark is that the deformation is extremly localized in shear bands, resulting locally in more than 200% accumulated strain, while the global strain is less than 15%. This was expected since perfect plasticity (without hardening) is considered. The first case ($\theta_0 = 0°$, top panels of Figure 9) shows the classical two crossed shear bands at $\phi = 45°$ and $-45°$ (as presented in Figure 8), while the two other cases ($\theta_0 = 54.7°/2$ and $65°$, middle and bottom panels of Figure 9) are mostly driven by a unique shear band, showing a high influence of the initial crystalline orientation on the strain field.

Initial observations show that all voids are distorted under large strain, but that, depending on the initial crystalline orientation, the void either tends to collapse or to grow. This is consistent with the numerical results of [15], which observed a shape evolution dependency in both crystallographic orientation and stress triaxiality. At high triaxialities, voids tend to grow and coalesce, while at lower triaxialities, distortion, localized shearing, and collapse are more prevalent. Although the effect of stress triaxiality is not analyzed in detail here, it is important to note that the low triaxiality in the uniaxial
<figure><img src="image_10.png" /><figcaption>Figure 9: Accumulated plastic strain $\epsilon_a^p$ ($\dot{\epsilon}_a^p = |D|$) evolution of a FCC pillar under tensile loading with three initial crystallographic orientations $\theta_0 = 0^\circ$ (top), $\theta_0 = 54.7^\circ/2$ (middle) and $\theta_0 = 65^\circ$ (bottom) at four levels of engineering strains $\epsilon^{eng} = 0, 0.05, 0.1$ and $0.15$.</figcaption></figure>
tensile test favors void distortion, collapse, or closure, while, in the radial growth case discussed earlier, higher triaxiality leads to void growth rather than collapse.

Continuing from Fig. 9, the scenario of void collapse is seen in the middle panels of the figure, which correspond to an initial angle of $\theta_0 = 54.7^\circ/2$. Due to the single localized shear band, the circular void progressively flattens and becomes more elongated with continued plastic deformation, eventually forming an almost linear ellipse, with a significantly reduced height (high eccentricity). This scenario is consistent with experimental results obtained with polycrystalline material by Nemcko et al. [32], where similar shear band localization is induced by anisotropic voids distribution.

In contrast to the flattening behavior, the void growth, mainly seen in the top panels of the figure ($\theta_0 = 0^\circ$), occurs through two crossed shear bands, stretching the void along two directions, preventing one dimension from becoming disproportionately small relative to the other. The double shearing is also observed in the bottom panels ($\theta_0 = 65^\circ$), but with a dominant shear band, leading to an anisotropic growth of the void, resulting in an elongated void shape, but without collapsing.

These first observations suggest that void shape evolution under uniaxial loading is mainly driven by the presence or absence of different shear bands around the void. Let's analyse in detail the observed shear bands.

First let us remark that in a slip band in traction the velocity gradient is $L = 2dc_1 \otimes c_2$ where $c_1$ and $c_2$ are the directions associated to the slip band (see Figure 8) and $d > 0$ for $\varphi > 0$ and $d < 0$ for $\varphi < 0$. From (2.5) we get that $D = d(c_1 \otimes c_2 + c_2 \otimes c_1)$ and $W = d(c_1 \otimes c_2 - c_2 \otimes c_1)$, which means that the rotations are in the clock sense for $\varphi > 0$ and in the trigonometric sense for $\varphi < 0$. The difference between the bands with positive or negative orientation comes from the induced rotations, which therefore are expected to have a huge influence on the evolution of the slip activity. We expect that the lattice in the band will rotate, in the sense given by $W$ (i.e. by $\varphi$), till one slip direction is parallel to the slip band, such that, after a while, only one slip system is active. Let us notice that this scenario is valid for one straight shear band with a given constant (in time) orientation $c_1$. This is not always the case, as a complicate process in which the sample shape, the void position, the boundary conditions and the crystalline orientation may interplay to explain the existence of one or several slip bands and their orientations.

For the case $\theta_0 = 0^\circ$ from Figure 9, we deal at the beginning with two and at the end with four shear bands. The angles $\varphi$ of all these bands are close to $\pi/4$ or $-\pi/4$. For $\varphi = \pi/4$, the rotation is in the clock sense, hence the first system parallel with the slip band is $s = 2$. This can be seen in Figure 11 top, with a negative lattice orientation in
the slip band. Moreover, we remark in Figure 10 top that the single active slip system is $s = 2$. For $\varphi = -\pi/4$, the crystal rotate in the trigonometric sense ($\theta$ positive in the band) and the active slip system is $s = 3$ (see Figures 10 and 11 top). At the beginning of the process the presence of the void gives the localization of the two slip bands in the center of the "X", but the void enlarges a lot during the deformation. This larger void will let another slip band system appear parallel with one of the incipient slip bands, but at the end a fourth slip band occurs and the shape of the void became symmetric.

The case $\theta_0 = 54.7^\circ/2$ does not follow the schematic representation with four rigid blocks depicted in Figure 8. The deformation is mostly driven by an unique shear band, with an orientation varying from $\varphi = -60^\circ$ to $\varphi = -70^\circ$, approximately. Simultaneously, significant body rotations occur, as can be seen by the rotations of the left and right boundaries of the domain, and confirmed by lattice rotation observation in the middle panels of Figure 11. The combined influence of shear band deformation and rigid body rotation leads to a complex evolution, hardly explainable by the kind of analysis provided for the $0^\circ$ case, as it cannot be reduced to a simple band deformation. Moreover, one can observe in Figure 9 that the shear band and the void follow the rigid body rotation, and therefore the slip band does not remain orientated at $\varphi = -60^\circ$ during the evolution of the considered geometry. This competition between shear band deformation and rigid body motion prevents from the stabilization of the lattice orientation on a slip band, leading to multi-slip deformation during all the tensile test, with mostly slip activity from the two systems 2 and 3 (Figure 10, middle panels). One might remark that the rigid body rotation might be a consequence of the imposed null lateral displacement at the top and at the bottom of the pillar, making a unique shear band unable to perform the necessary deformation alone, as it induces non null lateral displacement, as schematized in Figure 8 (velocity vectors $V_l$ or $V_r$). Different boundary conditions may lead to different shear band localization.

Finally, a configuration closed to the $54.7^\circ/2$ case is obtained for $\theta_0 = 65^\circ$, with a dominant shear band initially oriented at approximately $\varphi = 50^\circ$ and significant rigid body rotation. Again, this leads to multi-slip activity, with dominant activity of the systems 2 and 3 (bottom panels of Figure 10), and no stabilization of the lattice rotation (bottom panels of Figure 11). Contrary to the case $\theta_0 = 54.7^\circ/2$, a second shear band with lower intensity appears at approximately $\varphi = -50^\circ$. These combined shear band
deformations lead to a scattering of the dominant shear band, whose past plastic deformation is transported by the second shear band deformation, leading to the presence of layers of past shear band deformations (Figure 9, bottom panel d). This scattering looks like the trace of plastic slips experimentally observed by Crépin et al. [7] (see Figure 6), and suggests sequential activation of the two shear bands.

# 5 Conclusions

In this paper, the large deformation growth of voids in single crystals (HCP and FCC) was investigated. To this end, an Eulerian crystal plasticity 2D model with three slip systems was employed and numerically solved in finite element simulations. The voids' deformation was analyzed under two types of loading: (i) radial and (ii) uni-axial. In both cases, the voids, which have an initial circular shape, undergoes significant shape distortion, and shear bands appear, where the accumulated plastic strain is very large (> 200%), even for moderate global engineering strains (< 15%).

In case (i), the onset of deformation shows the formation of a slip and kink bands network (already observed in [33]). This network evolves under large strain, with vanishing of the kink bands around the voids, leading only the slip bands drive the evolution of the shape of the void. It ultimately leads to an *hexagonal shape* of the void, consistently with experimental observations (see for instance [7]). In the meantime, significant lattice rotations are observed around the void, in accordance with the kink band pattern, explaining the loss of symmetry between kink and slip bands in this vicinity. The final orientation of the crystal around the void is given with 6 lines of discontinuity, with rotation up to 30°, while, far from the void, accumulated deformations are not large enough to induce significant lattice rotations.

In case (ii) we have considered three initial orientations of the crystal with respect to the loading axis. We get two different scenarios of the void evolution. Only one scenario corresponds to the classical one involving four blocks separated by slip bands with only one active slip system in each band. In this scenario the initial circular void has a final rhomboidal shape. The other one exhibits only one shear band with a variable (in time) orientation, and with two active slip systems. The final shape of the void is a thin rectangle oriented almost parallel with the shear band. In this scenario we expect that the void will collapse in a crack and an auto-contact will occur.

These results demonstrate the ability of the Eulerian approach to model crystal plasticity under very large strain, with shear band localization and significant crystalline
<figure><img src="image_11.png" /><figcaption>Figure 10: The spatial distribution of the accumulate slip $\gamma_a^s$ ($\dot{\gamma}_a^s = |\dot{\gamma}^s|$) on each system at the end ($\epsilon^{engn} = 0.15$) of the loading of a FCC pillar under tensile loading with three initial crystallographic orientations $\theta_0 = 0^\circ$ (top), $\theta_0 = 54.7^\circ/2$ (middle) and $\theta_0 = 65^\circ$ (bottom). The left to right panels respectively give the accumulate slip on system 1, 2 and 3.</figcaption></figure>
<figure><img src="image_12.png" /><figcaption>Figure 11: Lattice orientation evolution of a FCC porous under tensile loading with three initial crystallographic orientations $\theta_0 = 0^\circ$ (top), $\theta_0 = 54.7^\circ/2$ (middle) and $\theta_0 = 65^\circ$ (bottom) at four levels of engineering strains $\epsilon^{engn} = 0, 0.05, 0.1$ and $\epsilon^{engn} = 0.15$.</figcaption></figure>
rotations. Moreover, the obtained results provide a better understanding of the micro mechanisms leading to the evolution of the porosity under large strain, and the importance of the complex crystal plasticity related phenomenon on the evolution of the shape of the voids at the crystalline scale, which are of significant importance for the onset of ductile fracture.

## References

[1] R. J. Asaro. Crystal plasticity. *Journal of Applied Mechanics*, 50(4b):921-934, 1983.

[2] R. J. Asaro and A. Needleman. Overview no. 42 texture development and strain hardening in rate dependent polycrystals. *Acta metallurgica*, 33(6):923-953, 1985.

[3] R. J. Asaro and J. R. Rice. Strain localization in ductile single crystals. *Journal of the Mechanics and Physics of Solids*, 25(5):309-338, 1977.

[4] U. Borg and J. W. Kysar. Strain gradient crystal plasticity analysis of a single crystal containing a cylindrical void. *International Journal of Solids and Structures*, 44(20):6382-6397, 2007.

[5] O. Cazacu and I. R. Ionescu. Augmented lagrangian method for eulerian modeling of viscoplastic crystals. *Computer methods in applied mechanics and engineering*, 199(9-12):689-699, 2010.

[6] O. Cazacu and I. R. Ionescu. Dynamic crystal plasticity: An eulerian approach. *Journal of the Mechanics and Physics of Solids*, 58(6):844-859, 2010.

[7] J. Crépin, T. Bretheau, and D. Caldemaison. Cavity growth and rupture of β-treated zirconium: A crystallographic model. *Acta materialia*, 44(12):4927-4935, 1996.

[8] F. Duarte, R. Gormaz, and S. Natesan. Arbitrary lagrangian-eulerian method for navier-stokes equations with moving boundaries. *Computer Methods in Applied Mechanics and Engineering*, 193(45-47):4819-4836, 2004.

[9] S. Forest. Modeling slip, kink and shear banding in classical and generalized single crystal plasticity. *Acta materialia*, 46(9):3265-3281, 1998.

[10] M. Fortin and R. Glowinski. *Méthodes de lagrangien augmenté: applications à la résolution numérique de problèmes aux limites*. Dunod, 1982.
[11] Y. X. Gan and J. W. Kysar. Cylindrical void in a rigid-ideally plastic single crystal iii: Hexagonal close-packed crystal. *International journal of plasticity*, 23(4):592–619, 2007.

[12] Y. X. Gan, J. W. Kysar, and T. L. Morse. Cylindrical void in a rigid-ideally plastic single crystal ii: Experiments and simulations. *International Journal of Plasticity*, 22(1):39–72, 2006.

[13] R. Glowinski and P. Le Tallec. *Augmented Lagrangian and operator-splitting methods in nonlinear mechanics*. Society for Industrial and Applied Mathematics, 1989.

[14] R. Glowinski and O. Pironneau. Finite element methods for navier-stokes equations. *Annual review of fluid mechanics*, 24(1):167–204, 1992.

[15] H. J. Guo, C. Ling, E. P. Busso, Z. Zhong, and D. F. Li. Crystal plasticity based investigation of micro-void evolution under multi-axial loading conditions. *International Journal of Plasticity*, 129:102673, 2020.

[16] A. L. Gurson. Continuum Theory of Ductile Rupture by Void Nucleation and Growth: Part I—Yield Criteria and Flow Rules for Porous Ductile Media. *Journal of Engineering Materials and Technology*, 99(1):2–15, 1977.

[17] X. Han, J. Besson, S. Forest, B. Tanguy, and S. Bugat. A yield function for single crystals containing voids. *International Journal of Solids and Structures*, 50(14-15):2115–2131, 2013.

[18] F. Hecht. New development in freefem++. *Journal of Numerical Mathematics*, 20(3-4):251–265, 2012.

[19] R. Hill. Generalized constitutive relations for incremental deformation of metal crystals by multislip. *Journal of the Mechanics and Physics of Solids*, 14(2):95–102, 1966.

[20] R. Hill and J. R. Rice. Constitutive analysis of elastic-plastic crystals at arbitrary strain. *Journal of the Mechanics and Physics of Solids*, 20(6):401–413, 1972.

[21] T. J. Hughes, W. K. Liu, and T. K. Zimmermann. Lagrangian-eulerian finite element formulation for incompressible viscous flows. *Computer methods in applied mechanics and engineering*, 29(3):329–349, 1981.
[22] J. W. Hutchinson. Bounds and self-consistent estimates for creep of polycrystalline materials. *Proceedings of the Royal Society of London. A. Mathematical and Physical Sciences*, 348(1652):101-127, 1976.

[23] S. Kok, A. J. Beaudoin, and D. A. Tortorelli. A polycrystal plasticity model based on the mechanical threshold. *International Journal of Plasticity*, 18(5-6):715-741, 2002.

[24] J. W. Kysar, Y. X. Gan, and G. Mendez-Arzuza. Cylindrical void in a rigid-ideally plastic single crystal. part i: Anisotropic slip line theory solution for face-centered cubic crystals. *International Journal of Plasticity*, 21(8):1481-1520, 2005.

[25] R. A. Lebensohn and C. N. Tomé. A self-consistent anisotropic approach for the simulation of plastic deformation and texture development of polycrystals: application to zirconium alloys. *Acta metallurgica et materialia*, 41(9):2611-2624, 1993.

[26] C. Ling, J. Besson, S. Forest, B. Tanguy, F. Latourte, and E. Bosso. An elastoviscoplastic model for porous single crystals at finite strains and its assessment based on unit cell simulations. *International Journal of Plasticity*, 84:58-87, 2016.

[27] J. Mandel. Plasticité classique et viscoplasticité. *CISM Lectures Notes*, 97, 1971.

[28] V. Maronnier, M. Picasso, and J. Rappaz. Numerical simulation of three-dimensional free surface flows. *International journal for numerical methods in fluids*, 42(7):697-716, 2003.

[29] B. Maury. Characteristics ale method for the unsteady 3d navier-stokes equations with a free surface. *International Journal of Computational Fluid Dynamics*, 6(3):175-188, 1996.

[30] A. Mbiakop, A. Constantinescu, and K. Danas. An analytical model for porous single crystals with ellipsoidal voids. *Journal of the Mechanics and Physics of Solids*, 84:436-467, 2015.

[31] P. E. McHugh. Introduction to crystal plasticity theory. In *Mechanics of Microstructured Materials*, pages 125-171. Springer, 2004.

[32] M. J. Nemcko, J. Li, and D. S. Wilkinson. Effects of void band orientation and crystallographic anisotropy on void growth and coalescence. *Journal of the Mechanics and Physics of Solids*, 95:270-283, 2016.
[33] J. Paux, L. Morin, and R. Brenner. A model of porous plastic single crystals based on fractal slip lines distribution. *Journal of the Mechanics and Physics of Solids*, 167:104948, 2022.

[34] J. Paux, L. Morin, R. Brenner, and D. Kondo. An approximate yield criterion for porous single crystals. *European Journal of Mechanics-A/Solids*, 51:1–10, 2015.

[35] D. Peirce, R. J. Asaro, and A. Needleman. Material rate dependence and localized deformation in crystalline solids. *Acta metallurgica*, 31(12):1951–1976, 1983.

[36] K. H. Prakash and K. M. Shyam. Effective yield criterion for a porous single crystal accounting for void size effects. *Mechanics of Materials*, 160:103950, 2021.

[37] H. Resk, L. Delannay, M. Bernacki, T. Coupez, and R. Logé. Adaptive mesh refinement and automatic remeshing in crystal plasticity finite element simulations. *Modelling and Simulation in Materials Science and Engineering*, 17(7):075012, 2009.

[38] J. R. Rice. Inelastic constitutive relations for solids: an internal-variable theory and its application to metal plasticity. *Journal of the Mechanics and Physics of Solids*, 19(6):433–455, 1971.

[39] F. Roters, P. Eisenlohr, L. Hantcherli, D. D. Tjahjanto, T. R. Bieler, and D. Raabe. Overview of constitutive laws, kinematics, homogenization and multiscale methods in crystal plasticity finite-element modeling: Theory, experiments, applications. *Acta materialia*, 58(4):1152–1211, 2010.

[40] R. Sartorti, W. Garhuom, and A. Düster. Remeshing and eigenvalue stabilization in the finite cell method for structures undergoing large elastoplastic deformations. *Archive of Applied Mechanics*, 94(9):2745–2768, 2024.

[41] K. Sedighiani, V. Shah, K. Traka, M. Diehl, F. Roters, J. Sietsma, and D. Raabe. Large-deformation crystal plasticity simulation of microstructure and microstructure evolution through adaptive remeshing. *International Journal of Plasticity*, 146:103078, 2021.

[42] J. C. Simo and T. J. Hughes. *Computational inelasticity*, volume 7. Springer Science & Business Media, 2006.

[43] J. Smiri. *Large plastic deformation of crystals : stability analysis and attractors*. Theses, Université Paris-Nord - Paris XIII, Dec. 2024.
[44] J. Smiri, O. U. Salman, and I. R. Ionescu. Orientation attractors in velocity gradient driven processes for large plastic deformations of crystals. *Preprint, submitted*, 2025.

[45] D. Song and P. P. Castañeda. A finite-strain homogenization model for viscoplastic porous single crystals: I—theory. *Journal of the Mechanics and Physics of Solids*, 107:560–579, 2017.

[46] A. Srivastava, B. Revil-Baudard, O. Cazacu, and A. Needleman. A model for creep of porous crystals with cubic symmetry. *International Journal of Solids and Structures*, 110:67–79, 2017.

[47] R. Temam. *Navier-Stokes equations: theory and numerical analysis*, volume 343. American Mathematical Society, 2024.

[48] C. Teodosiu and F. Sidoroff. A theory of finite elastoviscoplasticity of single crystals. *International Journal of Engineering Science*, 14(2):165–176, 1976.

[49] M. Zecevic, R. A. Lebensohn, and L. Capolungo. New large-strain fft-based formulation and its application to model strain localization in nano-metallic laminates and other strongly anisotropic crystalline materials. *Mechanics of Materials*, 166:104208, 2022.

# Appendix

## A Numerical scheme used for the simulations

The main goal of this section is to recall the Eulerian numerical strategy proposed in [6] for the simplified 2D rigid-visco-plastic crystal model used in this paper.

The use of a time implicit (backward) Euler scheme for time discretization gives a set of nonlinear equations for the velocities $v$ and lattice orientation $\theta$. At each iteration in time, an iterative algorithm is developed to solve these nonlinear equations. Specifically, a mixed finite-element (FE) and Galerkin discontinuous (GD) strategy is proposed. The variational formulation for the velocity field is discretized using the finite element method, while a Galerkin discontinuous method with an upwind choice of the flux is adopted for solving the hyperbolic equations that describe the evolution of the lattice orientation. It is to be noted that in the case of the rigid-viscoplastic model studied in this work, additional
difficulties arise from the non-differentiability of the plastic terms. That means that we
cannot simply make use of the finite element techniques developed for Navier-Stokes fluids
(see for instance [14, 47]). To overcome these difficulties the iterative decomposition-
coordination formulation coupled with the augmented Lagrangian method (introduced in
[10, 13]) was modified. The reason for this modification is that, in the crystal model there
is non co-axiality between the stress deviator and the rate of deformation in contrast with
the von-Mises model for which the original method was proposed (see [10, 13]). This type
of algorithm permits also to solve alternatively, at each iteration, the equations for the
velocity field and for the unit vectors that define the lattice orientation.

For vanishing viscosity, the adopted visco-plastic model contains as a limit case the
inviscid Schmid law. Even that Schmid model is very stiff, for small viscosities (as for
metals) and moderate strain rates the iterative decomposition coordination formulation
coupled with the augmented Lagrangian method works very well and no instabilities are
presents. However, for small viscosities and large velocities (more then 100m/s) with
large Reynolds numbers (Re > 1000) the above algorithm is converging very slowly. In
this last case, other techniques, coming from gas dynamics, have to be included.

If the Eulerian domain $\mathcal{D}$ has time variations then the above algorithm has been
adapted to an ALE (Arbitrary Eulerian-Lagrangian) description of the crystal evolution.
In this case is more convenient to have the same finite FE and GD meshes. This avoid
the interpolation of the lattice orientation on the deformed mesh. As a matter of fact,
the numerical algorithm proposed here deals only with a Stokes-type problem at each
time step and the implementation of the Navier-Stokes equations in an ALE formulation
is rather standard (see for instance [8, 21, 28, 29]).

For the sake of simplicity we shall present here the numerical scheme for the simplified
$N = 3$ slips systems in 2D, described in section 2.

## A.1 Time discretization

Having in mind the initial and boundary value problem, stated in Section BVP, let $\Delta t$ be
the time step and let us denote by $v^k, \sigma^k, p^k$ and $\theta^k$ the values of the unknowns $v(k\Delta t)$,
$\sigma(k\Delta t), p(k\Delta t)$, and $\theta(k\Delta t)$. Suppose that we have computed all these variables at time
$t = (k - 1)\Delta t$. Let also denote by $V^k = V(k\Delta t), S^k = S(k\Delta t), \theta_{in}^k = \theta_{in}(k\Delta t)$ the
boundary conditions at $t = k\Delta t$.

The time implicit (backward) Euler scheme for the field equations of the initial and
boundary value problem gives the following nonlinear equations for $v^k$, $\sigma^k$, $\theta_s^k$ and $\rho^{sk}$

$$ \rho^{\text{mass}} \frac{v^k - v^{k-1}}{\Delta t} + v^k \cdot \nabla v^k - \text{div} \, \sigma'^k + \nabla p^k = \rho \mathbf{f} \quad \text{in } \mathcal{D}, \qquad (A.1) $$

$$ \text{div}(v^k) = 0 \quad \text{in } \mathcal{D}, \qquad (A.2) $$

$$ \frac{\theta^k - \theta^{k-1}}{\Delta t} + v^k \cdot \nabla \theta^k = \frac{1}{2} \left( \sum_{s=1}^{N} \dot{\gamma}_s^k - \left( \frac{\partial v_1^k}{\partial x_2} - \frac{\partial v_2^k}{\partial x_1} \right) \right), \qquad (A.3) $$

$$ D(v^k) = \sum_{s=1}^{N} \dot{\gamma}_s^k M_s^k, \qquad \dot{\gamma}_s^k = \frac{1}{\eta_s} \left[ 1 - \frac{\tau_0^{sk}}{|\boldsymbol{\sigma}^k : \boldsymbol{M}_s^k|} \right]_+ (\boldsymbol{\sigma}^k : \boldsymbol{M}_s^k), \qquad (A.4) $$

while the boundary conditions read

$$ v^k = V^k \text{ on } \Gamma_v, \quad \sigma^k n = S^k \text{ on } \Gamma_s, \quad \theta_s^k = \theta_{in}^k, \quad \text{on } \partial_{in} \mathcal{D}. \qquad (A.5) $$

## A.2 The algorithm at each time step

Let us fix the iteration in time, $k$. An iterative decomposition-coordination formulation coupled with the augmented method (see [10, 13]) will be adapted here for the crystal plasticity model. This type of algorithm permits to solve alternatively, at each iteration $n$, equations (A.1-A.4) for the velocity field and (A.3) for the lattice orientation. The convergence is achieved when the difference between $v^{k,n}$, $\sigma^{k,n}$, $\theta_s^{k,n}$ and $v^{k,n-1}$, $\sigma^{k,n-1}$, $\theta_s^{k,n-1}$ is small enough.

In order to describe the algorithm let $r > 0$ be the augmented Lagrangian step and let $\mathcal{D}$ be discretized by using a family of triangulations $(\mathcal{T}_h)_h$ made of finite elements (here $h > 0$ is the discretization parameter representing the greatest diameter of a triangle in $\mathcal{T}_h$). We denote by $V_h$ the FE space for the velocity field $v^k$, by $W_h$ the FE space for the pressures field $p$ and by $Q_h$ the Galerkin discontinuous space for the stresses deviators $\boldsymbol{\sigma}'$, for lattice orientations $\theta$ and dislocation densities $\rho_s$).

We put $v^{k,0} = v^{k-1}$, $\sigma'^{k,0} =: \sigma'^{k-1}$, $\theta_s^{k,0} = \theta_s^{k-1}$, $\dot{\delta}_s^{k,0} = \dot{\gamma}_s^{k-1}$ and we suppose that $v^{k,n-1}$, $\sigma'^{k,n-1}$ and $\theta_s^{k,n-1}$, $\dot{\delta}_s^{k,n-1}$ are known.

**Step 1.)** The first step consists in solving the following linear equation of Stokes type for the velocity field $v^{k,n}$ and the pressure $p^{k,n}$:

$$ \rho^{\text{mass}} \left( \frac{v^{k,n} - v^{k-1}}{\Delta t} + v^{k,n-1} \cdot \nabla v^{k,n-1} \right) - \text{div} \left( r \mathcal{D}(v^{k,n}) \right) + \nabla p^{k,n} = \\ \text{div} \tilde{\boldsymbol{\sigma}}^{k,n-1} + \rho^{\text{mass}} \mathbf{f}, \qquad \text{div}(v^{k,n}) = 0, \qquad (A.6) $$
with the boundary conditions

$$ \mathbf{v}^{kn} = \mathbf{V}^k \text{ on } \Gamma_v, \quad (r\mathbf{D}(\mathbf{v}^{k,n}) - p^{k,n}\mathbf{I} + \tilde{\boldsymbol{\sigma}}^{k,n-1})\mathbf{n} = \mathbf{S}^k \text{ on } \Gamma_s. $$

where we have denoted by

$$ \tilde{\boldsymbol{\sigma}}^{k,n-1} = \boldsymbol{\sigma}'^{k,n-1} - r \sum_{s=1}^{N} \dot{\delta}_{s}^{k,n-1} \mathbf{M}_{s}^{k,n-1}. $$

The above problem is a standard one in fluid mechanics and there exists many techniques to solve it. In all the computations presented in this work, we have used a Lagrangian formulation with a [continuous P2, continuous P1] choice for finite element spaces ($V_h, W_h$) associated to the velocities and pressures.

**Step 2.)** The second step consists in finding the decomposition of the global rate of deformation $\mathbf{D}(\mathbf{v}^{k,n})$, into the slip rates $\dot{\gamma}_s^{k,n} \in Q_h$, according to (A.4). This can be done using the analytic formula, which provide directly the slip rates $\dot{\gamma}_s^{k,n}$ from the expression of $\mathbf{D}(\mathbf{v}^{k,n})$ and of the yield limit $\tau_0^{sk,n-1}$, see [43].

Note that the finite element spaces for the velocity fields $\mathbf{v}^{k,n}$ and that for the slip rates $\dot{\gamma}_s^{k,n}$ cannot be chosen independently. For instance if $V_h = [\text{continuous P2}]$ then we have to choose $Q_h=[\text{discontinuous P1}]$. This is the choice in all simulations presented in this work.

**Step 3.)** We introduce now the slip rate multipliers $\dot{\delta}_s^{k,n} : \mathcal{D} \to \mathbb{R}$, belonging to $Q_h$, computed according to the decomposition-coordination formulation coupled with the augmented method for each slip system:

$$ \dot{\delta}_s^{k,n} = \frac{1}{\eta_s + r} \left[ 1 - \frac{\tau_0^{sk,n-1}}{|\boldsymbol{\sigma}^{k,n-1} : \mathbf{M}_s^{k,n-1} + r\dot{\gamma}_s^{k,n}|} \right]_+ (\boldsymbol{\sigma}^{k,n-1} : \mathbf{M}_s^{k,n-1} + r\dot{\gamma}_s^{k,n}). \quad (A.7) $$

Then the deviator stress field is updated:

$$ \boldsymbol{\sigma}'^{k,n} = \boldsymbol{\sigma}'^{k,n-1} + r \left( \mathbf{D}(\mathbf{v}^{k,n}) - \sum_{s=1}^{N} \dot{\delta}_{s}^{k,n} \mathbf{M}_{s}^{k,n-1} \right). \quad (A.8) $$

**Step 4.)** In this step, we compute the lattice orientation from the linear hyperbolic
equation for $\theta^{k,n}$

$$ \frac{\theta^{k,n} - \theta^{k-1}}{\Delta t} + \mathbf{v}^{k,n} \cdot \nabla \theta^{k,n} = \frac{1}{2} \left( \sum_{s=1}^{N} \dot{\gamma}_{s}^{k,n} - \left( \frac{\partial v_{1}^{k,n}}{\partial x_{2}} - \frac{\partial v_{2}^{k,n}}{\partial x_{1}} \right) \right), \quad (A.9) $$

with the boundary conditions

$$ \theta^{k,n} = \theta_{in}^{k}, \quad \text{on } \partial_{in}\mathcal{D}. $$

To solve the linear system (A.9) we have adopted here a GD strategy with an "upwind" choice of the flux. In the numerical applications presented in this paper we have chosen the finite volume mesh to be the finite element triangulation $\mathcal{T}_h$, and the finite volume space to be $Q_h=[\text{discontinuous P1}]$.

### A.3 The algorithm in an ALE method

If the domain $\mathcal{D}$ occupied by the single crystal (or poly-crystals) varies in time, then an arbitrary Eulerian-Lagrangian (ALE) description was adopted. We want to point out here how the above algorithm have to be changed if it is used coupled with an ALE method. For that we have to have in mind that the passage from time iteration $k-1$ to $k$ involves the frame velocity $\mathbf{v}_{fr}^{k-1}$. Since in the ALE formulation $\mathbf{v}_{fr}^{k-1} \cdot \mathbf{n} = \mathbf{v}^{k-1} \cdot \mathbf{n}$ the income boundary $\partial_{in}\mathcal{D}(t)$ is always empty. There are only three equations from the above algorithm which have to be changed: equation (A.6) with

$$ \rho^{\text{mass}} \left( \frac{\mathbf{v}^{k,n} - \mathbf{v}^{k-1}}{\Delta t} + (\mathbf{v}^{k,n-1} - \mathbf{v}_{fr}^{k-1}) \cdot \nabla \mathbf{v}^{k,n-1} \right) - \text{div} (r \mathcal{D}(\mathbf{v}^{k,n})) + \nabla p^{k,n} = \\ \text{div} \tilde{\sigma}^{k,n-1} + \rho^{\text{mass}} \mathbf{f}, \quad \text{div}(\mathbf{v}^{k,n}) = 0, \quad (A.10) $$

equation (A.9) with

$$ \frac{\theta^{k,n} - \theta^{k-1}}{\Delta t} + (\mathbf{v}^{k,n} - \mathbf{v}_{fr}^{k-1}) \cdot \nabla \theta^{k,n} = \frac{1}{2} \left( \sum_{s=1}^{N} \dot{\gamma}_{s}^{kn} - \left( \frac{\partial v_{1}^{kn}}{\partial x_{2}} - \frac{\partial v_{2}^{kn}}{\partial x_{1}} \right) \right). \quad (A.11) $$

### B Remeshing procedure

The quality of the simulation results in mesh based approaches very much depends upon the characteristics of the mesh. A poor mesh quality impacts the computational efficiency, increases the computational time and may lead to unstable solutions[18]. Meshing can be either uniform or non-uniform. Adaptive meshing is a type of non-uniform mesh scheme
that is widely used in FEM based approaches see [41]. It is characterized by a mesh
density that varies across different regions. We employ a remeshing technique to address
mesh distortion issues in crystal plasticity simulations. This method involves replacing
distorted meshes with new, undistorted ones [40]. The variables from the deformed con-
figuration are transferred to the new mesh using a nearest-neighbor mapping algorithm.
The simulation is then restarted with the initial state set based on the most recent de-
formation state. During large deformations, the aspect ratio of the elements, defined
as the ratio between element size in the stretching versus compression directions, can
become excessively large. This can introduce errors and hinder strain localization. To
mitigate this, a multi-step mesh refinement strategy is applied at each remeshing stage,
adjusting the mesh density as the deformation progresses. The key idea is to maintain a
constant number of elements in the compression direction while increasing the number of
elements in the stretching direction to preserve a near-cubic element shape. This grad-
ual mesh refinement enhances simulation resolution during deformation while minimizing
information loss and avoiding a significant increase in computational cost (see [37]).

In our adaptive meshing approach for single crystal and polycrystal simulations, we
apply three combined criteria: strain gradient, orientation gradient, and accumulative
plastic strain.

*   **Velocity Gradient Criterion** : The strain gradient serves as a crucial indicator of regions experiencing significant deformation. Areas with steep strain gradients often correspond to locations of localized plasticity, such as near grain boundaries or in regions undergoing shear band formation. By refining the mesh in these high-gradient areas, we ensure that the simulation can accurately model the material's response to applied loads and capture critical features of the deformation process.

*   **Orientation Gradient Criterion:** The orientation gradient is equally important in polycrystal simulations, as it reflects the variation in crystallographic orientation across the material. Changes in orientation can significantly influence mechanical behavior, particularly in materials exhibiting anisotropic properties. By monitoring the orientation gradient, we can identify regions where the crystal orientations change rapidly, necessitating a finer mesh to resolve the interactions between grains and accurately capture the evolution of microstructural features.

*   **Accumulative Plastic Strain Criterion:** This additional criterion targets areas of high plastic strain to minimize mesh distortion during ongoing deformation. Finer meshing in these zones helps maintain element integrity, reduce numerical artifacts,
<figure><img src="image_13.png" /><figcaption>Figure 12: Adaptive meshing evolution for a disc unit cell with a single circular void at its center (monocrystal HPC) under radial loading simulations (section 3). The mesh distributions correspond to (a) $t = T/80$ ($\epsilon^{eng} = 0.0165\%$), (b) $t = T/10$ ($\epsilon^{eng} = 0.132\%$), (c) $t = T/5$ ($\epsilon^{eng} = 0.264\%$), $t = 2T/5$ ($\epsilon^{eng} = 0.528\%$) and (d) $t = 4T/5$ ($\epsilon^{eng} = 1.056\%$).</figcaption></figure>
and enhance stability.

Combining these criteria allows for a more nuanced adaptive meshing strategy. In regions where both the strain and orientation gradients are high, the mesh density is significantly increased, ensuring precise resolution of complex interactions and behaviors. Conversely, in areas where both gradients are low, the mesh can be coarsened to optimize computational efficiency without sacrificing accuracy.

These criteria adaptive meshing approach not only enhances the fidelity of our simulations but also minimizes computational costs by focusing resources where they are most needed. As the simulation progresses, the mesh dynamically adapts to the evolving deformation field and microstructural characteristics, resulting in improved simulation accuracy for crystalline materials under large deformations.

While this multi-criteria strategy generally captures complex deformation behaviors in polycrystalline materials and improves accuracy, using all three criteria simultaneously is not always necessary or efficient. For instance, when orientation gradients are weak, a dual-criterion approach using strain gradient and accumulative plastic strain may suffice. Additionally, when strain gradient and accumulative plastic strain are both significant within the same region, the strain gradient criterion alone may be used to capture the essential deformation characteristics.

To illustrate the effectiveness of our adaptive meshing strategy, mesh refinement patterns derived from various simulations discussed in the preceding sections are presented. The figures illustrate the adaptive meshing process across different simulations of single-crystal structures under various loading conditions and strain levels. Each figure clearly demonstrates how the mesh adjusts to capture deformation characteristics specific to the crystal type, initial orientation, and applied strain.

In Figure 12 which corresponds to the section 3, the mesh begins with a dense refinement around the void, forming two overlapping star-like structures: one defined by kink bands and the other by slip bands. As deformation progresses, the kink-band structure gradually diminishes, leaving only the star pattern formed by shear bands. This evolution in mesh refinement effectively captures the localized deformation and the transition in dominant band structures surrounding the void.

Figures 13 illustrates meshing in an FCC porous under tensile loading with different crystallographic orientations ($\theta_0 = 0°$, $\theta_0 = \frac{54.7°}{2}$, and $\theta_0 = 65°$). Each orientation exhibits unique deformation responses, with mesh adaptation reflecting evolving strain distribution at each strain level.
<figure><img src="image_14.png" /><figcaption>Figure 13: Adaptive meshing evolution of a FCC porous under tensile loading (see section 4), with initial crystallographic orientations $\theta_0 = 0^\circ$ (top), $\theta_0 = \frac{54.7}{2}^\circ$ (middle) and $\theta_0 = 65^\circ$ (bottom) at three levels of engineering strains : (a) $\epsilon^{eng} = 0.05$, (b) $\epsilon^{eng} = 0.1$ and (c) $\epsilon^{eng} = 0.15$.</figcaption></figure>