# Coupling the Phase Field Method for diffusive transformations with dislocation density-based crystal plasticity: Application to Ni-based superalloys

M. Cottura, B. Appolaire, A. Finel, Y. Le Bouar

► To cite this version:

M. Cottura, B. Appolaire, A. Finel, Y. Le Bouar. Coupling the Phase Field Method for diffusive transformations with dislocation density-based crystal plasticity: Application to Ni-based superalloys. Journal of the Mechanics and Physics of Solids, 2016, 94, pp.473-489. 10.1016/j.jmps.2016.05.016. hal-01407741

HAL Id: hal-01407741

https://hal.science/hal-01407741

Submitted on 2 Jul 2021

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
# Coupling the Phase Field Method for diffusive transformations with dislocation density-based crystal plasticity: Application to Ni-based superalloys

M. Cottura¹, B. Appolaire, A. Finel, Y. Le Bouar

Laboratoire d'Etude des Microstructures, CNRS/Onera, BP72, 92322 Châtillon Cedex, France

¹Present address: CEA, DEN, Serv. Rech. Met. Phys., F-91191 Gif-sur-Yvette, France

## Abstract

A phase field model is coupled to strain gradient crystal plasticity based on dislocation densities. The resulting model includes anisotropic plasticity and the size-dependence of plastic activity, required when plasticity is confined in region below few microns in size. These two features are important for handling microstructure evolutions during diffusive phase transformations that involve plastic deformation occurring in confined areas such as Ni-based superalloys undergoing rafting. The model also uses a storage-recovery law for the evolution of the dislocation density of each glide system and a hardening matrix to account for the short-range interactions between dislocations. First, it is shown that the unstable modes during the morphological destabilization of a growing misfitting circular precipitate are selected by the anisotropy of plasticity. Then, the rafting of $\gamma'$ precipitates in a Ni-based superalloy is investigated during [100] creep loadings. Our model includes most of the important physical phenomena accounted for during the microstructure evolution, such as the presence of different crystallographic $\gamma'$ variants, their misfit with the $\gamma$ matrix, the elastic inhomogeneity and anisotropy, the hardening, anisotropy and viscosity of plasticity. In agreement with experiments, the model predicts that rafting proceeds perpendicularly to the tensile loading axis and it is shown that plasticity slows down significantly the evolution of the rafts.

**Keywords:** phase transformation, crystal plasticity, phase field modeling, superalloys

## 1. Introduction

Materials properties, especially mechanical, are strongly dependent on the microstructures, most often involving several phases. For these properties to remain permanent, the stability of the microstructures must be ensured. This requirement is particularly challenging when the materials are exposed to severe conditions, such as high temperatures, loading, or both in the worst cases. Indeed, under such conditions, microstructures are likely to evolve driven by the diffusion of alloying species and by plasticity, most often in a coupled manner.

This is the case of Ni-based superalloys used in components of engines or gas turbines submitted to high temperature. The properties of these alloys are inherited from their particular microstructure which consists of a high volume fraction of strengthening $\gamma'$ precipitates (L1₂ ordered structure) embedded in a face-centered cubic (fcc) solid-solution $\gamma$ matrix. They are optimized using thermal treatments promoting the formation of

<footnote>*Corresponding author: yann.lebouar@onera.fr</footnote>
# 1 INTRODUCTION

an initial microstructure consisting of a high volume fraction of cuboidal $\gamma'$ precipitates aligned and defining $\gamma$ channels [1] (Figure 1 - a). During service, the directional coarsening (rafting) of $\gamma'$ precipitates occurs under creep loading together with plastic deformation mainly localized within the $\gamma$ channels. This microstructure evolution is followed by a topological inversion, detrimental to the creep properties. Therefore, beside the technological issue of improving the stability of their microstructures, Ni-based superalloys appear to be a most relevant candidate for investigating strong couplings between diffusion-controlled phase transformations and plasticity.

Despite some early attempts [2, 3, 4, 5], this coupling has not been investigated extensively from a theoretical point of view, because it requires approaches with the ability to predict morphological evolutions accurately. For more than 15 years, phase field models have proved to be the most suited method for that purpose. Indeed, a number of puzzling microstructural features have been explained using this approach. Focusing on Ni-based alloys, it has brought valuable insights into the roles of elastic anisotropy [6, 7, 8, 9, 10], and chemical ordering [11, 12] in microstructure evolutions. In parallel, several works aimed at describing plasticity at the scale of individual dislocations within the phase field approach [13, 14, 15], exploiting the analogy between a dislocation loop and a platelet precipitate with a specific eigenstrain. Although these models incorporate automatically the elastic interactions between dislocations and other microstructural features, they require subnanometer grids or sophisticated technics [16] to describe properly the short-range interactions involving the properties of the dislocation cores. Moreover, it must be stressed that they lack high temperature mechanisms, such as cross-slip and climb, despite recent advances [17, 18]. To reach larger length scales, plasticity in evolving microstructures has been addressed by relying on phenomenological models. Works along this route have been only very recently proposed by several groups using mesoscale plasticity models differing by their descriptions of hardening, viscosity and plastic anisotropy. The first attempts to couple a diffuse interface model with an isotropic plasticity model have been proposed in the references [19, 20]. Similar approaches including isotropic plasticity models have been developed in [21, 22, 23, 24, 25]. In the context of rafting in Ni-based superalloys, a few works have extended the PFM to include strain gradient plasticity [26] and anisotropic plasticity, either approximately [27] or relying on a more complete and well tested crystal plasticity framework [28]. Despite some valuable achievements, for calibrating their parameters the phenomenological laws used in these models require experiments that are unfortunately more complex than usual when considering evolving multi-phase materials. Hence, it is highly desirable to resort to plasticity models with firmer physical grounds, i.e. relying on dislocation densities.

This is the aim of the present work to settle the coupling between a phase field model for diffusion-controlled transformations with a dislocation density based crystal plasticity model, in view of improving the prediction of microstructure evolutions during diffusive phase transformations, especially in Ni-based superalloys. A particular emphasis is put on two ingredients that are generally neglected but which are essential for a reliable prediction of the mechanical behavior, in particular hardening, and thus of the coupling. First, it is important to account properly for the forest hardening due to the short-range interactions between dislocations, as shown in [29]. Second, at the micrometer scale, the hardening associated with the confinement of the plastic domains can be significant, and likely to change during microstructural evolutions [26]. Thus, it is also mandatory to account for this effect in the same consistent framework, such as in [30].

The paper is structured as follows. First, the phase field model, the dislocation density based crystal plasticity model, as well as their coupling will be described in §2. At the end of this section, data and input parameters specific to Ni-based superalloys will be carefully defined and specified with respect to available experiments. Then, in §3, we will discuss calculations illustrating the hardening as described by the plasticity model associated with (i) short-range interactions of dislocations, and (ii) the back-stress introducing a size effect. In Section §4, we will investigate the concomitant evolutions of misfitting precipitates and dislocation densities in microstructures of increasing complexities: first, in the case of a single growing precipitate to address the role of plastic anisotropy on the shape stability; second, during the creep loading of a typical $\gamma/\gamma'$ microstructure in Ni-based superalloys, exhibiting rafting.
<figure><img src="image_2.png" /><figcaption>Figure 1: AM1 alloy observed by Scanning Electron Microscopy: γ' precipitates (dark) in γ matrix (light) after a standard heat treatment: a - before and b - after creep loading under 150 MPa at 950°C during 100 h.</figcaption></figure>

## 2. Model

### 2.1. Phase Field Model

At equilibrium, Ni-based superalloys feature both the disordered γ and ordered γ' phases. Following [28, 26, 10], the superalloy is modeled as an effective binary alloy. In that case, in addition to the local concentration field c(r,t), three non-conservative structural fields ηᵢ=1,3(r,t) are introduced to account for the degeneracy of the low temperature γ' phase. The four translational variants of γ' are described by the following long-range order parameters: {η₁, η₂, η₃} = η₀{1, 1, 1}, η₀{1, 1, 1}, η₀{1, 1, 1}, η₀{1, 1, 1}.

The main ingredient of phase field models is a mesoscopic free energy functional F decomposed into a chemical Fₚh (§2.1.2) and an elastic Fₚl (§2.1.3) contributions.

$$F = F_{ch}(c, \{\eta_i\}) + F_{el}(c, \{\eta_i\}, \epsilon^{el}) \quad (1)$$

where ε^{el} is the elastic strain tensor.

#### 2.1.1. Kinetic equations

The evolutions of the concentration and order parameters are governed by kinetic equations relating time derivatives to the corresponding driving forces, defined as functional derivatives of F with respect to the fields. Assuming linear constitutive relationships, a Cahn-Hilliard equation is used for the conserved concentration field and Allen-Cahn equations for the non-conserved order parameters:

$$\frac{\partial c}{\partial t}(\underline{r}, t) = M \nabla^2 \frac{\delta F}{\delta c(\underline{r}, t)} \quad (2)$$

$$\frac{\partial \eta_i}{\partial t}(\underline{r}, t) = -L \frac{\delta F}{\delta \eta_i(\underline{r}, t)} \quad (3)$$

In the present work, the kinetic coefficients M and L, related to diffusion and structural relaxation respectively, are assumed constant.
### 2.1.2. Ginzburg-Landau free energy

The chemical free energy accounts for the volume free energy associated with phase transformation and interface energies. This free energy of any volume $V$ is given by a standard Ginzburg-Landau functional:

$$F_{ch}(c, \{\eta_i\}) = \int_V f_{hom}(c, \{\eta_i\}) + \frac{\lambda}{2} |\nabla c|^2 + \frac{\beta}{2} \sum_i |\nabla \eta_i|^2 \, dV \quad (4)$$

$f_{hom}(c, \{\eta_i\})$ is the free energy density of an homogeneous system characterized by the concentration $c$ and order parameters $\eta_i$.

The gradient terms are chosen isotropic to comply with isotropic interface energies as suggested by the spherical shapes of small precipitates observed in AM1 superalloys. The coefficients $\lambda$ and $\beta$ of the gradient terms, related to the energies of $\gamma/\gamma'$ interfaces and $\gamma'/\gamma'$ anti-phase boundaries, have been prescribed constant values. As usual, the free energy density of an homogeneous system $f_{hom}(c, \eta_i)$ is approximated by a Landau polynomial expansion with respect to the order parameters. Its form is dictated by the symmetry loss associated to the $\gamma \to \gamma'$ phase transformation. Following [10], we have chosen the lowest possible order of the expansion, as well as a quadratic dependence on $c(\underline{r})$:

$$f_{hom}(c, \{\eta_i\}) = \Delta f \left[ \frac{1}{2}(c - c_{\gamma}^{0})^2 + \frac{\mathcal{B}}{6}(c_2 - c) \sum_{i=1,3} \eta_i^2 - \frac{\mathcal{C}}{3}\eta_1\eta_2\eta_3 + \frac{\mathcal{D}}{12}\sum_{i=1,3} \eta_i^4 \right] \quad (5)$$

where $\Delta f$ is an energy density scale and $c_2$ an arbitrary concentration chosen between the equilibrium concentrations $c_{\gamma}^0$ and $c_{\gamma'}^0$ of the coexisting phases. $\mathcal{B}$, $\mathcal{C}$ and $\mathcal{D}$ are constants related to $c_2$, $c_{\gamma}^0$, $c_{\gamma'}^0$, and to the equilibrium long-range order parameter $\eta_0$ (§2.3.1).

### 2.1.3. Elastic energy

In the framework of linear elasticity, the potential elastic energy reads:

$$F_{el}(\varepsilon^{el}) = F_{el}^{a}(\bar{\varepsilon}) + \frac{1}{2} \int_V \varepsilon^{el} : \lambda : \varepsilon^{el} \, dV \quad (6)$$

where $\lambda$ is the 4th order tensor of local elastic constants and where $\varepsilon^{el}$ is the second order tensor of elastic strain. $F_{el}^a(\bar{\varepsilon})$ is an homogeneous term, function of the average value $\bar{\varepsilon}$ of the total strain $\varepsilon(\underline{r})$, and depending on the loading conditions. In the present work where an applied stress $\sigma^a$ has been considered, $F_{el}^a = -V\sigma^a : \bar{\varepsilon}$. Assuming that the local concentration is the relevant field for discriminating the elastic properties, $\lambda$ is assumed to depend linearly on $c(\underline{r})$ and is thus space dependent [10]. In the small deformation framework, $\varepsilon(\underline{r})$ can be split into three contributions as follows:

$$\varepsilon(\underline{r}) = \varepsilon^{el}(\underline{r}) + \varepsilon^{0}(\underline{r}) + \varepsilon^{p}(\underline{r}). \quad (7)$$

where $\varepsilon^0(\underline{r})$ is the stress-free strain associated with the change in lattice parameter accompanying the $\gamma \to \gamma'$ transformation, and $\varepsilon^p(\underline{r})$ the plastic strain.

Assuming Vegard's law, $\varepsilon^0(\underline{r}) = \varepsilon^T \Delta c(\underline{r}) \mathbf{1}$ where $\mathbf{1}$ is the identity matrix and $\varepsilon^T = \delta/(c_{\gamma'}^0 - c_{\gamma}^0)$. The misfit $\delta$ is related to the lattice parameters $a_{\gamma}$ and $a_{\gamma'}$ of the stress-free $\gamma$ and $\gamma'$ phases $\delta' = 2(a_{\gamma'} - a_{\gamma})/(a_{\gamma'} + a_{\gamma})$.

In diffusive phase transformations, static mechanical equilibrium can safely be assumed because the relaxation of the elastic waves is by orders of magnitude faster than the evolution of $c$ and $\eta_i$. In the case of homogeneous elasticity, this problem can be solved analytically in Fourier space [31]. Otherwise, a fixed-point algorithm is used at each time step to numerically solve mechanical equilibrium [10].
## 2.2. Mesoscopic Dislocation Density-based Crystal Plasticity Model (MPM)

### 2.2.1. Plastic strain

In crystalline materials, plastic deformation usually results from the motion of dislocations. In many crystals, such as in fcc crystals, dislocations mainly glide in well-defined planes and are therefore conveniently divided into several slip systems characterized by the normal of the glide plane $\underline{n}^{\alpha}$ and the slip direction $\underline{m}^{\alpha}$. In fcc crystals, the 12 favorable (octahedral) slip systems are defined by the 4 {111} slip planes, each one containing 3 (110) slip directions.

In the small deformation framework, the plastic strain $\varepsilon^p$ is the sum of the crystallographic slips $\gamma^{\alpha}$ in each slip system $\alpha$, such that:

$$ \varepsilon^p = \sum_{\alpha} \gamma^{\alpha} \underline{P}_{\alpha} \qquad (8) $$

where $\underline{P}_{\alpha}$ is the symmetric Schmid tensor defined by $\underline{P}_{\alpha} = \frac{1}{2}(\underline{m}^{\alpha} \otimes \underline{n}^{\alpha} + \underline{n}^{\alpha} \otimes \underline{m}^{\alpha})$.

### 2.2.2. Resolved shear stresses

Glide of dislocations belonging to the system $\alpha$ is driven by the resolved shear stress $\tau^{\alpha} = \underline{\sigma} : \underline{P}^{\alpha}$ where $\underline{\sigma}$ is the stress tensor. More precisely, dislocations will glide only if the resolved stress exceeds a certain threshold, the critical stress $\tau_c^{\alpha}$, which depends on the Peierls stress and on short-range interactions with dislocations belonging to the same or other glide systems. In fcc crystals, the dissociation of the dislocation core leads to a very low Peierls stress and this quantity is usually neglected.

The critical stress $\tau_c^{\alpha}$ is expressed using a generalized Taylor formula which accounts for the anisotropy of the interactions between slip systems in fcc crystals [32]:

$$ \tau_c^{\alpha} = \mu b \sqrt{\sum_{\beta} a^{\alpha\beta} \rho^{\beta}} \qquad (9) $$

with b the norm of the Burgers vector, $\mu$ an average shear modulus and $\rho^{\beta}$ the total dislocation density on the slip system $\beta$. $a^{\alpha\beta}$ is the hardening matrix: $\sqrt{a^{\alpha\beta}}$ is the average strength of the interaction between the two slip systems $\alpha$ and $\beta$. In a fcc crystal, the matrix has $12^2$ coefficients that can be reduced to 6 independent coefficients using symmetries. These coefficients are associated with six types of interaction [33]: dislocations gliding on parallel slip systems, involve the self-interaction ($a_0$) and coplanar ($a_1^{copl}$) interaction coefficients. Three other coefficients account for the interactions between non-coplanar slip systems, forming Hirth locks ($a_1^{orth}$), glissile junctions ($a_2^{glis}$) and Lomer-Cottrell locks ($a_3^{LC}$). The last coefficient, $a_1^{col}$, is related to collinear interactions that generate a strong hardening effect due to the partial annihilation of gliding collinear dislocation segments and the resulting generation of highly curved dislocations that require a higher stress to recover a flowing state. We have employed the values for fcc crystals inferred in [34] from dislocation dynamics simulations and given in §2.3.2.

### 2.2.3. Stockage and recovery law

The evolution of the dislocation density is given by a standard storage-recovery balance equation [35, 36]. Following Teodosiu et al. [37], the time derivative of the dislocation density is expressed as a function of the plastic slip rates $\dot{\gamma}^{\alpha}$ using Orowan's relation as follows:

$$ \dot{\rho}^{\alpha} = \frac{1}{b} \left( \frac{1}{L^{\alpha}} - 2 y_c \rho^{\alpha} \right) |\dot{\gamma}^{\alpha}| \qquad (10) $$

where $L^{\alpha}$ is the dislocation mean free path on slip system $\alpha$ and $y_c$ is a characteristic length associated with the annihilation process. $L^{\alpha}$ depends on the density of obstacles encountered by dislocations, i.e. on
the density of dislocations belonging to other slip systems. Following Kubin et al. [38], $L^\alpha$ is assumed to be inversely proportional to the critical stress on the system:

$$L^\alpha = K / \sqrt{\sum_\beta a^{\alpha\beta} \rho^\beta} \qquad (11)$$

where K is a material parameter.

### 2.2.4. Plastic flow rule

Following a standard approach, the single crystal kinematics is related to the underlying dislocation density evolution with a phenomenological viscoplastic power law relating the plastic slip rates $\dot{\gamma}^\alpha$ to the resolved shear stress $\tau^\alpha$, the backstress $\tau_b^\alpha$ and the slip hardening $\tau_c^\alpha$ given by [39]:

$$\dot{\gamma}^{\alpha} = \dot{\gamma}_0 \left( \frac{|\tau^{\alpha} - \tau_b^{\alpha}|}{\tau_c^{\alpha}} \right)^{\frac{1}{m}} \exp \left( -\frac{\Delta G^{\alpha}}{k_B T} \right) \mathrm{sign} (\tau^{\alpha} - \tau_b^{\alpha}) \qquad (12)$$

where $\dot{\gamma}_0$ and $m$ are the reference plastic shear rate and rate sensitivity exponent, respectively. T is the absolute temperature and $k_B$ Boltzmann's constant. The resistance to slip in system $\alpha$ is reproduced by the critical resolved shear stress $\tau_c^\alpha$ discussed above, as well as by some back stress $\tau_b^\alpha$ to account for a size dependency of plastic flow as discussed below. The exponential term corresponds to the thermally-activated dislocation motions where the activation energy $\Delta G^\alpha$ is defined by:

$$\Delta G^{\alpha} = G^{0} \left( 1 - \frac{|\tau^{\alpha} - \tau_{b}^{\alpha}|}{\tau_{c}^{\alpha}} \right) \qquad (13)$$

$G_0$ is the free energy, usually assumed constant, that is required for a dislocation to overcome a short-range barrier without any external work.

### 2.2.5. Back-stress

Back-stress terms are usually involved in continuous formulation of plasticity. These terms, as well as the local critical stress, are reminiscent of the short-range dislocation-dislocation correlations that are not resolved at the scale at which the continuous model is elaborated [40, 41, 42]. In other words, the back and critical stresses $\tau_b^\alpha$ and $\tau_c^\alpha$ should emerge from the coarse graining procedure that is needed to make a link between the microscopic and mesoscopic scales. Such a procedure has already been elaborated in the simplest situation where a single two-dimensional glide system is considered [42]. However, this coarse graining procedure has still not been performed in the general situation where multiple three dimensional glide systems are simultaneously active. Therefore, following Ref. [30], we use for the present work the following back-stress expression:

$$\tau_b^\alpha = \frac{\mu b R^2}{8} \left[ \frac{1}{(1-\nu)} \nabla \rho_{G,e}^\alpha \cdot \underline{m}^\alpha - 2 \sum_\beta \delta_s^{\alpha\beta} \left( \nabla \rho_{G,s}^\beta \cdot \underline{p}^\alpha \right) \right] \qquad (14)$$

with

$$\delta_s^{\alpha i} = \begin{cases} 1 & \text{for } (\alpha \ i) = (4,13), (6,18), (8,17), (9,15), (10,16), (11,14) \\ (1,16), (2,17), (3,18), (5,14), (7,13), (12,15) & \\ 0 & \text{otherwise} \end{cases} \qquad (15)$$

where $\mu$ and $\nu$ are the shear modulus and Poisson ratio, respectively, for isotropic materials, and where $\underline{p}^\alpha = \underline{m}^\alpha \times \underline{n}^\alpha$ (edge dislocations are then parallel to $\underline{p}$ and screw to $\underline{m} = \underline{b}/b$. $\rho_{G,s/e}^\alpha$ corresponds to the GND density with a screw/edge character on slip system $\alpha$.
<table><caption>Table 1: Description of glide systems in fcc crystals.</caption><thead><tr><th>Glide system α</th><th>Dislocation i</th><th>Type</th><th>√2 <u>m</u></th><th>√3 <u>n</u></th><th>Schmidt & Boas</th></tr></thead><tbody><tr><td>1</td><td>1</td><td>edge</td><td>[011]</td><td>[111]</td><td>A2</td></tr><tr><td>2</td><td>2</td><td>edge</td><td>[101]</td><td>[111]</td><td>A3</td></tr><tr><td>3</td><td>3</td><td>edge</td><td>[110]</td><td>[111]</td><td>A6</td></tr><tr><td>4</td><td>4</td><td>edge</td><td>[011]</td><td>[111]</td><td>B2</td></tr><tr><td>5</td><td>5</td><td>edge</td><td>[101]</td><td>[111]</td><td>B4</td></tr><tr><td>6</td><td>6</td><td>edge</td><td>[110]</td><td>[111]</td><td>B5</td></tr><tr><td>7</td><td>7</td><td>edge</td><td>[011]</td><td>[111]</td><td>C1</td></tr><tr><td>8</td><td>8</td><td>edge</td><td>[101]</td><td>[111]</td><td>C3</td></tr><tr><td>9</td><td>9</td><td>edge</td><td>[110]</td><td>[111]</td><td>C5</td></tr><tr><td>10</td><td>10</td><td>edge</td><td>[011]</td><td>[111]</td><td>D1</td></tr><tr><td>11</td><td>11</td><td>edge</td><td>[101]</td><td>[111]</td><td>D4</td></tr><tr><td>12</td><td>12</td><td>edge</td><td>[110]</td><td>[111]</td><td>D6</td></tr><tr><td>4 or 7</td><td>13</td><td>screw</td><td>[110]</td><td>[111] or [111]</td><td>A6 or D6</td></tr><tr><td>5 or 11</td><td>14</td><td>screw</td><td>[101]</td><td>[111] or [111]</td><td>A3 or C3</td></tr><tr><td>9 or 12</td><td>15</td><td>screw</td><td>[011]</td><td>[111] or [111]</td><td>C1 or D1</td></tr><tr><td>1 or 10</td><td>16</td><td>screw</td><td>[110]</td><td>[111] or [111]</td><td>B5 or C5</td></tr><tr><td>2 or 8</td><td>17</td><td>screw</td><td>[010]</td><td>[111] or [111]</td><td>B4 or D4</td></tr><tr><td>3 or 6</td><td>18</td><td>screw</td><td>[011]</td><td>[111] or [111]</td><td>A2 or B2</td></tr></tbody></table>

In the spirit of Ref. [30], the length scale *R* should be of the same order as the discretisation scale *d* used in the implementation of the model. Consequently we have chosen *R* = 1.5*d*

The GND density is evaluated from the plastic slip gradients assuming that all dislocations are screw or edge (no mixed dislocations). The 18 types of GND are listed in Table 1. Concerning screw or edge dislocations, the GND densities are related to the plastic slip gradients throughout:

$$ \begin{align} \rho_{G,e}^{\alpha} b &= - \nabla \gamma^{\alpha} \cdot \underline{m}^{\alpha} \nonumber \\ \rho_{G,s}^{\alpha} b &= \nabla \gamma^{\alpha} \cdot \underline{p}^{\alpha} \tag{16} \end{align} $$

Finally, using equations (14) and (16), the calculation of the back-stress involves the gradients of the crystallographic plastic slip:

$$ \tau_b^\alpha = -\frac{\mu R^2}{8} \left[ \frac{1}{1-\nu} (\underline{m}^\alpha \otimes \underline{m}^\alpha) : (\nabla \otimes \nabla \gamma^\alpha) + 2 \sum_\beta \left( \sum_i \delta_s^{\alpha i} \delta_s^{\beta i} \underline{p}^\beta \otimes \underline{p}^\alpha \right) : (\nabla \otimes \nabla \gamma^\beta) \right] \quad (17) $$

### 2.2.6. Coupling PFM & MPM

The MPM is coupled to the PFM presented in §2.1 through the total strain field:

$$ \varepsilon(\underline{r}) = \varepsilon^{el}(\underline{r}) + \varepsilon^{0}(\underline{r}) + \varepsilon^{p}(\underline{r}) \quad (18) $$
<table><caption>Table 2: Physical parameters for AM1 superalloys at T = 950°C.</caption><thead><tr><th>D (m<sup>2</sup>.s<sup>-1</sup>)</th><th>σ<sub>exp</sub> (mJ.m<sup>-2</sup>)</th><th>δ</th><th>C<sub>11</sub> (GPa)</th><th>C<sub>12</sub> (GPa)</th><th>C<sub>44</sub> (GPa)</th></tr></thead><tbody><tr><td>3 10<sup>-19</sup></td><td>4 [43, 9]</td><td>-0.1% [44, 45]</td><td>γ</td><td>197</td><td>144</td><td>90</td></tr><tr><td></td><td></td><td></td><td>γ'</td><td>193</td><td>118</td><td>124</td></tr></tbody></table>

where $\varepsilon^p(\underline{r})$ is the plastic strain field evolving according to Eqs. (8) and (12).

In heterogeneous materials, the viscoplastic parameters $a_{\alpha\beta}$, $G_0$, $m$, $\dot{\gamma}_0$ and $R$ are space dependent. In the ranges of temperature and stress considered, only $\gamma$ undergoes plastic strain while $\gamma'$ behaves elastically. To reproduce this behavior, all viscoplastic parameters have been set at their values in the $\gamma$ phase and the governing equations $\rho^\alpha$ and $\gamma^\alpha$ are only activated inside $\gamma$. This is achieved by ponderating the rates $\dot{\rho}^\alpha$ and $\dot{\gamma}^\alpha$ in Eqs. (9) and (12) by a sigmoid function $j(\phi)$ of the order parameters:

$$j(\phi) = \frac{1}{2} \left( 1 - \tanh \left[ \theta \left( \phi - \frac{1}{2} \right) \right] \right) \qquad (19)$$

with $\phi = \frac{1}{3\eta_0^2} \sum_{i=1,3} \eta_i^2$ and $\theta$ a parameter controlling plasticity in the interface.

## 2.3. Input parameters

### 2.3.1. Physical properties and phase field parameters

The physical properties and elastic constants used in the subsequent calculations are reported in Table 2. The inhomogeneous and anisotropic elastic constants have been retrieved from [46], where it has been shown that a large inhomogeneity of the shear modulus $C' = (C_{11} - C_{12})/2$ is required to generate well-aligned cuboidal precipitates with aspect ratios close to the experimental observations of AM1

Following [28], $M$ is such as to recover the interdiffusion coefficient in $\gamma$: $D = D_0 \exp(-\Delta U/k_B T)$ with $D_0 = 1.45 \ 10^{-4} \ m^2s^{-1}$ and $\Delta U = 2.8 \ eV$ [47]. We have used $L = 100 \ Md^{-2}$, where $d$ is the grid spacing, to ensure that kinetics is much faster for the order parameters than for the concentration field. We have used the equilibrium concentrations $c_\gamma^0 = 0.15$ and $c_{\gamma'}^0 = 0.231$ identified on the Ni-Al phase diagram at 950°C and set $c_2 = 0.18$. Following [10], an equilibrium order parameter $\eta_0$ set to 1 gives $B = 0.162$, $C = 0.01458$ and $D = 0.022842$ in Eq. (5). The non-dimensional gradient coefficients $\tilde{\lambda} = \lambda/(\Delta fd^2)$ and $\tilde{\beta} = \beta/(\Delta fd^2)$ are chosen as $\tilde{\lambda} = 0.21$ and $\tilde{\beta} = 9.75 \ 10^{-4}$. This choice ensures that: (i) the interface width is equal to 6d, i.e. large enough to avoid pinning on the numerical grid, (ii) antiphase boundaries in the $\gamma'$ phase are wetted by the $\gamma$ matrix, and (iii) the interface energy contributions due to the gradients of order parameter and concentration fields are equal. Then, using 1D stress-free simulations, we compute the numerical interface energy $\tilde{\sigma} = 0.66 \ 10^{-3}$ from the non-dimensional profile $f_{hom}/\Delta f$, $\tilde{\lambda}$ and $\tilde{\beta}$. $\tilde{\sigma}$ is related to the experimental interface energy through:

$$\sigma_{exp} = \tilde{\sigma} \, d \, \Delta f \qquad (20)$$

For a given $d$, $\Delta f$ can be deduced from Eq. (20).

### 2.3.2. Viscoplastic parameters

The coefficients of the hardening matrix $a_{\alpha\beta}$ for superalloys are given in reference [34]. The values have been determined from dislocation dynamic simulations in a fcc crystal [33]. In Eq. (14), the shear modulus $\mu$ has been taken equal to $\sqrt{C'C_{44}}$ where $C_{44}$ and $C'$ are the two shear moduli of the $\gamma$ phase.

The viscoplastic parameters are identified on the macroscopic mechanical behavior of the bulk $\gamma$ phase using the experimental data of Espié et al. [48]. For that purpose, we have considered [100] tensile tests performed
## 3 MECHANICAL BEHAVIOR OF FROZEN TWO-PHASE MICROSTRUCTURES

at 950°C and two strain rates $\dot{\epsilon}_{11} = 9 \times 10^{-4}$ and $9 \times 10^{-3} \text{ s}^{-1}$ on bulk $\gamma$ phase with a composition close to the $\gamma$ phase in AM1. The values of $y_c$, $\dot{\gamma}_0$ and $m$ reported in Table 3 result from the best fit with the calculations considering the activation energy $G_0 = 4.54 \times 10^{-20} \text{ J}$ proposed in [30] and the same initial dislocation density on all slip systems $\bar{\rho}_0 = 2.5 \times 10^{12} \text{ m}^{-2}$ consistent with experimental assessments [49, 50].

Finally, to confine viscoplasticity in the $\gamma$ phase only, we have set $\theta = 100$. With this value, the viscoplastic parameters vary over one grid spacing $d$, a distance smaller than the interface width.

<table><caption>Table 3: Viscoplastic parameters of the $\gamma$ phase in the AM1 superalloys.</caption><thead><tr><th>b (nm)</th><th>y<sub>c</sub> (nm)</th><th>K</th><th>a<sub>0</sub></th><th>a<sub>1</sub><sup>copl</sup></th><th>a<sub>1</sub><sup>orth</sup></th><th>a<sub>1</sub><sup>coll</sup></th><th>a<sub>2</sub><sup>glis</sup></th><th>a<sub>3</sub><sup>LC</sup></th><th>μ (GPa)</th><th>ν</th><th>γ<sub>0</sub> (s<sup>-1</sup>)</th><th>m</th><th>G<sub>0</sub> (J)</th></tr></thead><tbody><tr><td>0.256</td><td>1.6</td><td>80</td><td>0.12</td><td>0.12</td><td>0.06</td><td>0.62</td><td>0.11</td><td>0.12</td><td>48</td><td>0.31</td><td>0.001</td><td>0.04</td><td>4.54 10<sup>-20</sup></td></tr></tbody></table>

## 3. Mechanical behavior of frozen two-phase microstructures

Before analyzing the consequences of the full coupling between PFM and MPM on diffusion-controlled microstructure evolutions, it is important to have a clear idea of the effect of the two ingredients usually neglected in existing PFM-MPM couplings, i.e. the hardening matrix and back stress. For that purpose, we illustrate sequentially their roles in situations simple enough to highlight the main features likely to be encountered in the more complex situations of §4.

### 3.1. Hardening matrix

We consider the simple case of a homogeneous $\gamma$ crystal under tensile loading with a deformation rate $\dot{\epsilon}_{11} = 9 \times 10^{-3}$. First, all glide systems are attributed the same initial homogeneous dislocation density $\rho^\alpha = \bar{\rho}_0 = 3.8 \times 10^{10} \text{ m}^{-2}$. Under a [100] load, 8 glide systems are activated with the same Schmid factor (A3, A6, B4, B5, C3, C5, D4, D6).

<figure><img src="image_3.png" /><figcaption>Figure 2: Left: dislocation densities vs. plastic strain $\epsilon_{11}^p$ during tensile [100] loading for homogeneous (continuous lines without symbol) and random (dashed lines + continuous lines with symbols) initial distributions of dislocation densities on the different glide systems (Tab. 1). Right: corresponding stress-strain curves $\sigma_{11}$ vs. $\epsilon_{11}$.</figcaption></figure>

In this particular configuration, the dislocation densities evolve identically on the activated slip systems: their value increases from the average $\bar{\rho}_0$ up to $2.4 \times 10^{11} \text{ m}^{-2}$ during 0.7 s (Fig. 2 - left, full black line). In
a real material, it is unlikely that all slip systems feature the same initial dislocation density. Thus, the different slip systems are attributed initial dislocation densities randomly picked in a normal distribution with mean value $\bar{\rho}_0 = 3.8 \ 10^{10} \ \text{m}^{-2}$ and a relative standard deviation of 0.2. The evolutions of the dislocation densities in all glide systems are shown in Figure 2 (left) during the tensile loading.

As mentioned previously, 4 glide systems are not activated (horizontal dashed lines, C1 and D1 are superposed), but the evolution of the densities in the 8 activated glide systems are not equivalent anymore. The densities in C3, A6, D4 and B5 (bullets) increase more rapidly than in their corresponding collinear systems A3, D6, B4 and C5 (squares). Moreover, the more active a system is, the less its collinear system is. This behavior is in agreement with [51], where fcc single crystals have been investigated by both dislocation dynamics simulations and continuum modelling. It is due to the strong short range interactions between dislocations on collinear systems, as mentioned above. In fcc crystals, this interaction has been shown to be the strongest of all interactions between slip systems [29]. In our model, these interactions are accounted for with $a_1^{\text{coll}}$ in the hardening matrix appearing in Eq. (9). Finally, the impact of the initial dislocation densities on the mechanical behavior can clearly be observed on the tensile stress-strain curves displayed in Figure 2 (right): the case with random initial dislocation densities features a yield stress by a few percents higher than the one with identical densities, as well as a larger hardening that increases the stress difference between both cases during the tensile loading.

## 3.2. Role of the back-stress

<figure><img src="image_4.png" /><figcaption>Figure 3: Simple 1D configuration.</figcaption></figure>

Following [26, 52], the role of the back-stress is illustrated in a simple 1D static laminate made of layers of $\gamma$ and $\gamma'$ (Fig. 3). To comply with Ni-based superalloys during low-stress high temperature creep [53], it is assumed that only $\gamma$ phase can deform plastically $\gamma'$ responding only elastically.

<figure><img src="image_5.png" /><figcaption>Figure 4: Left: Stress-strain curve $\sigma_{11}$ vs. $\varepsilon_{11}$. Right: profiles of $\varepsilon_{11}^p$ at $\bar{\varepsilon}_{11}^p = 0.2\%$ for laminates with different $\gamma'$ channel widths $s$.</figcaption></figure>

We have investigated the mechanical behavior of a 1D laminate constituted of $\gamma$ (channels) and $\gamma'$ (precipitates) with interfaces normal to some cubic directions of the crystal lattices (Fig. 3). We have considered several laminates with different width $s$ ranging from 10 to 200 nm for the $\gamma$ channels, keeping constant the volume fraction of $\gamma' \tau_{\gamma'} = 0.7$. They all have been discretized with 128 nodes, and initialized with an homogeneous dislocation density $\bar{\rho}_0 = 2.5 \ 10^{12} \ \text{m}^{-2}$ on all slip systems. Moreover, they have been submitted
to tensile strain normal to the layers with rate $\dot{\varepsilon}_{11} = 10^{-3} \text{ s}^{-1}$. For each channel size $s$, calculations without and with back-stress have been carried out. The corresponding tensile stress-strain curves are plotted in Figure 4 (left), together with the purely elastic cases (black dashed line). As expected, whatever the size of the $\gamma$ channel, the elastic response is the same because the volume fraction is kept constant.

When neglecting back-stress (black line), the behavior is the same for a constant $\tau_{\gamma'}$: the yield stress and hardening do not depend on $s$. On the contrary, when considering back-stress, the plastic regime depends on $s$, and the two-phase materials is strengthened as the channel width decreases: the yield stress as well as the hardening increase when $s$ decreases. The 0.2% offset yield stress increases from 480 MPa for an infinite width (or no backstress), to 530 MPa for $s = 200$ nm and 620 MPa for $s = 20$ nm.

The differences in mechanical behavior originate from the distribution of plastic deformation confined in the $\gamma$ channel. Typical profiles of $\varepsilon_{11}^p$ are shown in Figure 4 (right) for different values of $s$ at the same average plastic strain $\varepsilon_{11}^p = 0.2\%$. When $\tau_b = 0$, $\varepsilon_{11}^p$ is homogeneous in the whole $\gamma$ channel for every values of $s$. When $\tau_b \neq 0$, $\varepsilon_{11}^p$ is inhomogeneous and increases from 0 at the $\gamma/\gamma'$ interface to a maximum value at the center of the channel over a plastic length $\xi_P$. When the channel width is much larger than $\xi_P$ (e.g. $s = 100$ nm), the plastic deformation remains almost constant over most of the $\gamma$ channel. When these two lengths are of the same order of magnitude (e.g. $s = 10$ nm), $\varepsilon_{11}^p$ exhibits profiles nearly parabolic. According to Eq. (17), this kind of inhomogeneous distribution increases the back-stress and contributes to the strengthening of two-phase materials. It is worth stressing that although the present model delivers outcomes similar to the isotropic strain gradient plasticity model used in [26] in the present symmetrical configuration, more complex configurations (e.g. loadings) would induce different behaviors. As a conclusion, our model correctly reproduces the strengthening of the material when the size of the plastic regions decreases. When studying evolving microstructures coupled with plasticity, this effect has to be included because the size of the plastic regions may evolve with time. This is in particular the case when studying microstructure evolutions in single crystal nickel based superalloys.

## 4. Evolving two-phase microstructures

### 4.1. Growth of a single misfitting precipitate in a plastic matrix

First, we investigate in what respect anisotropic plasticity may influence morphological evolutions. Indeed, it has been shown recently that the influence of isotropic plasticity on phase transformations is very dependent on the eigenstrain [54], morphology [25], loading [55] and whether plastic strain is inherited by the growing phase or not [56]. Hence, it is of valuable interest to examine how the anisotropic plasticity model introduced above relaxes stresses generated by the transformation itself in a simple case, before handling the case of rafting involving complex initial microstructures and an external loading.

For that purpose, we consider a misfitting $\gamma'$ precipitate growing in a supersaturated $\gamma$ matrix. To emphasize the role of the anisotropy of plastic deformation, we have changed the following parameters with respect to Tab. 2 such that (i) the elastic energy is isotropic when plasticity is discarded and such that (ii) the misfit stress promotes significant plastic flow: homogeneous isotropic elastic constants with Young's modulus $E = 93$ GPa and Poisson's ratio $\nu = 0.37$ ; and an isotropic eigenstrain $\varepsilon^T = 0.2\%$.

A $3.2 \times 3.2 \ \mu\text{m}^2$ stress-free periodic box aligned with the [100] and [010] cubic axes of the fcc lattice and corresponding to a 3D system invariant along [001] is discretized into $512^2$ nodes ($d = 6$ nm). A circular precipitate with radius $r_0 = 30$ d, and with the equilibrium concentration of the phase diagram $c_{\gamma'}^0 = 0.231$, is included in a matrix with 18.5% of supersaturation ($c_\infty = 0.165 > c_\gamma^0$). To start with smooth relaxed interfaces, we let the microstructure evolve without any plasticity during 0.1 h. Then, dislocation densities are homogeneously distributed on all glide systems with $\bar{\rho}_0 = 2.0 \ 10^{12} \ \text{m}^{-2}$.

The evolutions in three cases are compared in Fig. 5 where the concentration fields are plotted at 4 successive times. The first row displays the evolution when both phases are purely elastic. The evolution of the second
<figure><img src="image_6.png" /><figcaption>Figure 5: Concentration field c at t = 0.1, 0.75, 1.5 and 3 h obtained with the elastic (top) and elasto-viscoplastic model with: K = 80 (middle) and K = 150 (bottom).</figcaption></figure>

row ensues from the viscoplastic model with the viscoplastic parameters in Tab. 3, in particular K = 80. In the third row, the plastic parameter K is increased up to 150 (rather than 80 in Tab. 3) in order to increase the dislocation mean free path (Eq. 11) and therefore favor plastic flow. In all cases, the chemical driving force induced by the supersaturation leads to the growth of the precipitate. In the elastic case (top row), the precipitate remains circular as it grows (r = 102 d after 3 h), consistently with the isotropy of interface energy, elasticity and diffusion. In the viscoplastic case with K = 80 (middle row), the growth is slightly enhanced along the axes and diagonals of the box to form approximately an octagon with facets close to {120} planes. In the viscoplastic case with K = 150 (bottom row), the anisotropic growth is more pronounced. The shape adopts rapidly an octagonal shape (after 0.75 h rather than 3 h in the previous case), before evolving towards a nearly square shape with sides aligned along the {110} planes. These simulations clearly demonstrate that plastic anisotropy is able to induce a shape change during the growth of a precipitate. They also show that the anisotropy of plasticity influences significantly the anisotropy of the resulting shape in a non trivial manner as we will discuss now.

Plasticity of the fcc matrix results from the anisotropic activity of the twelve slip systems. As shown in Figure 6, the primary slip systems are A6, D6, B5, C5, and the associated dislocation densities increase along the cubic axes to give a cross-like pattern. The eight secondary slip systems are mainly active along two opposite directions, at angles θ with respect to [100] direction: θ is close to π/8 for B4 and C3, 3π/8 for B2 and C1, 5π/8 for A2 and D1 and 7π/8 for A3 and D4. In addition, these systems are weakly active at θ + π/2.

The very first moments of the activation of the fcc slip systems can be predicted from the analysis of the stress fields generated by a cylindrical precipitate of radius $r_i$ in a cylindrical matrix of external radius $r_e$ with finite boundaries. Assuming isotropic elastic constants, the matrix features the following non-zero stress components in the cylindrical frame $(r, \theta, z)$:

$$ 
\begin{align} 
\sigma_{rr} &= (a^2 - b^2) p_e \nonumber \\ 
\sigma_{\theta\theta} &= (a^2 + b^2) p_e \tag{21} \\ 
\sigma_{zz} &= 2 a^2 p_e \nonumber 
\end{align} 
$$
where $a = r_i/r_e$, $b = r_i/r$ and $p_e = E \varepsilon^T/(2(1 - \nu))$. It appears that the maximum reduced shear stress $2p_e/\sqrt{6}$ is obtained for the primary slip systems A6, D6, B5, C5 along the cubic directions ($\theta = 0 \mod \pi/2$). At the interface, the reduced stresses in all secondary slip systems feature two local maxima equal to either $(\sqrt{2} + a^2)p_e/\sqrt{6}$ or $(\sqrt{2} - a^2)p_e/\sqrt{6}$. A careful analysis reveals that the angle $\theta$ corresponding to the global maximum is $\pi/8$ for B4 and C3, $3\pi/8$ for B2 and C1, $5\pi/8$ for A2 and D1 and $7\pi/8$ for A3 and D4. These directions sketched in Figure 6 (bottom right) and are fully consistent with the patterns of the corresponding dislocation densities.

As shown in Figure 5, the influence of viscoplasticity on the phase fraction is small. However, a significant increase of the growth velocity in the cubic directions is observed when viscoplasticity is included. As shown in Figure 6, this rapid growth is related to the activation of the four primary slip systems A6, D6, B5 and C5. Along $\theta = 0 \mod \pi/2$, the strain field generated by a single circular precipitate in the matrix has only two non-zero components $\epsilon_{11} = -\epsilon_{22}$. Thus full plastic relaxation cannot be achieved by a single primary slip system because the associated plastic strain $\epsilon_{ij}^p = (b_i n_j + b_j n_i)/2$ contains non diagonal components. For example, the plastic strain generated by the slip systems A6 and D6 are respectively:

$$ \varepsilon^p(A6) = \frac{2}{\sqrt{6}} \begin{pmatrix} -1 & 1/2 & 1/2 \\ 1/2 & 1 & 0 \\ 1/2 & 0 & 0 \end{pmatrix} ; \quad \varepsilon^p(D6) = \frac{2}{\sqrt{6}} \begin{pmatrix} 1 & 1/2 & 1/2 \\ 1/2 & -1 & 0 \\ 1/2 & 0 & 0 \end{pmatrix} \qquad (22) $$

However, if the systems A6 and D6 are activated with the same intensity but different signs, the resulting total plastic strain is diagonal with only two non-zero components $\epsilon_{11}^p = -\epsilon_{22}^p$ to achieve full plastic relaxation. This explains the equivalent activation of the collinear systems A6 and D6 despite the strong short range interactions between them. A similar conclusion can be drawn for the systems B5 and C5.

<figure><img src="image_7.png" /><figcaption>Figure 6: Dislocation densities $\rho - \bar{\rho}_0$ in slip systems A2, A3, B2, B4 (secondary) and A6 (primary) at $t = 3$ h with $K = 150$. The dashed white lines correspond to the interface position indicated by the isovalue $\eta_{i=1,3} = 0.5$. Bottom right: orientations of the slip systems (Tab. 1) featuring maximum reduced shear stresses (primary in black, secondary in blue and red).</figcaption></figure>

The activation of the primary slip systems induces plastic strain whose major components are $\epsilon_{11}^p = -\epsilon_{22}^p$. The components of the plastic strain display the same cross pattern as the primary dislocations densities
(Fig. 7 - a). Within the cross pattern, three zones can be distinguished as shown in Fig. 7 - b where the profile of $\epsilon_{22}^p$ along the [100] axis (dashed dotted line in Fig. 7 - a) is plotted in blue. Zone 1 corresponds to the circular precipitate which has grown before plasticity has been activated (delimited by the inner dashed line in Fig. 7 - a, and the light grey slab in Fig. 7 - b). Zone 2 spans over the branches of the cross in Fig. 7 - a inside the precipitate and features an increase of plastic strain along the radius to achieve a maximum at the interface (outer dashed line in Fig. 7 - a, and darker grey slab in Fig. 7 - b). Zone 3 lies in the matrix and displays a decrease of plastic strain on a short distance (~20 d), with longer tails observed in Fig. 7 - a at about 45° from the cubic directions. The profile at a time equal to 1.5 h is added in Fig. 7 - b (dashed line) to illustrate the process responsible for the increase of plastic strain in zone 2, although the precipitate does not undergo any plastic flow. As soon as plasticity is activated, plastic flow is the most pronounced at the interface in the matrix because the resolved shear stress is the higher there. During growth, the interface sweeps the border of zone 3 and the plastic strain is included in the precipitate. Once included, the plastic strain does not evolve any longer but its progressive build-up along the precipitate radius generates increasing shear stresses in the matrix at the interface which promote more plastic flow.

Even if the precipitate does not deform plastically, inhomogeneous plastic strain fields are present within the precipitate because of the inheritance of plasticity in the present model. This also leads to inhomogeneous distribution of stresses within the precipitate. Finally, the amount of plastic strain inherited by the $\gamma'$ phase from the matrix can be seen to increase when the particle size increases.

This is in line with the calculations of Ammar et al. [56] when isotropic plasticity is inherited by the growing phase. However, we clearly put into evidence that the anisotropic pattern of the plastic strain and dislocation densities can play a significant role in the morphological evolution of the precipitates. It is worth stressing that this role strongly relies on the spatial pattern of plastic strain and dislocation densities, and as such depends significantly on their inheritance through the moving interface.

<figure><img src="image_8.png" /><figcaption>Figure 7: Viscoplastic case with $K = 150$. (a) Plastic strain $\epsilon_{22}^p$ at 3 h. Dashed lines are the level sets $\varphi = 0.5$ at 0.1 and 3 h; (b) Profiles of $\epsilon_{22}^p$ along the light-blue line in (a) at 1.5 h (dashed line) and at 3 h (thick blue line). Grey slabs indicate the positions of the diffuse interface at 0.1 and 3 h from left to right respectively.</figcaption></figure>

In summary, we have analyzed the influence of the anisotropic viscoplastic driving force accounting for the back stress on the diffusive growth of a precipitate in the presence of coherency stresses. The generated plastic deformation alters the microstructural evolution compared to the elastic case. A morphological transition is observed e.g. from circular to square in the presence of plasticity. When compared to the purely elastic case, this shape change results from an increase of the growth rate along the cubic direction and a decrease along the (110) directions. This result extends the previous results based on isotropic plasticity models stating that, depending on the precipitate geometry, the diffusive growth may either be increased (planar interface [57]) or decreased (sphere [25, 56]).
4.2. Rafting of Ni-based superalloys during creep

Next, we consider another situation exhibiting a strong coupling between phase transformation and plasticity, where an external load drives the plastic flow in the matrix as well as the evolution of the two-phase microstructure: the rafting (directional coarsening) of $\gamma'$ precipitates in Ni-based superalloys under creep loading (Fig. 1), where the anisotropic plastic flow is likely to interfere with the anisotropic microstructure evolution. It is worth stressing that using the present model constitutes a significant improvement over previous works on this phenomenon [26, 58] by accounting simultaneously for all together the ingredients illustrated in the previous sections and based on a better physical description of plasticity.

<figure><img src="image_9.png" /><figcaption>Figure 8: Snapshots of the microstructure (colors indicate $\gamma'$ variants) in a 2.3 × 2.3 µm² periodic box at $t = 0, 1, 5, 15$ and 45 h (from left to right) in the elastic (top) and elasto-viscoplastic (bottom) cases. The arrows point out different types of coarsening mechanisms (see text for details).</figcaption></figure>

First, the cuboidal microstructure is computed considering only elasticity with stress-free conditions, using the parameters in Table 2. The 2.3 × 2.3 µm² periodic box is discretized with 1024² nodes. It has been shown that such discretisation is sufficient to correctly describe the microstructure evolution [59, 46]. Starting from an initial disordered $\gamma$ phase, the alloy is aged during 16 h at 950°C. The resulting microstructure (Fig. 8 - left) consists of $\gamma'$ cuboidal precipitates, aligned along the cubic directions, separated by well-defined $\gamma$ channels (dark regions). The colors of the precipitates identify the four translational variants of the $\gamma'$ structure. The volume fraction of the $\gamma'$ precipitates is 0.63 and their average size is 160 nm. The $\gamma$ channels are 78 nm wide on average. This microstructure is in good qualitative agreement with microstructures observed in AM1 (Fig. 1), even if the precipitate size is smaller than in the standard heating procedure. In the second step, a constant uniaxial stress $\sigma^a=150$ MPa is applied along the [100] direction while holding the temperature at 950°C (Fig. 8). In these conditions, plasticity is only active inside the $\gamma$ channels (using the data in Tab. 3), whereas the $\gamma'$ phase behaves elastically. In order to highlight the influence of the anisotropic plastic driving force on the microstructure evolution, two different creep simulations have been performed: a purely elastic case (Fig. 8 - top) and an elasto-viscoplastic case (Fig. 8 - bottom).

In both cases, the precipitates self-organize into rafts that are aligned along the direction perpendicular to the tensile axis. This is in agreement with previous studies of creep showing that, for a [100] loading axis, the elastic and plastic driving forces lead to the same raft orientation [10, 60, 26]. The comparison of elastic and plastic simulations reveals that viscoplasticity has consequences at different scales. First, it changes the kinetics of rafting and the final shape of the rafts. Indeed, although the rafts are already formed after 4h with similar wavy shapes in both simulations, the microstructure continues to evolve towards straight rafts at 45 h in the elastic case, whereas the wavy rafts are almost frozen after 5h when viscoplasticity is activated in the $\gamma$ phase. This last observation is in agreement with the very slow evolution of the rafts observed in this alloy during the secondary creep regime [61]. It is also consistent with previous works based
on simplified plasticity models [28, 26]. Second, viscoplasticity changes the shape of the precipitates which become more curvy. This indicates that the anisotropic coherency stresses at the origin of the cuboidal shapes are partially relaxed. Inside a raft, at the junctions between two variants, the $γ/γ'$ interfaces are almost aligned with [110] directions (Fig. 8 bottom right, zoom). Moreover, it can be observed that plasticity may induce different behaviors depending on the local configuration. Most of the time, viscoplasticity in the channel slows down the coagulation process of neighboring identical variants (red arrows in Fig. 8). This increases the number of small precipitates when compared to the elastic simulation. However, in a few local configurations, coagulation events are observed to be either accelerated (yellow arrows) or qualitatively modified with the formation of large domains (white arrows).

As expected, the dislocation densities are only positive inside the $γ$ channels with typical spatial distributions illustrated in Figure 9. It appears that the slip systems A6, D6, B5 and C5 are active only in the $γ$ channels perpendicular to the loading axis, whereas A3, C3, D4 and B4 are active in all $γ$ channels. On the contrary, A2, B2, C1 and D1 are not activated during the simulation. To understand these observations, the activation of the glide systems can be analyzed by assessing the elastic stresses in the cuboidal microstructure. Indeed, when the $γ$ channels are very thin because of the high volume fraction of precipitates, the stress tensors $\underline{\sigma}^{\perp}$ and $\underline{\sigma}^{\parallel}$ inside the channels, respectively perpendicular and parallel to the loading axis, can be assumed constant and expressed as:

$$ \underline{\sigma}^{\perp} = \begin{pmatrix} \sigma_1 + \sigma^a & 0 & 0 \\ 0 & \sigma_2 & 0 \\ 0 & 0 & \sigma_2 \end{pmatrix}; \quad \underline{\sigma}^{\parallel} = \begin{pmatrix} \sigma_2 + \sigma^a & 0 & 0 \\ 0 & \sigma_1 & 0 \\ 0 & 0 & \sigma_2 \end{pmatrix} \qquad (23) $$

A rough estimate from the calculations gives $\sigma_1 \approx 0$ and $\sigma_2 \approx -50$ MPa. In the channels perpendicular to the loading axis, the primary slip systems are B4, A3, D4 and C3 as well as A6, D6, B5, and C5. They are associated with the resolved shear stress $(\sigma^a + \sigma_1 - \sigma_2)/\sqrt{6}$. It can be reminded that these glide systems are also those activated during a uniaxial loading of an homogeneous system (§ 3.1). On the contrary, the resolved shear stress in B2, D1, A2 and C1 is zero, explaining why the corresponding dislocation densities remain null inside the channels. In the channels parallel to the loading axis, the primary slip systems B4, A3, D4 and C3 are associated with the resolved shear stress $\sigma^a/\sqrt{6}$. The resolved shear stress in the other slip systems are lower: $(\sigma^a - \sigma_1 + \sigma_2)/\sqrt{6}$ and $(\sigma_2 - \sigma_1)/\sqrt{6}$ for A6, D6, B5 and C5, as well as B2, D1, A2 and C1, respectively. Based on this simple analysis, the slip systems expected to be activated in the two kinds of $γ$ channel are sketched on the right hand side of Figure 9. This analysis qualitatively provides a simple explanation for the distributions featured by the different dislocation densities.

The analysis can be pursued further by comparing the spatial distributions of the dislocation densities in the activated collinear slip systems. The spatial distributions in A3 and C3 systems differ strongly for two reasons. First, near the precipitate corners the resolved shear stresses differ because of the non zero value of the $\sigma_{12}$ component. Second, the collinear interactions in the hardening matrix hinder the activity in a given system where the dislocation density in its collinear system is high. As a consequence, the spatial distributions in A3 and C3 are almost mutually exclusive. On the contrary, the collinear systems A6 and D6 are similarly populated inside the $γ$ channels perpendicular to the loading axis. Indeed, their reduced shear stresses are equal because $\sigma_{13}$ and $\sigma_{23}$ are null due to the 2D geometry of the simulation. Our predictions concerning the distribution of the activated slip systems among the different channels are in agreement with recent calculations using dislocation dynamics [62], although a few small scale features have not been captured, such as the formation of dislocations networks at the interfaces. However, the present model offers two major advantages over dislocation dynamics. First, the microstructure is not static but evolving as a response to the applied stress and plastic deformation. Second, we are able to address realistic strain rates, a crucial point when considering a microstructure evolution resulting from the coupling between diffusion and plasticity. Hence, the present model will be used to address issues still pending such as creep loadings along [111] or [110] directions for which small misorientations can change the raft direction [58, 63].
<figure><img src="image_10.png" /><figcaption>Figure 9: Left: dislocation densities in the slip systems A3, C3, A6 and D6 at 5 h. Right: sketch of the predicted active slip systems in the $\gamma$ channels (see text for details).</figcaption></figure>

5. Conclusion

We have proposed a Phase Field model coupled to a dislocation density based crystal plasticity model. The continuous anisotropic viscoplastic model accounts for several important aspects of plasticity in heterogeneous materials. First, the plastic behavior depends locally on the phase. Second, the forest hardening resulting from the short-range interactions between dislocations is included through a hardening matrix in the Taylor relation whose components are identified on dislocation dynamics calculations taken from the literature. Third, the size effect of the plastic behavior, i.e. the hardening (resp. softening) induced by the decrease (resp. increase) of the size of the plastic regions is also included in the model. This effect results from a back-stress generated by gradients of GND densities. This effect may be rather important in evolving heterogeneous materials at the mesoscale because the sizes and shapes of the coexisting plastic regions may differ and, in addition, evolve during thermo-mechanical treatments.

Static configurations have been considered to emphasize that the model accounts for the forest hardening and plastic size effect. Then, two dimensional simulations have been performed to investigate the influence of plastic deformation on the microstructural evolution. It has been clearly shown that, during its diffusion controlled growth, a misfitting precipitate may undergo a morphological transition induced by the crystal viscoplasticity, as already noted in previous works [56, 55], resulting in an anisotropic shape related to the anisotropy of plastic slip.

Finally, the model has been used to address the microstructural evolution in Ni-based superalloys in creep conditions during which plasticity only proceeds inside the $\gamma$ phase. The parameters of the model have been carefully selected to mimic the behavior of the AM1 monocrystalline superalloy. Most of the important physical phenomena, such as the presence of different crystallographic $\gamma'$ variants, their misfit with the $\gamma$
matrix, the elastic inhomogeneity and anisotropy, the hardening, anisotropy and viscosity of plasticity are included in the model. In agreement with experiments, the model predicts that rafting proceeds perpendicularly to the [100] tensile loading axis and it is shown that plasticity slows down significantly the evolution of the rafts. It is also shown that, depending on the local precipitate behavior, viscoplasticity may either speed up, slow down, or even qualitatively change the coagulation processes that contribute to rafting.

Concerning Ni-based superalloys, the present model could be used to analyse the microstructure evolution during creep loadings along complex axis such as [111] or [110], in particular the origin of the sensitivity of the microstructure evolution to small misorientations of the loading axis [63, 58]. More generally, the work presented in this paper provides a general framework to analyse microstructure evolutions in other alloys in which microstructure evolution and plasticity are strongly coupled.

## 6. Acknowledgments

The authors acknowledge the financial support of the French Agence Nationale de la Recherche (ANR) under reference ANR-BLAN08-1_321567 (project Couphin). The authors would like to thank E. Berthier, A. Gaubert and D. Boivin (ONERA) for providing the SEM images.

## 7. Annex: fcc slip systems

The reduced stress as well as the associated plastic strain are detailed in Table 4 for each glide system of the fcc lattice using the Schmidt and Boas notations [64].

## 8. References

[1] T. Khan, P. Caron, Effect of heat treatment on the creep behavior of a Ni-base single crystal superalloy, in: J. Bilde-Sorensen, N. Hansen, A. Horsewell, T. Leffers, H. Lilholt (Eds.), Deformation of Multi-Phase and Particle Containing Materials, Riso National Laboratory, 1983, pp. 333-338.

[2] J. Ganghoffer, K. Simonsson, S. Denis, E. Gautier, S. Sjöström, A. Simon, Martensitic transformation plasticity simulations by finite elements, Journal de Physique IV 04 (C3) (1994) C3-215-C3-220.

[3] Y. Wen, Y., S. Denis, E. Gautier, Computer simulation of martensitic transformation under stress, J. Phys. IV France 06 (C1) (1996) C1-475-C1-483.

[4] J. Ganghoffer, S. Denis, E. Gautier, S. Sjöström, Micromechanical analysis of the finite element calculation of a diffusional transformation, Journal of Materials Science 32 (18) (1997) 4941-4955.

[5] T. J. Su, E. Aeby-Gautier, S. Denis, Morphology changes in bainite formed under stress, Scripta Materialia 54 (12) (2006) 2185-2189.

[6] H. Nishimori, A. Onuki, Pattern formation in phase-separating alloys with cubic symmetry, Physical Review B 42 (1) (1990) 980-983.

[7] J. Z. Zhu, T. Wang, A. J. Ardell, S. H. Zhou, Z. K. Liu, L. Q. Chen, Three-dimensional phase-field simulations of coarsening kinetics of γ' particles in binary Ni-Al alloys, Acta Materialia 52 (9) (2004) 2837-2845.

[8] T. Wang, G. Sheng, Z. K. Liu, L. Q. Chen, Coarsening kinetics of γ' precipitates in the Ni-Al-Mo system, Acta Materialia 56 (19) (2008) 5544-5551.

[9] G. Boussinot, A. Finel, Y. Le Bouar, Phase-field modeling of bimodal microstructures in nickel-based superalloys, Acta Materialia 57 (2009) 921.

[10] G. Boussinot, Y. Le Bouar, A. Finel, Phase-field simulations with inhomogeneous elasticity comparison with an atomic-scale method and application to superalloys, Acta Materialia 58 (2010) 4170-4181.

[11] Y. Wang, D. Banerjee, C. C. Su, A. G. Khachaturyan, Field kinetic model and computer simulation of precipitation of L1₂ ordered intermetallics from f.c.c. solid solution, Acta Materialia 46 (9) (1998) 2983-3001.

[12] V. Vaithyanathan, L. Q. Chen, Coarsening of ordered intermetallic precipitates with coherency stress, Acta Materialia 50 (2002) 4061-4073.

[13] D. Rodney, A. Finel, Phase field methods and dislocations, MRS Online Proceedings Library Archive, Volume 652, (2000) Y4.9, Published online by Cambridge University Press.

[14] Y. U. Wang, Y. M. Jin, A. M. Cuitino, A. G. Khachaturyan, Nanoscale phase field microelasticity theory of dislocations model and 3D simulations, Acta Materialia 49 (2001) 1847-1857.

[15] M. Koslowski, A. M. Cuitino, M. Ortiz, A phase-field theory of dislocation dynamics, strain hardening and hysteresis in ductile single crystals, Journal of the Mechanics and Physics of Solids 50 (2002) 2597.
[16] D. Rodney, Y. Le Bouar, A. Finel, Phase field methods and dislocations, Acta Materialia 51 (2003) 17–30.

[17] P.-A. Geslin, B. Appolaire, A. Finel, A phase field model for dislocation climb, Appl. Phys. Lett. 104 (2014) 011903.

[18] P.-A. Geslin, B. Appolaire, A. Finel, Multiscale Theory of Dislocation Climb, Physical Review Letters 115 (2015) 265501.

[19] X. H. Guo, S.-Q. Shi, Elastoplastic phase field model for microstructure evolution, Applied Physics Letters 87 (2005) 221910.

[20] R. L. J. M. Ubachs, Thermomechanical modelling of microstructure evolution in solder alloys, Ph.D. thesis, Eindhoven University (2005).

[21] T. Uehara, T. Tsujino, N. Ohno, Elasto-plastic simulation of stress evolution during grain growth using a phase field model, Journal of Crystal Growth (300) (2007) 530–537.

[22] A. Yamanaka, T. Takaki, Y. Tomita, Elastoplastic phase-field simulation of self- and plastic accommodations in cubic → tetragonal martensitic transformation, Materials Science and Engineering A 491 (2008) 378–384.

[23] A. Gaubert, A. Finel, Y. Le Bouar, G. Boussinot, Viscoplastic phase field modelling of rafting in Ni base superalloys, in: D. Jeulin, S. Forest (Eds.), Continuum Models and Discrete Systems CMDS11, Mines Paris Les Presses, 2008, pp. 161–166.

[24] K. Ammar, B. Appolaire, G. Cailletaud, S. Forest, Combining phase field approach and homogenization methods for modelling phase transformation in elastoplastic media, European Journal of Computational Mechanics 18 (2009) 485–523.

[25] K. Ammar, B. Appolaire, G. Cailletaud, S. Forest, Phase field modeling of elasto-plastic deformation induced by diffusion controlled growth of a misfitting spherical precipitate, Philosophical Magazine Letters 91 (2011) 164–172.

[26] M. Cottura, Y. Le Bouar, A. Finel, B. Appolaire, K. Ammar, S. Forest, A phase field model incorporating strain gradient viscoplasticity: Application to rafting in Ni-based superalloys, Journal of the Mechanics and Physics of Solids 60 (7) (2012) 1243–1256.

[27] N. Zhou, C. Shen, M. Mills, Y. Wang, Large-scale three-dimensional phase field simulation of $\gamma'$ rafting and creep deformation, Philosophical Magazine B. 90 (1-4) (2010).

[28] A. Gaubert, Y. Le Bouar, A. Finel, Coupling phase field and viscoplasticity to study rafting in Ni-based superalloys, Philosophical Magazine B 90 (1-4) (2010) 375–404.

[29] R. Madec, B. Devincre, L. Kubin, T. Hoc, D. Rodney, The role of collinear interaction in dislocation-induced hardening, Science 301 (5641) (2003) 1879–1882.

[30] L. P. Evers, W. A. M. Brekelmans, M. G. D. Geers, Non-local crystal plasticity model with intrinsic SSD and GND effects, Journal of the Mechanics and Physics of Solids 52 (2004) 2379–2401.

[31] A. G. Khachaturyan, The Theory of Structural Transformations in Solids, Wiley, New York, 1983.

[32] P. Franciosi, M. Berveiller, A. Zaoui, Latent hardening in copper and aluminium single crystals, Acta Metallurgica 28 (1980) 273–283.

[33] B. Devincre, L. Kubin, T. Hoc, Physical analyses of crystal plasticity by D.D. simulations, Scripta Materialia 54 (2006) 741–746.

[34] A. Vattré, Durcissement des superalliages monocristallins des mécanismes physiques à la modélisation continue, Ph.D. thesis, École Nationale Supérieure des Mines de Paris (2009).

[35] H. Mecking, U. F. Kocks, Kinetics of flow and strain-hardening, Acta Met., 29 (1981) 1865–1875.

[36] U. F. Kocks, H. Mecking, Physics and phenomenology of strain hardening the fcc case, Progress in Materials Science 48 (2003) 171–273.

[37] C. Teodosiu, J. L. Raphanel, Finite element simulation of the large elastoplastic deformation, in: C. Teodosiu, J. L. Raphanel, F. Sidoroff (Eds.), Large Plastic Deformations - Fundamentals Aspects and Applications to Metal Forming, A. A. Balkema, 1993, pp. 153–168.

[38] L. Kubin, B. Devincre, T. Hoc, Modeling dislocation storage rates and mean free paths in face-centered cubic crystals, Acta Mat. 56 (2008) 6040–6049.

[39] L. Nemat-Nasser, S. Ni, T. Okinaka, A constitutive model for fcc crystals with application to polycrystalline OFHC copper, Mechanics of Materials 30 (4) (1998) 325–341.

[40] I. Groma, P. Balogh, Link between the individual and continuum approaches of the description of the collective behavior of dislocations, Materials Science and Engineering A 234–236 (1997) 249–252.

[41] I. Groma, F. F. Csikor, M. Zaiser, Spatial correlations and higher-order gradient terms in a continuum description of dislocation dynamics, Acta Materialia 51 (5) (2003) 1271–1281.

[42] P. L. Valdenaire, Y. Le Bouar, B. Appolaire, A. Finel, Density-based crystal plasticity: from the discrete to the continuum http://arxiv.org/abs/1602.08695 Submitted to Physical Review B.

[43] A. J. Ardell, V. Ozolins, Trans-interface diffusion-controlled coarsening, Nature Materials 4 (2005) 309–316.

[44] A. Royer, P. Bastie, M. Véron, In situ determination of $\gamma'$ phase volume fraction and of relations between lattice parameters and precipitate morphology in Ni-based single crystal superalloy, Acta Materialia 46 (1998) 5357–5368.

[45] F. Diologent, P. Caron, T. d'Almeida, A. Jacques, P. Bastie, The $\gamma - \gamma'$ mismatch in Ni based superalloys in situ measurements during a creep test, Nuclear Instruments and Methods in Physics Research B 200 (2003) 346–351.

[46] M. Cottura, B. Appolaire, Y. Le Bouar, A. Finel, Rôle of elastic inhomogeneity in the development of cuboidal microstructures in Ni-based superalloys, Acta Mater. 94 (2015) 15–25.

[47] F. Fujiwara, M. Watanabe, Z. H. N. Nemoto, K. Noumi, T. Simozaki, Measurements of intrinsic diffusion coefficients in Ni₃Al, in: K. O. M. Koiwa, T. Miyazaki (Eds.), International Conference on Solid-Solid Transformations, in JIMIC-3, 1999, pp. 481–484.

[48] L. Espié, Etude expérimentale et modélisation numérique du comportement de monocristaux de superalliages, PhD Thesis, École Nationale Supérieure des Mines de Paris (1996).

[49] A. Fredholm, Monocristaux d'alliages base nickel relation entre composition, microstructure et comportement en fluage à
8 REFERENCES

haute température, Ph.D. thesis, École Nationale Supérieure des Mines de Paris (1987).

[50] D. Ayrault, Fluage à haute température de superalliages base nickel monocristallins, Ph.D. thesis, École Nationale Supérieure des Mines de Paris (1989).

[51] B. Devincre, T. Hoc, L. P. Kubin, Collinear interactions of dislocations and slip systems, Materials Science and Engineering A 400-401 (0) (2005) 182-185.

[52] N. M. Cordero, A. Gaubert, S. Forest, E. Busso, F. Gallerneau, S. Kruch, Size effects in generalised continuum crystal plasticity for two-phase laminates, Journal of the Mechanics and Physics of Solids 58 (2010) 1963-1994.

[53] X.G. Wang, J.L. Liu, T. Jin, X.F. Sun, Z.Q. Hu, J.H. Do, B.G. Choi, I.S. Kim, C.Y. Jo, Dislocation motion during high-temperature low-stress creep in Ru-free and Ru-containing single-crystal superalloys, Materials and Design 67 (2015) 543-55.

[54] M. Cottura, B. Appolaire, A. Finel, Y. Le Bouar, Plastic relaxation during growth of Widmanstätten plates, Scripta Mat. 108 (2015) 117-121.

[55] V. de Rancourt, K. Ammar, B. Appolaire, S. Forest, Homogenization of viscoplastic constitutive laws within a phase field approach, Journal of the Mechanics and Physics of Solids (2016) in Press.

[56] K. Ammar, B. Appolaire, S. Forest, M. Cottura and Y. Le Bouar and A. Finel, Modelling inheritance of plastic deformation during migration of phase boundaries using a phase field method, Meccanica 49 (2014) 2699-2717.

[57] K. Ammar, Modélisation et simulation du couplage changement de phases-mécanique par la méthode des champs de phases, Ph.D. thesis, École nationale supérieure des mines de Paris (2010).

[58] A. Gaubert, M. Jouiad and J. Cormier and Y. Le Bouar, J. Ghighi, 3-D imaging and phase field simulations of the microstructure evolution during creep tests of < 011 > oriented Ni-based superalloys, Acta Mater. 84 (2015) 237-255.

[59] M. Cottura, Modélisation Champ de Phase du couplage entre evolution microstructurale et propriétés mécaniques: Cube, radeaux et pointes, PhD thesis, Université Paris VI (2013).

[60] N. Zhou, C. Shen, M. Mills, Y. Wang, Contributions from elastic inhomogeneity and from plasticity to γ' rafting in single-crystal Ni-Al, Acta Materialia 56 (2008) 6156-6173.

[61] F. Diologent, P. Caron, On the creep behavior at 1033 K of new generation single-crystal superalloys, Materials Science and Engineering A 385 (2004) 245-257.

[62] A. Vattré, B. Devincre, A. Roos, Orientation dependence of plastic deformation in nickel-based single crystal superalloys: Discrete-continuous model simulations, Acta Mater. 58 (2010) 1938-1951.

[63] L. Agudo Jacome, P. Nortershuser, J.-K. Heyer, A. Lahni, J. Frenzel, A. Dlouhy, C. Somsen, G. Eggeler, High-temperature and low-stress creep anisotropy of single-crystal superalloys, Acta Materialia 61 (8) (2013) 2926-2943.

[64] L.P. Kubin, Dislocations, Mesoscale Simulations and Plastic Flow, Oxford Series On Materials Modelling, Oxford University Press 2013.
<table><caption>Table 4: Reduced stress $\sigma_r$ and plastic strain $\epsilon^p$ associated to the glide of the fcc slip systems. $\mathbf{m}$ and $\mathbf{n}$ are the slip direction and the normal to the glide plane, respectively.</caption><thead><tr><th>System</th><th>√2 <b>m</b></th><th>√3 <b>n</b></th><th>√6 σ<sub>r</sub></th><th>√6 ε<sup>p</sup></th></tr></thead><tbody><tr><td>A2</td><td>[011]</td><td>[111]</td><td>σ<sub>33</sub> - σ<sub>22</sub> + σ<sub>12</sub> - σ<sub>13</sub></td><td>( 0 1/2 1/2<br>1/2 1 0<br>1/2 0 1 )</td></tr><tr><td>A3</td><td>[101]</td><td>[111]</td><td>σ<sub>33</sub> - σ<sub>11</sub> + σ<sub>12</sub> + σ<sub>23</sub></td><td>( 1 1/2 0<br>1/2 0 1/2<br>0 1/2 1 )</td></tr><tr><td>A6</td><td>[110]</td><td>[111]</td><td>σ<sub>22</sub> - σ<sub>11</sub> + σ<sub>13</sub> + σ<sub>23</sub></td><td>( 1 0 1/2<br>0 1 1/2<br>1/2 1/2 0 )</td></tr><tr><td>B2</td><td>[011]</td><td>[111]</td><td>σ<sub>33</sub> - σ<sub>22</sub> - σ<sub>12</sub> + σ<sub>13</sub></td><td>( 0 1/2 1/2<br>1/2 1 0<br>1/2 0 1 )</td></tr><tr><td>B4</td><td>[101]</td><td>[111]</td><td>σ<sub>33</sub> - σ<sub>11</sub> - σ<sub>12</sub> + σ<sub>23</sub></td><td>( 1 1/2 0<br>1/2 0 1/2<br>0 1/2 1 )</td></tr><tr><td>B5</td><td>[110]</td><td>[111]</td><td>σ<sub>22</sub> - σ<sub>11</sub> - σ<sub>13</sub> + σ<sub>23</sub></td><td>( 1 0 1/2<br>0 1 1/2<br>1/2 1/2 0 )</td></tr><tr><td>C1</td><td>[011]</td><td>[111]</td><td>σ<sub>33</sub> - σ<sub>22</sub> - σ<sub>12</sub> - σ<sub>13</sub></td><td>( 0 1/2 1/2<br>1/2 1 0<br>1/2 0 1 )</td></tr><tr><td>C3</td><td>[101]</td><td>[111]</td><td>σ<sub>33</sub> - σ<sub>11</sub> - σ<sub>12</sub> - σ<sub>23</sub></td><td>( 1 1/2 0<br>1/2 0 1/2<br>0 1/2 1 )</td></tr><tr><td>C5</td><td>[110]</td><td>[111]</td><td>σ<sub>11</sub> - σ<sub>22</sub> - σ<sub>13</sub> + σ<sub>23</sub></td><td>( 1 0 1/2<br>0 1 1/2<br>1/2 1/2 0 )</td></tr><tr><td>D1</td><td>[011]</td><td>[111]</td><td>σ<sub>33</sub> - σ<sub>22</sub> + σ<sub>12</sub> + σ<sub>13</sub></td><td>( 0 1/2 1/2<br>1/2 1 0<br>1/2 0 1 )</td></tr><tr><td>D4</td><td>[101]</td><td>[111]</td><td>σ<sub>33</sub> - σ<sub>11</sub> + σ<sub>12</sub> - σ<sub>23</sub></td><td>( 1 1/2 0<br>1/2 0 1/2<br>0 1/2 1 )</td></tr><tr><td>D6</td><td>[110]</td><td>[111]</td><td>σ<sub>11</sub> - σ<sub>22</sub> + σ<sub>13</sub> + σ<sub>23</sub></td><td>( 1 0 1/2<br>0 1 1/2<br>1/2 1/2 0 )</td></tr></tbody></table>