# Thermo-elastic micromechanical modeling of tetragonal ZrO_2 with a herringbone microstructure inherited from the cubic phase

Lukas Petrich, Katell Derrien, Volker Schmidt, René Guinebretière, Hervé Moulinec, Olivier Castelnau

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

## ► To cite this version:

Lukas Petrich, Katell Derrien, Volker Schmidt, René Guinebretière, Hervé Moulinec, et al., Thermo-elastic micromechanical modeling of tetragonal ZrO_2 with a herringbone microstructure inherited from the cubic phase. Materials & Design, 2025, 257, pp.114425. 10.1016/j.matdes.2025.114425. hal-05195270

HAL Id: hal-05195270

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

Submitted on 1 Aug 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_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_3.png" /></figure>

# Thermo-elastic micromechanical modeling of tetragonal ZrO₂ with a herringbone microstructure inherited from the cubic phase

Lukas Petrich<sup>a, ID,*</sup>, Katell Derrien<sup>b, ID</sup>, Volker Schmidt<sup>a</sup>, René Guinebretière<sup>c, ID</sup>, Hervé Moulinec<sup>d</sup>, Olivier Castelnau<sup>b, ID,*</sup>

## ARTICLE INFO

### Keywords:

Zirconia
Three-dimensional herringbone microstructure
Phase transition
Full-field residual stress
Thermo-elasticity
Stochastic geometry

## ABSTRACT

The distribution of residual stresses at the crystal scale is investigated in several random polycrystalline aggregates of zirconia. The underlying stochastic model generates tetragonal crystals forming three-dimensional herringbone microstructures, originating from cooling an initial single cubic crystal at high temperatures followed by a solid-state phase transition. The obtained microstructures are constructed following crystallographic constraints in terms of variant selection, as well as twin and band boundary orientations. This stochastic modeling approach allows generating microstructures with twin domains of various aspect ratios, embedded in several Voronoi cells as observed when the phase transition starts simultaneously from different locations in the cubic crystal. Considering the transformation strain and the anisotropic elastic and dilation properties at the crystal scale, these microstructures are solved with the spectral (FFT) full-field method for pure thermal loading (cooling). Thermal dilation has a limited effect on the residual stress field, about one order of magnitude smaller than the transformation strain, even for a cooling of 1000 °C. Normal stresses along the tetragonal crystal axes are in the order of a few GPa and proportional to the tetragonality of the zirconia crystal, except for a specific microstructure where stresses even vanish for infinite aspect ratios.

## 1. Introduction

Zirconia is a well-known oxide used in numerous applications which is renowned for its dual benefits: excellent thermo-mechanical properties (e.g. used as a refractory material in high-temperature applications) and high ionic conductivity. Dense zirconia-based materials are usually obtained through a sintering process at high temperatures or fuse cast from the liquid state. In both cases, the final polycrystalline material is out-of-equilibrium, with the presence of a dense network of nanocracks [4,10,30]. Although crack formation induces local decreases in the stress level, residual mean stresses in the range of a few tens to a few hundreds of MPa are most often present [33]. Such evaluations of mean stresses through conventional x-ray or neutron diffraction experiments give values that are averaged over a very large number of crystals that are in the probed volume. Raman spectroscopy can detect stresses up to 1 GPa [7]. Our recent experimental work performed on pure zirconia polycrystals and based on Laue microdiffraction demonstrates that the

local residual stresses lie in the GPa range [32]. In this last case, the size of the probed volume for each measurement was of a few hundreds of nanometers, whereas it was micrometric for the Raman measurements and millimetric for conventional x-rays or neutron experiments.

The coupling between solid-state phase transitions (SPT), residual stresses and mechanical response is complex. Nevertheless, for large variations of temperature, the mechanical response is mainly guided by the formation of local heterogeneous stress fields, due to the intrinsic highly anisotropic elastic and thermal expansion properties of the crystals constituting the material, the transformation strain arising during SPT, and mechanical interactions between neighboring twin domains (so-called “microstructural effects”) in order to accommodate the incompatible thermal and transformation deformations.

Due to the ability of zirconium to accommodate different atomic coordinations, zirconia crystallizes in different phases depending on pressure and temperature. The different phase stability domains can

<footnote><sup>a</sup> Institute of Stochastics, Faculty of Mathematics and Economics, Ulm University, Helmholtzstraße 18, Ulm, 89069, Germany</footnote>

<footnote><sup>b</sup> PIMM, Arts et Métiers Institute of Technology, CNRS, Cnam, 151 boulevard de l'hopital, Paris, 75013, France</footnote>

<footnote><sup>c</sup> IRCER, Université de Limoges, CNRS, Centre Européen de la Céramique, Limoges, 87068, France</footnote>

<footnote><sup>d</sup> Aix-Marseille Université, CNRS, Centrale Marseille, LMA UMR 7031, 4 impasse Nikola Tesla, Marseille, 13453, France</footnote>

<footnote>* Corresponding authors.
E-mail addresses: lukas.petrich@alumni.uni-ulm.de (L. Petrich), olivier.castelnau@ensam.eu (O. Castelnau).</footnote>
vary with the nature and the amount of doping elements [18], but the present study is restricted to pure zirconia (ZrO₂). At atmospheric pressure, stress-free pure zirconia solidifies into a cubic structure (denoted hereafter by ‘c’, space group *Fm$\overline{3}$m) at about 2700°C, transforms to tetragonal one (denoted by ‘t’, space group *P4$_2$/nmc) upon cooling to 2300°C and becomes monoclinic (denoted by ‘m’, space group *P2$_1$/c) at 1170°C. In a pure zirconia single crystal, the *t* → *m* SPT is of first order and it is associated with a large volume expansion of about 4% [47]. It has been described in detail in, e.g., Simha [46] and Kelly and Francis Rose [18]. In polycrystals, this *t* → *m* transformation is at the origin of huge internal stresses (GPa range) promoting the development of a microcrack network into bulk specimens [14,31,32]. On the contrary, the *c* → *t* SPT is a continuous second-order phase transition that can be described by an alternative positive and negative shift of the oxygen atoms along the *c*-axis of the cubic cell. It means that along this axis, the oxygen coordinate is no longer than 0.25 but equal to 0.25 – *η*, where *η* is the order parameter in the Landau sense [21,50]. Because the *c* → *t* SPT corresponds to a decrease of symmetry, a single cubic crystal induces the formation of several adjacent crystals (twins) upon transition. In Section 2, we describe the microstructures inherited from this SPT process according to these crystallography rules.

The *c* → *t* SPT induces an anisotropic stress-free strain. Further cooling of *t* crystals leads to anisotropic dilation. These two incompatible strains must be accommodated by a significant elastic strain (in the order of 10⁻²) and huge stress (in the GPa range) resulting from the anisotropic elastic stiffness at the crystal scale. The aim of the present paper is to make the link between the stress distribution at the crystal scale and the material microstructure, accounting for the crystallographic constraints of the *c* → *t* SPT.

For the sake of modeling the mechanical response of a polycrystalline aggregate using a full-field numerical scheme, representative microstructures are required. One possibility for this are phase field simulations (see, e.g., [27,36,37,49]), which stand as a powerful computational tool for capturing the complex interplay of thermodynamics and kinetics. Applications of modeling the evolution of the crack network and its influence on the zirconia transformation essentially focus on the *t* → *m* SPT, e.g. Mamivand et al. [26], Zhu and Luo [52]. We are not aware of any publications dealing with the phase-field modeling of the *c* → *t* SPT. However, these simulations usually require a lot of computing power and time due to their inherent complexity and the need for fine spatial-temporal resolution. In pursuit of a more expedient microstructure generation for our statistical analysis, our research adopts more direct and faster modeling (a few minutes of runtime compared to typically several hours). Indeed, another well-established approach for generating virtual, yet realistic, microstructures of polycrystalline materials is given by using tessellations, i.e. a geometric partitioning of a given sampling window into individual cells (i.e. twin domains) without gaps. These can be obtained by employing methods from spatial stochastics such as random point processes or random tessellations, see e.g. [11,12,43,44]. Alternatively, it is possible to start with a random, but over-simplistic, microstructure, e.g., from a packing algorithm, and modify it with an optimization procedure that utilizes some experimentally observed characteristic, in order to acquire a more realistic structure, see Quey and Renversade [42], Serrao et al. [45]. However, all these approaches rely on three-dimensional experimental data, which is not readily available for the cubic-tetragonal phase transition of zirconia due to the high temperatures involved. In the present paper, we therefore employ crystallographic laws associated with the phase transition process to obtain realistic microstructures without the need of experimental data.

The structure of the present paper is as follows: In Section 2 we provide a detailed description of the herringbone microstructure of *t* ZrO₂, which originated from a *c* ZrO₂ single crystal and was subjected to the SPT. Then, in Section 3, we present the microstructure model that has been developed to generate 3D herringbone microstructures of the *t* zirconia phase. Details of the numerical implementation of the spectral

method used to compute the thermo-mechanical response of the microstructures are presented in Section 4. Full-field results are shown in Section 5 and compared with the analytical solution for lamella with infinite aspect ratio recalled in Appendix A. The paper ends with some concluding remarks in Section 6.

## 2. Herringbone microstructure inherited from *c* → *t* phase transformation

In this section, we provide the crystallographic elements necessary for the random generation of synthetic herringbone microstructures of *t* zirconia (Section 3), as the 3D geometry of crystal orientations and shapes is strongly related to crystallographic relationships. This study is restricted to the case of the SPT occurring during the cooling of a *c* single crystal, ending with a herringbone structure made of many *t* crystals.

In the following, by $(hk\ell)_c$ and $(hk\ell)_t$ we denote the Miller indexes of lattice planes when expressed with respect to the *c* and *t* lattices, respectively. A similar convention is used for directions, i.e. $[uvw]_c$ and $[uvw]_t$. As mentioned above the tetragonal lattice is a primitive one and the $a_t$ and $b_t$ lattice vectors are rotated by 45° with respect to the corresponding cubic lattice vectors. In order to avoid additional rotations, all along the present article, the *t* structure will be described with respect to the multiple “tetragonal face centered” cell where the axes of the cubic and tetragonal cells are in parallel of each other. The obtained tetragonal microstructure results from the association of all the rotational variants allowed by the *c* and *t* point group relationships. It is noteworthy that the 4/mmm point group is a subgroup of the $m\overline{3}m$ one [2]. These two point groups do not belong to the same crystal system, but instead are the holohedral class of their crystal lattice. This transition is thus a ferroelastic one, i.e. it is associated to a spontaneous change of the crystal lattice and of a variation of the number of independent nonzero components of the symmetric second rank tensor associated to spontaneous strain. According to Aizu [1] it can be noted by $m\overline{3}mF4/mmm$ and the number of “spontaneous nonzero strain components” (called “stress-free strain components” in the following) is equal to 3 in two cases, while the number of those that are independent is 1 and 2 for the $m\overline{3}m$ and 4/mmm, respectively. The number of variants and their relative orientations are directly induced by the loss of symmetry associated with the transition from the $m\overline{3}m$ to 4/mmm point group. It is obvious that this is associated with the loss of the 3-fold axis and thus it generates three single crystalline states (variants). In the cubic state, the 3-fold axis is around the [111] direction (or any equivalent direction), that is collinear to the [111]* direction (here the notation ‘.*’ stands for the reciprocal space). This reciprocal space direction is orthogonal to any $\langle h\overline{h}0 \rangle^*$ directions. Thus, the 3-fold axis around the [111] direction lies in the $\{h\overline{h}0\}$ plane families that are oriented at 120° from each other. Considering that the $c_t/a_t$ ratio in the tetragonal state is close to 1, the twin domains appearing through the $c_t$ → *t* SPT are separated by twin walls belonging to a {110}, {101}, or {011} planes. These planes are rotated by 120° with each other around [111] that is the common zone axis.

From the experimental point of view, the formation of *t* microstructures in doped zirconia has been described in detail by Hayakawa et al. [16]. Typical examples of herringbone microstructures are shown also in Hayakawa et al. [15]. There are three trivial possibilities to grow a *t* child phase into a *c* parent single crystal. The two sets of lattice vectors, $a_t, b_t, c_t$ and $a_c, b_c, c_c$, stay parallel with each other, but the $c_t$ vector can align either with $a_c$ or $b_c$ or $c_c$ due to the cubic and tetragonal symmetries. According to the International Table For Crystallography [2], Hayakawa et al. [16] denotes the orientation relationship for which $c_t \parallel a_c$ (with $a_t \parallel b_c$ and $b_t \parallel c_c$) as *x* variant. Similarly, a *y* variant has $c_t \parallel b_c$ (with $a_t \parallel c_c$ and $b_t \parallel a_c$), and a *z* variant has $c_t \parallel c_c$ (with $a_t \parallel a_c$ and $b_t \parallel b_c$). The herringbone structure consists of bands containing either *x* – *y*, *x* – *z* or *y* – *z* crystal sequences.

A schematic representation of such a microstructure is provided in Fig. 1. The *x* and *z* twins in the *x* – *z* band are separated by a (010)$_t$ plane, the *y* and *z* twins in the *y* – *z* band by (011)$_t$, and the two bands
<figure><img src="image_4.png" /><figcaption>Fig. 1. Schematic illustration of the herringbone structure comprising many t twins denoted x, y or z according to the convention indicated in the text. The twin walls (in blue) separate x − z or y − z twin pairs that constitute the two dif- ferent kinds of bands. Band boundaries (in red) separate x − y twins. The figure is a section of the microstructure normal to the [111] direction. It corresponds to the stacking sequence denoted x − z ↔ y − z . The loss of the 3-fold symme- try is clearly illustrated by the association of all twins.</figcaption></figure>

by the (1 ¯ 10)$_{t}$ plane. Note the continuity of z twins across the (1 ¯ 10)$_{t}$ band boundaries, whereas x and y twins alternate. As the three possible vari- ants have an identical probability to grow, due to the symmetries of the c structure, their volume fractions in the final herringbone structure should be 1/3, as also observed experimentally. Therefore, the width of z twins in the configuration (called “stacking sequence” in the follow- ing) x − z ↔ y − z shown in Fig. 1 should be half that of x and y ones as in the figure (and similarly for other stacking sequences).

The plane at the interface between two crystals, called twin wall, can be either (1 ¯ 10)$_{t}$ or (110)$_{t}$ when there is a succession of x and y crystals. Similarly, it is either (101)$_{t}$ or (1 ¯ 10)$_{t}$ for a succession of x and z orientations, and (011)$_{t}$ or (0 ¯ 11)$_{t}$ for a succession of y and z orientations. In general, a single twin wall is selected all along a band, but it can appear that the plane changes within the same band, e.g. the plane is (101)$_{t}$ in some part of the band but (1 ¯ 10)$_{t}$ further away, e.g. see Figure 6 of Hayakawa et al. [16], leading to more complex microstructures that are not considered in the present work.

3. Generation of synthetic 3D herringbone microstructures

In this section, we detail the model that has been developed for the generation of synthetic herringbone microstructures. Due to the lack of experimental data of pure zirconia at the very high temperatures ( ∼ 2300 $^{◦}$C) of the phase transition, it is infeasible to directly calibrate the output of a microstructure generator. As an alternative approach, we use the well-established crystallographic rules described in Section 2. In order to keep the generation code simple, we have slightly relaxed the strict crystallographic rules mentioned above, but rectify this during the later analysis: We still consider that only three t variants can form from the c parent crystal. Sequences of x − y , x − z and y − z twins are considered, organized in bands, but we take more degrees of freedom concerning the twin planes and band boundaries by unfolding the crystal symmetries - e.g. we consider all possible planes families of the { 110 } form.

The 3D t microstructures are generated in two steps. First some seed points are simulated in a sampling window filled with a cubic single crystal at high temperature. Then, the region of influence for each seed point is subdivided by inserting herringbone structures generated in ac- cordance with the crystallographic rules described in Section 2. Fig. 2 illustrates the whole microstructure generation process.

More formally, consider the sampling window W = [0 , 1] 3 contain- ing a c single crystal with lattice parameter a$_{c}$ and whose crystal axes a$_{c}$ , b$_{c}$ and c$_{c}$ are oriented along the axes of the used coordinate sys- tem. Note that the generated microstructures do not have a physical length scale attached, which would be superfluous as the subsequent

mechanical analysis is independent of the length scale. Based on the model parameter n Vo = 1 , 2 , . . . , which specifies the number of nucle- ation (i.e. seed) points for the phase transition, in total n Vo seed points p Vo 1 , . . . , p Vo n ∈ W are drawn uniformly and independently in W (i.e. a realization of a binomial point process, see [6]). With these, a Voronoi tessellation [6] with the cells C Vo 1 , . . . , C Vo n ⊂ W is built. Here, the i -th cell C i Vo comprises all points in the sampling window W that are closer to the i -th seed point p i Vo than to any other seed point. In other words, C i Vo is given by

$$
C_i^{\text{Vo}} = \{ \mathbf{x} \in \mathcal{W} : \|\mathbf{x} - \mathbf{p}_i^{\text{Vo}}\| \le \|\mathbf{x} - \mathbf{p}_j^{\text{Vo}}\| \text{ for all } j \ne i \} \quad (1)
$$

for all i = 1 , . . . , n Vo where ‖·‖ denotes the Euclidean distance. In each Voronoi cell C Vo i , a t herringbone structure is then constructed as fol- lows: As the first step, to disentangle the crystal symmetries, a coordi- nate system is chosen at random from the 24 equivalent possibilities that are aligned with the cubic crystal of the sampling window (6 possibilities to align the new x -axis with the axes of the coordinate system, then 4 for the y -axis and only 1 remaining for the z -axis to ensure a right-handed system). This leads to (orthonormal) basis vectors b$_{1}$ , b$_{2}$ , b$_{3}$ ∈ R 3 and the lattice L Vo = ( { a$_{c}$i$_{1}$b$_{1}$ , a$_{c}$i$_{2}$b$_{2}$ , a$_{c}$i$_{3}$b$_{3}$ } : i$_{1}$, i$_{2}$, i$_{3}$ integers), which is used in the following.

Next, a herringbone configuration as described in Section 2 is se- lected at random consisting of a stacking sequence s (either y − x ↔ z − x , x − y ↔ z − y , or x − z ↔ y − z ), twin wall boundaries and band boundaries. Let n bn be the unit normal vector for the band boundaries and n t w A , n t w B be the unit normal vectors of the twin walls for crystals in the two bands (denoted A and B), respectively. With this, parallel planes P bn K min , . . . , P bn K max ⊂ C Vo that define individual band boundaries are con- structed such that they are orthogonal to n bn and the distances to their nearest neighboring planes are equal to w bn for some width parameter w $^{bn}$> 0 and integers K bn min ≤ K bn max . More specifically, for the k -th band boundary plane, consider the points p bn k ∈ P bn k given by

$$
\mathbf{p}_k^{\text{bn}} = \left[ \frac{k+1}{2} \right] \mathbf{v}_A + \left[ \frac{k}{2} \right] \mathbf{v}_B + \mathbf{p}^{\text{Vo}} \quad (2)
$$

with the shift vector

$$
\mathbf{v}_A = w^{\text{bn}} (\mathbf{n}^{\text{bn}} - \langle \mathbf{n}^{\text{bn}}, \mathbf{n}_A^{\text{t}w} \rangle \mathbf{n}_A^{\text{t}w}) \quad (3)
$$

for the A-bands and

$$
\mathbf{v}_B = w^{\text{bn}} (\mathbf{n}^{\text{bn}} - \langle \mathbf{n}^{\text{bn}}, \mathbf{n}_B^{\text{t}w} \rangle \mathbf{n}_B^{\text{t}w}) \quad (4)
$$

for the B-bands where [ · ] denotes the floor function, 〈· , ·〉 denotes the scalar product and K bn min ≤ k ≤ K bn max . The corresponding band boundary plane P bn k is then given by

$$
P_k^{\text{bn}} = \{ \mathbf{x} \in C^{\text{Vo}} : \langle (\mathbf{x} - \mathbf{p}_k^{\text{bn}}), \mathbf{n}^{\text{bn}} \rangle = 0 \}. \quad (5)
$$

In the following, the band C bn k ⊂ C Vo is the region enclosed by P bn k − 1 and P $^{bn}$. The values of K bn min and K bn max are selected such that C bn k − 1 and C bn K max have a volume larger or equal to some minimum volume v $^{min}$> 0, respectively. In the case where the Voronoi cell C $^{Vo}$has already a smaller volume than 2 v $^{min}$, no band boundary is inserted and C bn 0 = C $^{Vo}$. Otherwise, P $^{bn}$0 is certainly present because the seed point p $^{Vo}$0 = p $^{bn}$0 is always included in the Voronoi cell C $^{Vo}$. As the last step for constructing the bands, a type is assigned to each C $^{bn}$k : if k is odd the C $^{bn}$k is of type A, and of type B otherwise.

Similarly to the Voronoi cells above, every band is subsequently sub- divided into individual t crystals by inserting twin walls P t w k ,L t w min , . . . , P t w k ,L t w max ⊂ C $^{bn}$k , with some integers L t w min ≤ L t w max , whose normal vector is given by the considered herringbone configuration. The twin wall normal vector n t w A if the k -th band is of type A, and n t w B if it is of type B. Each band consists of crystals from two dif-
<figure><img src="image_5.png" /><figcaption>Fig. 2. 2D schematic of the microstructure generation process. (a) draw seed points ($n^{Vo} = 2$) for the phase transformation in the sampling window corresponding to a cubic single crystal; (b) compute the Voronoi tessellation for the seed points and assign each cell a random herringbone configuration; based on this configuration generate band boundaries (c) and twin wall boundaries with corresponding crystal orientations (d).</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Fig. 3. 3D visualizations of two realizations from the dataset $D^{hs}$ comprising a unique Voronoi cell ($n^{Vo} = 1$) with (a) aspect ratio $\alpha_{AR} = 1$ and (b) $\alpha_{AR} = 10$, as well as two realizations drawn from $D^{sm}$ where (c) $n^{Vo} = 4$ with $\alpha_{AR} = 10$, and (d) $n^{Vo} = 20$ with $\alpha_{AR} = 10$.</figcaption></figure>

ferent crystal variants – one that is unique for this type of band, and another one that is present in both types, as described in the previous section. Their assignment to each band is dictated by the stacking sequence $s$ of the herringbone configuration, e.g., for $s = x-z \leftrightarrow y-z$ as in Fig. 1, the $z$-variant is shared between all bands, the $y$-variant belongs to the A-bands and the $z$-variant to the B-bands. The (maximum) width of a crystal is determined by whether it is shared or unique. In the first case, twin walls $P_{k,\ell}^{tw}$ are separated by a distance $u_s^{tw} > 0$ and by $u_u^{tw} > 0$ in the second case. For our simulations, we chose $u_u^{tw} = 2u_s^{tw}$ and $u_s^{tw} = \alpha_{AR}u_s^{bn}$ for some aspect ratio parameter $\alpha_{AR} > 0$. The $\ell$-th twin plane contains the point

$$ \mathbf{p}_{k,\ell}^{\text{tw}} = \left( \left[ \frac{\ell + 1}{2} \right] u_{\text{u}}^{\text{tw}} + \left[ \frac{\ell}{2} \right] u_{\text{s}}^{\text{tw}} \right) \mathbf{n}^{\text{tw}} + \mathbf{p}_{k}^{\text{bn}} \quad (6) $$

and is therefore given by

$$ P_{k,\ell}^{\text{tw}} = \{ \mathbf{x} \in C_k^{\text{bn}} : \langle \mathbf{x} - \mathbf{p}_{k,\ell}^{\text{tw}}, \mathbf{n}^{\text{tw}} \rangle = 0 \} \quad (7) $$

with the integer $L_{\min}^{\text{tw}} \le \ell \le L_{\max}^{\text{tw}}$. The crystal $G_{k,\ell}$ is given by the band boundaries $P_{k-1}^{\text{bn}}$ and $P_k^{\text{bn}}$ as well as the twin walls $P_{k,\ell-1}^{\text{tw}}$ and $P_{k,\ell}^{\text{tw}}$. The values of $L_{\min}^{\text{tw}}$ and $L_{\max}^{\text{tw}}$ are chosen such that $G_{K_{\min}^{\text{bn}}+1}^{\text{tw}}$ and $G_{L_{\max}^{\text{tw}}}^{\text{tw}}$ have a volume larger or equal to the minimum volume $v_{\min}^{\text{min}}$, respectively. For bands whose volume is lower than $2v_{\min}^{\text{min}}$, only one crystal exists, which is equal to the band. A lattice orientation is then assigned to the crystals, depending on the selected variants sequence.

Note that not all combinations of variant sequences, twin walls, and band boundaries lead to geometrically acceptable microstructures. For example, in a hypothetical $x-z \leftrightarrow y-z$ sequence with a (110) band boundaries, a (101) twin wall for the $x-z$ A band and a (011) twin wall for the $y-z$ B band, one ends up with an unrealistic microstructure with no continuity of the $z$ variant as it should be. Such microstructures are thus discarded in the present work.

In the following, two datasets are considered.

i) The first one, denoted by $D^{hs}$, focuses only on the herringbone structure inside a single Voronoi tessellation cell. More precisely,

it was obtained by forcing the Voronoi tessellation to have exactly one cell $C^{Vo} = W$ with seed point $\mathbf{p}^{Vo} = (0, 0, 0)$. This dataset comprises 100 realization for each value of the aspect ratio $\alpha_{AR} \in \{1, 10\}$. The minimum volume parameter was set to $v_{\min}^{\text{min}} = 3 \cdot 10^{-5}$.

ii) The second dataset $D^{sm}$ consists of 400 realizations of the (unmodified) stochastic model. For this, 100 microstructures were obtained for each parameter configuration given by a value for the number of Voronoi cells $n^{Vo} \in \{4, 20\}$ and for the aspect ratio $\alpha_{AR} \in \{1, 10\}$. Again, $v_{\min}^{\text{min}} = 3 \cdot 10^{-5}$ was chosen.

In order to be able to compute the mechanical response of those microstructures with the FFT method (see next section), all realizations were discretized into 256 × 256 × 256 regularly spaced voxels. Fig. 3 shows three typical 3D microstructures obtained with the above model, comprising 1, 4 and 20 Voronoi cells filled with herringbone structures. It can be observed that the obtained synthetic microstructures qualitatively resemble those observed by Hayakawa et al. [16] on doped zirconia at room temperature, showing twin lengths and twin widths in the micrometer and 50 nm ranges, respectively. Note however that a more quantitative comparison to experimental observations is not possible as (i) we don't know any in situ microstructure observations at the high temperatures considered for pure tetragonal zirconia, and (ii) dopants modify the tetragonality ratio and the internal stress level (see Section 5.2) and thus possibly also the microstructure resulting from the SPT.

In Fig. 4, some characteristics of $D^{sm}$ are presented. In Fig. 4(a) the mean number of twins in each realization of the stochastic microstructure model is shown. As might be expected, the number of crystal variants – whose type depends on the local stacking sequence – that are located in both types of band are approximately twice the number of variants that are only in one type of band. It is noteworthy that the aspect ratio $\alpha_{AR}$ has a stronger influence on the number of twins than the number (and therefore the size) of Voronoi cells $n^{Vo}$. This fact is further illustrated when considering the mean volume-equivalent diameters in Fig. 4(b).
<figure><img src="image_7.png" /><figcaption>Fig. 4. (a) Plot of the mean number of twins in each realization, and (b) box plot of the volume-equivalent diameter of the twins, for each variant type and model configuration in the dataset D<sup>sm</sup>.</figcaption></figure>

# 4. Full-field computation

## 4.1. Computational method

Significantly high internal stress is expected in *t* zirconia. On the one hand, this is due to the *c* → *t* phase transformation occurring at ~2300°C, which leads to changes in the crystal lattice parameters and resulting in stress-free deformations with a marked anisotropy. The second source of internal stresses is the anisotropic thermal dilation of the *t* phase, leading to anisotropic strain during subsequent cooling. These two sources of stress-free and incompatible deformation add up. In this work, only transformation, thermal and elastic strain are considered, no plasticity nor damage is introduced. The constitutive law of the material is thus similar to a thermo-elastic one, and reads

$$ \sigma(\mathbf{x}) = \mathbf{C}(\mathbf{x}) : \varepsilon^e(\mathbf{x}) = \mathbf{C}(\mathbf{x}) : (\varepsilon(\mathbf{x}) - \varepsilon^0(\mathbf{x})) \quad (8) $$

with

$$ \varepsilon^0(\mathbf{x}) = \varepsilon^r(\mathbf{x}) + \varepsilon^{th}(\mathbf{x}) \quad (9) $$

and *x* denoting some position inside the herringbone microstructure and $\varepsilon^e$, $\varepsilon^r$, $\varepsilon^{th}$ and $\varepsilon$ the elastic, transformation, thermal and total strain.

To solve for the stress fields, we use the full-field method based on the spectral method initially proposed by Moulinec and Suquet [29] and the associated open-source numerical code CRaFT [5]. In this method, the computational domain $W$ is discretized into 256 × 256 × 256 voxels (or Fourier points). This discretization determines a regular grid in the cartesian space and a corresponding grid in the Fourier space. The method also requires the selection of a linear reference medium of stiffness $C^0$. The heterogeneous problem of a polycrystal exhibiting a different elastic stiffness $C$ at each position $x$, as in Eq. (8), is rewritten equivalently as a homogeneous problem with the arbitrary homogeneous stiffness $C^0$ and an additional stress-free strain (or polarization) field $\varepsilon^0(x)$, which is unknown. The solution is given by a convolution of the Green tensor associated to $C^0$ with the polarization field of interest. In the Fourier space, this convolution turns into a direct product. The Fourier transform of the Green operator associated with the reference medium can be readily calculated. Due to highly optimized numerical implementations of FFTs, the numerical efficiency of this method is very high. An iterative scheme must be implemented to obtain, upon convergence, the compatible strain field associate with a kinematically admissible displacement field, that minimizes the average of the local strain energies, under the constraint imposed by stress field balance [48].

## 4.2. Material behavior and loading conditions

The spectral method described above is limited to computational domains submitted to periodic boundary conditions as discussed e.g. in

Lebensohn and Rollett [24] and Lucarini et al. [25]. It is only very recently that a method allowing non-periodic boundary conditions has been proposed [35]. Therefore, periodic microstructures are often used in combination with this FFT method, e.g. see Lebensohn et al. [23]. This prevents from having an impact of the applied boundary conditions on the mechanical fields, but at the same time this also limits the randomness of the studied microstructures. Periodicity of synthetic herringbone microstructures cannot be achieved. When computing the mechanical response of polycrystalline aggregate, as a rule of thumb, specific boundary conditions are often considered to have an impact on the outer shell of the computation domain over a thickness of typically one twin domain. Therefore, a layer of 20-voxels thickness at the outer surface of the computational domain has been removed for the statistical analysis presented below. Moreover, we have verified that the results are similar if a thicker outer shell is removed.

The results presented below are obtained for cooling. Zirconia crystallizes into a polycrystalline aggregate made of cubic crystals at about 2700 °C. As cubic crystals exhibit an isotropic thermal dilation, residual stresses due to cooling are not expected in the cubic domain down to 2300 °C, the temperature at which the $c \to t$ transformation is engaged. Here, we consider a simple scenario in which the herringbone $t$ microstructure is formed at once at 2300 °C in a parent $c$ stress-free single crystal. The $t$ zirconia is then further cooled down by $\Delta T = -1000$ °C, i.e. down to 1300 °C which lies still higher than the $t \to m$ transformation temperature for pure zirconia. All results are given at this target temperature, which corresponds more or less to the operating conditions of large zirconia-rich refractory blocks used at an industrial scale in the glass furnaces [34].

The residual stress field at 1300 °C depends on the Bain strain during the $c \to t$ transformation, the thermal dilation of the $t$ phase, the elastic behavior of $t$ crystals, and of course also on the specimen microstructure.

The herringbone microstructure for $t$ zirconia has been mostly observed for doped zirconia, as dopants such as yttrium stabilize the $t$ phase at lower temperatures. Lattice parameters of $c$ and $t$ phases depend on the dopant concentration. Here, we use the data from Mamiyand et al. [27]:

$$ a_c = 5.27 \AA, \quad a_t = b_t = 5.141 \AA, \quad c_t = 5.2609 \AA. \quad (10) $$

The transformation strain for a $z$-variant is therefore

$$ \varepsilon^{tr} = \begin{bmatrix} -2.448 \times 10^{-2} & 0 & 0 \\ 0 & -2.448 \times 10^{-2} & 0 \\ 0. & 0. & -1.727 \times 10^{-3} \end{bmatrix}. \quad (11) $$

As for the thermal dilation of the $t$ phase, several published values have been listed in Table 2 of Guinebretière et al. [14] for bulk specimens and powders. Here, the data of zirconia powders measured by Lang [22] between 1150 °C and 1700 °C was considered with three different compositions. It was found that the lattice parameters evolve almost linearly with temperature:
<caption>**Table 1**</caption>

Mean value of the normal stress $\sigma_c$, expressed in MPa, for various microstructures. Full-field results for the aspect ratios $\alpha_{AR} = 10$ and $\alpha_{AR} = 1$ and a single Voronoi cell. Exact results for rank-2 laminates, corresponding to infinite aspect ratio ($\alpha_{AR} = \infty$) are also provided. The color code indicates the microstructures with similar thermo-mechanical response. For example, Fig. 5a provides the mechanical fields for the red microstructures and Fig. 5b those for the green microstructures. Other microstructures look geometrically very similar but the stress distribution might differ significantly.

<table><thead><tr><th></th><th>y<sup>(110)</sup> x<sup>(101)</sup> x<sup>(101)</sup> z</th><th>x<sup>(110)</sup> y<sup>(101)</sup> y<sup>(101)</sup> z</th><th>x<sup>(110)</sup> z<sup>(110)</sup> z<sup>(101)</sup> y</th></tr></thead><tbody><tr><td>α<sub>AR</sub> = ∞</td><td>0</td><td>0</td><td>0</td></tr><tr><td>α<sub>AR</sub> = 10</td><td>-200</td><td>-32</td><td>-200</td></tr><tr><td>α<sub>AR</sub> = 1</td><td>-540</td><td>-700</td><td>-700</td></tr><tr><td></td><td>y<sup>(110)</sup> x<sup>(101)</sup> x<sup>(101)</sup> z</td><td>x<sup>(110)</sup> y<sup>(101)</sup> y<sup>(101)</sup> z</td><td>x<sup>(110)</sup> z<sup>(110)</sup> z<sup>(101)</sup> y</td></tr><tr><td>α<sub>AR</sub> = ∞</td><td>-1939.1</td><td>-130.1</td><td>-3841.7</td></tr><tr><td>α<sub>AR</sub> = 10</td><td>-1992</td><td>-130</td><td>-3514</td></tr><tr><td>α<sub>AR</sub> = 1</td><td>-2625</td><td>-1372</td><td>-1798</td></tr><tr><td></td><td>y<sup>(110)</sup> x<sup>(101)</sup> x<sup>(101)</sup> z</td><td>x<sup>(110)</sup> y<sup>(101)</sup> y<sup>(101)</sup> z</td><td>x<sup>(110)</sup> z<sup>(110)</sup> z<sup>(101)</sup> y</td></tr><tr><td>α<sub>AR</sub> = ∞</td><td>-1883.4</td><td>-402.8</td><td>-402.8</td></tr><tr><td>α<sub>AR</sub> = 10</td><td>-1850</td><td>-376</td><td>-1850</td></tr><tr><td>α<sub>AR</sub> = 1</td><td>-881</td><td>-735</td><td>-735</td></tr><tr><td></td><td>y<sup>(110)</sup> x<sup>(101)</sup> x<sup>(101)</sup> z</td><td>x<sup>(110)</sup> y<sup>(101)</sup> y<sup>(101)</sup> z</td><td>x<sup>(110)</sup> z<sup>(110)</sup> z<sup>(101)</sup> y</td></tr><tr><td>α<sub>AR</sub> = ∞</td><td>-2137.4</td><td>-618.0</td><td>-3415.4</td></tr><tr><td>α<sub>AR</sub> = 10</td><td>-2047</td><td>-755</td><td>-3166</td></tr><tr><td>α<sub>AR</sub> = 1</td><td>-1083</td><td>-1560</td><td>-1656</td></tr><tr><td></td><td>y<sup>(110)</sup> x<sup>(101)</sup> x<sup>(101)</sup> z</td><td>x<sup>(110)</sup> y<sup>(101)</sup> y<sup>(101)</sup> z</td><td>x<sup>(110)</sup> z<sup>(110)</sup> z<sup>(101)</sup> y</td></tr><tr><td>α<sub>AR</sub> = ∞</td><td>-2129.9</td><td>-3415.4</td><td>-618.0</td></tr><tr><td>α<sub>AR</sub> = 10</td><td>-2156</td><td>-3166</td><td>-2047</td></tr><tr><td>α<sub>AR</sub> = 1</td><td>-2630</td><td>-1656</td><td>-1083</td></tr><tr><td></td><td>y<sup>(110)</sup> x<sup>(101)</sup> x<sup>(101)</sup> z</td><td>x<sup>(110)</sup> y<sup>(101)</sup> y<sup>(101)</sup> z</td><td>x<sup>(110)</sup> z<sup>(110)</sup> z<sup>(101)</sup> y</td></tr><tr><td>α<sub>AR</sub> = ∞</td><td>-410.5</td><td>-3841.7</td><td>-130.1</td></tr><tr><td>α<sub>AR</sub> = 10</td><td>-529</td><td>-3514</td><td>-130</td></tr><tr><td>α<sub>AR</sub> = 1</td><td>-749</td><td>-1798</td><td>-1372</td></tr><tr><td></td><td>y<sup>(110)</sup> x<sup>(101)</sup> x<sup>(101)</sup> z</td><td>x<sup>(110)</sup> y<sup>(101)</sup> y<sup>(101)</sup> z</td><td>x<sup>(110)</sup> z<sup>(110)</sup> z<sup>(101)</sup> y</td></tr><tr><td>α<sub>AR</sub> = ∞</td><td>-1883.4</td><td>-402.8</td><td>-402.8</td></tr><tr><td>α<sub>AR</sub> = 10</td><td>-1850</td><td>-376</td><td>-1850</td></tr><tr><td>α<sub>AR</sub> = 1</td><td>-881</td><td>-735</td><td>-735</td></tr></tbody></table>

$$ \frac{da_t}{dT} = 6.36 \times 10^{-5} \AA K^{-1}, \quad \frac{dc_t}{dT} = 7.57 \times 10^{-5} \AA K^{-1}. \qquad (12) $$

From these values, one can already anticipate that the thermal strain, due to dilation, will have a limited effect on the residual stress field. Indeed, with $\Delta T = -1000^{\circ}C$, the difference in the longitudinal thermal strain along $a_t$ and $c_t$ directions is approximately one order of magnitude smaller than the difference for the transformation strain.

Concerning the elastic stiffness of $t$ crystals, several values have been proposed in the literature, see for example Table 2 of Fadda et al. [8]. The values of the components can vary significantly, for example $C_{12}$ and $C_{13}$ can differ by a factor ~4, and the largest component $C_{11}$ by 70%, depending on the author. Here, we use the values measured by Kisi and Howard [19] and also used as reference by Zhao et al. [51] and Mamivand et al. [27]:

$$ \begin{aligned} C_{1111} &= 327 \text{ GPa}, & C_{1122} &= 100 \text{ GPa}, & C_{1133} &= 62 \text{ GPa}, \\ C_{3333} &= 264 \text{ GPa}, & C_{1212} &= 59 \text{ GPa}, & C_{2323} &= 64 \text{ GPa}. \end{aligned} \qquad (13) $$

# 5. Results

## 5.1. Results for a single Voronoi cell

In this section, we start with the microstructures from the simpler dataset $D^{hs}$, i.e. those for which only a single Voronoi cell was used (i.e. $n^{Vo} = 1$) as in Figs. 3a and 3b. Each cell is defined by a number of geometrical/crystallographical features, namely a stacking sequence, the twin planes inside the bands and the band boundary plane. Three different stacking sequences are possible, named as $y-x \leftrightarrow x-z$, $x-y \leftrightarrow y-z$ and $x-z \leftrightarrow z-y$ where the repeated variant name corresponds to the variant shared by the two bands. To indicate the twin planes and the band boundary orientation, we use a notation such as $y^{(011)} x^{(110)} x^{(101)} z$ indicating a (011) twin planes between $y$ and $x$ variants in the first band, a (101) twin plane between variants $x$ and $z$ in the second band, and a plane (110) as band boundary. In total, 18 different combinations are possible, 3 possible stacking sequences and

for each stacking sequence 6 possible permutations of plane indices, all listed in Table 1. Evidently, as the initial $c$ crystal exhibits cubic symmetries, and since we are considering only isotropic loadings (thermal loading) of the herringbone $t$ microstructures, some of these 18 microstructures generate similar stress and strain fields that are only rotated by 90° with respect to one of the cubic axes. This is for example the case for $y^{(110)} x^{(011)} x^{(101)} z$ and $x^{(110)} y^{(101)} y^{(011)} z$. Accounting for these symmetries, we end up with 6 different microstructures and 6 different mechanical responses. Symmetrically equivalent microstructures are grouped by a color code in Table 1.

In the following, we discuss the results from the full-field computations, performed for two aspect ratios, $\alpha_{AR} = 10$ and $\alpha_{AR} = 1$. We expressed all stress tensors in a reference frame attached to the $t$ crystal lattice. Here, we will mostly concentrate on the normal stress components denoted $\sigma_a$, $\sigma_b$, $\sigma_c$. For example, $\sigma_c$ is the normal stress along the $c$-axis of a $t$-crystal, i.e. acting on a (001)$_t$ lattice plane along the direction [001]$_t$. For each generated microstructure, we have computed the stress distributions in twin domains with a given lattice orientation (i.e. same variant) and a given band, and estimated the mean and dispersion values of these distributions. The results for the mean normal stress $\sigma_c$ within each variant are reported in Table 1. Exact results for rank-2 laminates, corresponding to infinite aspect ratio ($\alpha_{AR} = \infty$) are also provided, see Appendix A.

(i) It can be remarked that the obtained mean stress values for $\sigma_c$ are huge. In all microstructures except for the ones indicated in the first line of the table (yellow cells), there is at least one variant for which $\sigma_c$ lies in the GPa range, and up to -3.8 GPa. In general, only non-positive values (compressive stress) are obtained.

(ii) All mean normal stress components are reported in Table 2 for the case $\alpha_{AR} = 10$ and for 6 representative microstructures (the 12 remaining ones can be obtained by symmetry operations as discussed above). Similar to $\sigma_c$, the values for $\sigma_a$ and $\sigma_b$ also lie
<table><caption>**Table 2**
Mean value of the normal stresses components $\sigma_a$, $\sigma_b$ and $\sigma_c$, expressed in MPa, for each set of equivalent microstructures. Full-field results for an aspect ratio $\alpha_{AR} = 10$ and a single Voronoi cell. The color code is similar than for Table 1.</caption><thead><tr><th></th><th>y <sup>(110)</sup> / x</th><th><sup>(011)</sup> &harr; x</th><th><sup>(110)</sup> / z</th></tr></thead><tbody><tr><td>σ<sub>a</sub></td><td>31</td><td>124</td><td>97</td><td>144</td></tr><tr><td>σ<sub>b</sub></td><td>144</td><td>97</td><td>124</td><td>31</td></tr><tr><td>σ<sub>c</sub></td><td>-200</td><td>-32</td><td>-32</td><td>-200</td></tr></tbody><thead><tr><th></th><th>y <sup>(110)</sup> / x</th><th><sup>(101)</sup> &harr; x</th><th><sup>(011)</sup> / z</th></tr></thead><tbody><tr><td>σ<sub>a</sub></td><td>307</td><td>-1838</td><td>2218</td><td>1593</td></tr><tr><td>σ<sub>b</sub></td><td>254</td><td>194</td><td>271</td><td>1803</td></tr><tr><td>σ<sub>c</sub></td><td>-1992</td><td>-130</td><td>-3514</td><td>-529</td></tr></tbody><thead><tr><th></th><th>y <sup>(101)</sup> / x</th><th><sup>(011)</sup> &harr; x</th><th><sup>(110)</sup> / z</th></tr></thead><tbody><tr><td>σ<sub>a</sub></td><td>257</td><td>3275</td><td>-142</td><td>188</td></tr><tr><td>σ<sub>b</sub></td><td>188</td><td>-142</td><td>3275</td><td>257</td></tr><tr><td>σ<sub>c</sub></td><td>-1850</td><td>-376</td><td>-376</td><td>-1850</td></tr></tbody><thead><tr><th></th><th>y <sup>(101)</sup> / x</th><th><sup>(110)</sup> &harr; x</th><th><sup>(011)</sup> / z</th></tr></thead><tbody><tr><td>σ<sub>a</sub></td><td>124</td><td>1450</td><td>-1445</td><td>87</td></tr><tr><td>σ<sub>b</sub></td><td>2144</td><td>3203</td><td>667</td><td>1844</td></tr><tr><td>σ<sub>c</sub></td><td>-2047</td><td>-755</td><td>-3166</td><td>-2156</td></tr></tbody><thead><tr><th></th><th>y <sup>(011)</sup> / x</th><th><sup>(101)</sup> &harr; x</th><th><sup>(110)</sup> / z</th></tr></thead><tbody><tr><td>σ<sub>a</sub></td><td>1844</td><td>667</td><td>3203</td><td>2144</td></tr><tr><td>σ<sub>b</sub></td><td>87</td><td>-1445</td><td>1450</td><td>124</td></tr><tr><td>σ<sub>c</sub></td><td>-2156</td><td>-3166</td><td>-755</td><td>-2047</td></tr></tbody><thead><tr><th></th><th>y <sup>(011)</sup> / x</th><th><sup>(110)</sup> &harr; x</th><th><sup>(011)</sup> / z</th></tr></thead><tbody><tr><td>σ<sub>a</sub></td><td>1803</td><td>271</td><td>194</td><td>254</td></tr><tr><td>σ<sub>b</sub></td><td>1593</td><td>2218</td><td>-1838</td><td>307</td></tr><tr><td>σ<sub>c</sub></td><td>-529</td><td>-3514</td><td>-130</td><td>-1992</td></tr></tbody></table>

in the GPa range, except for the yellow cells. This time, however, positive (tensile) stress values, up to +3.2 GPa, are also obtained.

(iii) For the sake of conciseness, shear stress values are not given explicitly here. They have been found to be significantly smaller (in absolute value) than normal stresses, generally in the range of tens of MPa. The largest obtained value for $\alpha_{AR} = 10$ is 450 MPa.

(iv) According to Table 1, the aspect ratio of the *t* lamellae plays an important role in the stress distribution. Results for $\alpha_{AR} = 1$ and $\alpha_{AR} = 10$ show some similarities but the obtained values for $\sigma_c$ can be quite different.

(v) Note that in Table 1 some microstructures exhibit similar distribution of $\sigma_c$ beyond the basic symmetry operations discussed above. This is, for example, the case for the gray and green cells, as well as for the orange and blue ones. For these microstructures, however, $\sigma_a$ and $\sigma_b$ do not match together (see Table 2) for a given variant, but we find the same stress values in other variants. The values are simply permuted. For the yellow and red cells, we find the same values for $\sigma_a$ and $\sigma_b$ belonging to the variant not shared by both bands. These values are simply swapped. We also find the same swapped sets of values for gray and green cells.

When looking into the details of the stress distributions within the various microstructures, it is apparent that there are basically three different cases:

1. For some microstructures, a bimodal stress distribution is obtained within the whole microstructure. This is the case for the yellow and red cells in Table 1, where the twin domains belonging to the variant shared by the two bands on one side, and belonging to the two other variants on the other side, exhibit the same mean stress value $\sigma_c$.

2. For other microstructures (blue and orange cells in Table 1), the mean stress in twin domains belonging to the variant shared by both

bands is significantly different, whereas in the other two variants, values of $\sigma_c$ are almost identical.

3. For the remaining configurations (gray and green), four significantly different mean values of $\sigma_c$ are obtained for each variant in each band. In the case of $\alpha_{AR} = 10$, for example in $y^{(110)} / x^{(101)} \leftrightarrow x^{(011)} / z$, the mean stress in twin domains belonging to variant *x* can be relatively small inside the band containing also the variant *y* (namely -0.13 GPa), but huge (in absolute value) inside the other band containing also the variant *z* (namely -3.5 GPa).

Note that the classification above depends on the stress component considered (Table 2). For example in $y^{(101)} / x^{(011)} \leftrightarrow x^{(110)} / z$, the same value for $\sigma_c$ is found within the variant *x* for both bands, but significantly different values for $\sigma_a$ and $\sigma_b$ in the same variant between both bands.

Fig. 5 illustrates the observations above. It shows the 3D field of $\sigma_c$ for a microstructure $y^{(101)} / x^{(011)} \leftrightarrow x^{(110)} / z$ corresponding to Case 1 (bimodal distribution of $\sigma_c$), and a microstructure $y^{(011)} / x^{(110)} \leftrightarrow x^{(101)} / z$ for which 4 different stress values arise (Case 3), both with $\alpha_{AR} = 10$. The probability densities of the residual stress are also provided. Here, the stress heterogeneity inside each twin domain is small, leading to well distinct and narrow peaks in the stress distribution. In $y^{(101)} / x^{(011)} \leftrightarrow x^{(110)} / z$, the *x* variants of both bands show the same stress level. This is also the case for variants *y* and *z*. Contrarily, for $y^{(011)} / x^{(110)} \leftrightarrow x^{(101)} / z$, the 4 variants exhibit different stress states, and in particular the stress in the *x* variant in both bands is significantly different (-3.5 GPa and -0.13 GPa). These results show that in such herringbone microstructures, the mean stress in twin domains with similar crystallographic orientation is not guided essentially by their lattice orientation as often encountered in polycrystalline aggregates e.g. see Gu et al. [13] and Purushottam Raj Purohit et al. [41], but crystal stresses are highly dependent on the local crystal arrangement.

To go one step further in the analysis of the statistical distribution inside the constituting crystals, in Fig. 6 the probability density of $\sigma_a$, $\sigma_b$ and $\sigma_c$ is shown for an aspect ratio $\alpha_{AR} = 10$ when considering all possible microstructures indicated in Table 1. The main findings are:

(i) The distributions of $\sigma_a$ and $\sigma_b$ turn out to be similar, due to the tetragonal symmetry of the crystal lattice.

(ii) The distributions of $\sigma_a$ and $\sigma_b$ show 4 peaks, whereas that of $\sigma_c$ has only 3.

(iii) These three distributions are extremely broad with widths between -4 GPa and 4 GPa.

These results highlight again the huge impact of the local crystallographic and morphologic environment on the stress state of a given crystal in the whole microstructure.

## 5.2. Influence of the aspect ratio of twin domains

We have also investigated the impact of the twin domain aspect ratio on the resulting distribution of residual stress in the various microstructures. Examples of microstructures with aspect ratios $\alpha_{AR} = 1$ and $\alpha_{AR} = 10$ are shown in Fig. 3a and 3b, respectively.

The probability density of $\sigma_c$ for the sequence $y^{(011)} / x^{(101)} \leftrightarrow x^{(110)} / z$ is shown in Fig. 7, for both aspect ratio. Increasing the aspect ratio from 1 to 10 very clearly leads to a narrowing of the stress distribution in all variants. Apparently, for $\alpha_{AR} = 1$, the distribution of $\sigma_c$ in the *y* and *z* variants is well separated (by ~2 GPa), but for $\alpha_{AR} = 10$ they are almost overlapping. Conversely, the *x* variants in both bands behave very similarly for $\alpha_{AR} = 1$, but are clearly distinct for $\alpha_{AR} = 10$.

The case $\alpha_{AR} = 1$ was investigated mostly to track the effects of drastic changes in microstructures. It probably departs from real microstructures due to exhibiting a very extreme morphology. For example,
<figure><img src="image_8.png" /><figcaption>**Fig. 5.** Distribution of $\sigma_c$ in microstructures $y \xrightarrow{\sigma_c^{(101)}} x \xrightarrow{\sigma_c^{(101)}} x \xrightarrow{\sigma_c^{(110)}} z$ (left) and $y \xrightarrow{\sigma_c^{(011)}} x \xrightarrow{\sigma_c^{(110)}} x \xrightarrow{\sigma_c^{(101)}} z$ (right) with one Voronoi cell ($n^{Vo} = 1$) and aspect ratio $\alpha_{AR} = 10$. Spatial 3D stress distributions (top) and corresponding probability distribution functions (bottom) showing 2 or 4 distinct peaks (the arrows indicate the variant for the distribution shown).</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Fig. 6.** Probability density functions averaged over the 18 possible microstructures with $\alpha_{AR} = 10$ and $n^{Vo} = 1$ for (a) $\sigma_a$, (b) $\sigma_b$, and (c) $\sigma_c$.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Fig. 7.** Probability density functions of $\sigma_c$ for the microstructure $y \xrightarrow{\sigma_c^{(011)}} x \xrightarrow{\sigma_c^{(101)}} x \xrightarrow{\sigma_c^{(110)}} z$ with $n^{Vo} = 1$ for the aspect ratios (a) $\alpha_{AR} = 1$ and (b) $\alpha_{AR} = 10$. The arrows indicate which variant the distribution corresponds to.</figcaption></figure>

ple, the experimental observations (electron micrograph) of Hayakawa et al. [16] show aspect ratio more in the range $\alpha_{AR} \sim 10 - 20$, i.e. similar to microstructures investigated in Section 5.1. On the other hand, this shows that the procedure for microstructure generation, described in Section 3, is robust and capable of constructing highly different microstructure types.

Herringbone microstructures with very large aspect ratio can hardly be studied numerically, as this would require a much larger computational domain for the FFT solver. However, microstructures with infinite aspect ratio ($\alpha_{AR} \to \infty$) are easy to handle through the analytical solution for thermo-elastic laminates. The microstructure of such a periodic composite is illustrated in Fig. 8. A rank-1 laminate is composed of a
<figure><img src="image_11.png" /><figcaption>**Fig. 8.** (a) Schematic representation of a rank-1 laminate with normal n, made of two materials with elastic stiffness and stress-free strain C and $\varepsilon^0$, respectively. (b) In the rank-2 laminate with normal n, each layer is itself a rank-1 laminate with effective stiffness and effective stress-free strain $\tilde{C}$ and $\tilde{\varepsilon}^0$.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>**Fig. 9.** Effect of the tetragonality $c_i/a_t$ on the stress component $\sigma_c$ in one of the x twin domains for the microstructure $y^{(110)} x^{(101)} x^{(011)} z$ (gray curves), accounting for SPT only (continuous line) and for SPT plus cooling (dashed line). For $y^{(110)} x^{(011)} x^{(101)} z$ (yellow curve), the value of $\sigma_c$ in all twin domains remains zero regardless of tetragonality, dilation coefficients of the *t* phase and cooling range. Results from the rank-2 analytical solution corresponding to $\alpha_{AR} \to \infty$ and $n^{VO} = 1$.</figcaption></figure>

succession of planar, parallel, infinite and perfectly bonded layers. Each layer has a different elastic stiffness and stress-free strain. The exact effective stiffness and effective stress-free strain of such a composite can be easily calculated. A rank-2 laminate can be constructed similarly by piling-up layers made of rank-1 laminates. Continuing the same process, it is possible to construct rank-*n* laminates, up to infinite order as in Francfort and Murat [9] and Idiart [17]. The elastic response of rank-1 laminates has been investigated first by Backus [3]. The thermo-elastic response has received much less attention. In Appendix A, we extend the analytical solution given in Milton [28] by providing a compact expression of the residual stress field. A well-known exact result concerning the thermo-elastic response of laminates is that the stress and strain fields inside each layer are *uniform*, and therefore the stress distributions shown above for finite $\alpha_{AR}$ values (Figs. 5, 6 and 7) become Dirac shaped for $\alpha_{AR} \to \infty$.

The interest for such a rank-2 laminate structure within the present study is twofold. First, it resembles the herringbone structure described in Section 2, with the finest lamellae representing the twinning process, and the largest (rank-2) lamination the large bands observed in herringbone structures. And second, if at each stacking step, one piles up a very large number of thin layers so that the physical scale of the whole laminate is much larger than the one of each layer, then the analytical solution provided in Appendix A is exact. This condition implies that the aspect ratio of the layers, at each stacking step, is infinite.

We have thus constructed rank-2 laminates for herringbone microstructures following the same crystallographic rules as described in Section 3. The values obtained for $\sigma_c$ are indicated in Table 1.

(i) One can observe that the difference in $\sigma_c$ between the cases $\alpha_{AR} = 10$ (mean values of $\sigma_c$) and $\alpha_{AR} \to \infty$ (exact values of $\sigma_c$) is rather small, generally of the order of a few tens of MPa, the largest difference being 0.33 GPa.

(ii) Recalling that the analytical solution for a rank-*n* laminate predicts a uniform stress in each twin variant when $\alpha_{AR} \to \infty$, the distribution shown in Figs. 5, 6, 7 would be Dirac shaped for infinite aspect ratios. One may understand the effect of the aspect ratio on the stress distribution as a deviation from this analytical solution, with strong stress heterogeneities arising at the end of each twin for shorter twins leading to broader stress distributions.

(iii) The 3 microstructures highlighted in yellow in Table 1 turn out to be very specific. Already in the cases $\alpha_{AR} = 1$ and $\alpha_{AR} = 10$, much lower $\sigma_c$ values were found compared to other microstructures. In the case $\alpha_{AR} \to \infty$, one obtains that $\sigma_c$ is strictly zero, and we have checked that all other stress components also vanish. These specific microstructures thus correspond to a case for which both SPT and anisotropic thermal dilation of the *t* phase do not lead to any residual stress, thus no stored elastic energy, and also no reason for the appearance of cracks. It can be checked that this corresponds to configuration for which both transformation strain and thermal strain are compatible, i.e. transformation strain and thermal strain alone would not lead to any displacement discontinuity between the twins, and therefore there is no need for additional elastic accommodation during the cooling process. Also, interestingly, these 3 microstructures correspond to the experimental ones described in Section 2, based on the possible twinning planes associated with the crystal structures.

(iv) We have verified with the analytical solution for the laminate that the sign of the Miller indices used to describe the twin planes and band boundaries has no effect on the stress field, i.e. calculating the response of a microstructure comprising e.g. a (110) plane or a (\#110) plane gives exactly the same results. The same behavior has also been found with the full-field numerical approach, even if the microstructures are geometrically different. This is why we have removed the sign of all Miller indices in Tables 1 and 2.

(v) Using the analytical solution for laminates, one can also easily estimate the relative importance of the transformation strain associated with the SPT and thermal dilation in the build-up of the residual stress field. Performing calculations accounting for the SPT either with cooling ($\Delta T = -1000°$ K) or without cooling ($\Delta T = 0°$ K), both resulting residual stress fields are very close. All stress components in the case that included cooling are slightly smaller by about 5% (in absolute value) than without cooling. The residual stress field is thus essentially governed by the transformation strain. As already introduced after equation (12), this comes from the fact that the deviatoric strain due to thermal dilation is much smaller than the deviatoric strain due to SPT, and therefore the internal stress level is mostly driven by the transformation strain.

(vi) Among all material parameters listed in Section 4.2, all results presented above do not depend directly on the lattice parameters, nor on the dilation coefficients. Only the tetragonality $c_i/a_t$ and the ratio of the dilation coefficients along $c_i$ and $a_t$ of the *t* phase come in play. The lattice parameter of the *c* phase does not have any influence on the stress field for the calculation conditions adopted here. The tetragonality is a particularly important parameter as it is routinely adjusted by adding dopants such as yttria in order to stabilize the *t* phase down to room temperature, see e.g. Krogstad
<figure><img src="image_13.png" /><figcaption>**Fig. 10.** Box plots describing the distributions of $\sigma_c$ in the Voronoi cells comprising the sequence $y = x^{(101)} \ x^{(011)} \ x^{(110)}$ for nine generated microstructures with $\alpha_{AR} = 10$ and $n^{Vo} = 20$ Voronoi cells. The colors correspond to the variants (i.e. x, y and z) as indicated by the arrows.</figcaption></figure>

et al. [20]. Fig. 9 shows the effect of the tetragonality on the value of $\sigma_c$ in one illustrative microstructure $y = x^{(110)} \ x^{(011)} \ x^{(011)} \ x^{(110)} \ z$ when the zirconia has been submitted to SPT only (thus forming the herringbone structure) or when it has been submitted to SPT and cooling ($\Delta T = -1000$ °C) as above. This figure illustrates the small effect of thermal cooling compared to the SPT, i.e. a constant difference of 216 MPa being found between both curves. For both cases, the residual stress field evolves linearly with $c_i/a_i$, and a similar linear response (not shown) is obtained when modifying the ratio of the dilation coefficients of the t phase. Again, the specific $y = x^{(110)} \ x^{(011)} \ x^{(101)} \ z$ microstructure (yellow in Tables 1 and 2) shows a very specific behavior with vanishing residual stress (all components in all twin domains), regardless of tetragonality, cooling range, and ratio of dilation coefficients for the t phase.

### 5.3. Microstructure comprising several Voronoi cells

Finally, we have investigated the effect of having several Voronoi cells in the computational domain, each of which contains a herringbone structure of a random configuration. These generated microstructures match the experimental ones more closely. The different Voronoi cells may reflect a process in which the phase transition starts simultaneously at different positions within a large parent c crystal at high temperature. Several herringbone microstructures thus grow concurrently and meet at some point where a new boundary is formed between the individual herringbone structures. These boundaries are modeled by the facets of the Voronoi cells described in Section 3.

We focus in the following on the stacking sequence $y = x^{(101)} \ x^{(011)} \ x^{(110)} \ z$ with an aspect ratio $\alpha_{AR} = 10$. In Fig. 10, nine generated microstructures are considered each comprising $n^{Vo} = 20$ Voronoi cells. For these, box plots of the distribution of $\sigma_c$ within the Voronoi cells containing the sequence $y = x^{(101)} \ x^{(011)} \ x^{(110)} \ z$ are shown. It is apparent that quite strong fluctuations occur in all the 3 variants from one microstructure to another, not only for the mean values of $\sigma_c$, but also for the dispersion of the distributions. For example, the standard deviation of $\sigma_c$ can vary from 100 MPa to 600 MPa in variant y, depending on the microstructure.

Fig. 11 shows the distributions of $\sigma_c$, for the same sequence, still in microstructures comprising 20 Voronoi cells, but now as a statistical average over many replications from the stochastic microstructure model. These distributions can be compared to that given in Fig. 5a which was calculated for a microstructure comprising a single Voronoi cell for the same sequence. Clearly, the presence of boundaries randomly delimit-

<figure><img src="image_14.png" /><figcaption>**Fig. 11.** Probability density functions of $\sigma_c$ in one Voronoi cell comprising the sequence $y = x^{(101)} \ x^{(011)} \ x^{(110)} \ z$ averaged over many generated microstructures with $\alpha_{AR} = 10$ and $n^{Vo} = 20$ Voronoi cells.</figcaption></figure>

ing two specific sequences of variants significantly increases the spread of the distributions of residual stress.

A general conclusion of this section is that when different herringbone structures grow simultaneously in the same initial c crystal, the residual stress field is globally of higher intensity, and stress distributions are broader compared to the case with a single Voronoi cell (i.e. a single nucleation point of the SPT). This is due to the additional strain incompatibilities between all herringbone structure, compared to the case where only one Voronoi cell is present.

However, there is still a very specific case with the microstructures $y = x^{(110)} \ x^{(011)} \ x^{(101)} \ z$ highlighted in yellow in Table 1. For $\alpha_{AR} \to \infty$, besides vanishing residual stresses for the case of a single herringbone structure as found in Section 5.2, the effective stress-free strain $\bar{\epsilon}^0$ of such microstructures turns out to be isotropic (strain value is found to be $-2.36 \times 10^{-2} \text{i}$, with i the second rank isotropic tensor). Therefore, when all Voronoi cells of the microstructure contain such a sequence, even oriented differently, the final residual stress field still vanishes, with no elastic energy stored, regardless of the position and orientation of the boundary between the different herringbone structures (i.e. facets of the Voronoi tessellation) because the effective stress-free strain inside each Voronoi cell is equal and isotropic, thus compatible. This might be the reason why the Laue microdiffraction patterns measured at high temperature in pure tetragonal zirconia in Purushottam Raj Purohit et al. [39] seem to indicate very little lattice distortion.

## 6. Concluding remarks

In this work, we have investigated the residual stress field that may arise during the $c \to t$ phase transition and thermal cooling in zirconia. For doing this, we have proposed a new scheme for the generation of three-dimensional herringbone microstructures that are representative of the tetragonal phase. These synthetic microstructures can contain one or more herringbone structures, each in a separate Voronoi cells, and the individual twin domains inside these structures can exhibit various aspect ratios. Such microstructures have been investigated numerically by a spectral method for two twin domain aspect ratios, $\alpha_{AR} = 1$ and $\alpha_{AR} = 10$. The exact behavior in the case of infinite twin domain aspect ratio has also been provided based on the analytical solution for rank-2 thermo-elastic laminates. The main results can be summarized as follows:

(i) 18 individual herringbone microstructures have been generated and investigated. They can be grouped into 6 different sets each containing 3 similar microstructures, which, due to the isotropic
thermal loading, differ only by symmetry operations of the initial c crystal.

(ii) In all investigated cases – except the specific ones described below – normal residual stresses in each type of twin domain reach huge values, with mean values lying in the range [−3.8 GPa : +3.8 GPa]. Shear stresses are found to be much smaller by one order or magnitude. The mean residual stress along the c axis, $\sigma_c$, is always negative, while the mean values of $\sigma_a$ and $\sigma_b$ are mostly positive, although some negative values are also found for some specific microstructures.

(iii) The stress distribution over the whole microstructure can have 2, 3 or 4 well-defined maxima, depending on the microstructure type. In particular, twin domains of the same variant can exhibit significantly different residual stress levels depending on their local environment, i.e. to which band they belong. A consequence of this result is that, if local stress measurements can be performed once, e.g., using high resolution EBSD [38], the interpretation of the results will require thorough knowledge of the local specimen microstructure. Such an experimental comparison would be highly valuable but data acquisition remains highly challenging due to the high temperature range. A possibility would be to use the high temperature Laue microdiffraction setup recently developed at the European Synchrotron Radiation Facility ESRF. The proof of concept for such experiment has been demonstrated on a pure zirconia polycrystal at 1250 °C, i.e. under its tetragonal state. Values of elastic strain of up to 1% have been recorded, which, according to the elastic properties of zirconia, correspond to stress levels of a few GPa [39].

(iv) The effect of thermal expansion is negligible compared to the one of the transformation strain (~5%) even for a significant cooling of $\Delta T = -1000$ °K.

(v) Microstructures with twin domains having a small aspect ratio ($\alpha_{AR} = 1$) globally exhibit similar features than the ones with higher aspect ratio ($\alpha_{AR} = 10$), but show a larger spread of the stress distributions.

(vi) When several herringbone structures live together in a single microstructure, represented in this work by several Voronoi cells, the distributions of the residual stress field remain globally similar to those for a single herringbone structure, but with broader stress distributions. We expect that the stress distributions in a specimen of t zirconia arising from a polycrystal of c zirconia (at high temperature), will be similar to those investigated here and comprising many Voronoi cells.

(vii) Finally, three very specific microstructures have been identified, namely $y(\frac{110}{101}) \ x(\frac{101}{101}) \ x(\frac{101}{101}) \ y(\frac{101}{101}) \ z(\frac{101}{101}) \ z(\frac{101}{101}) \ y(\frac{101}{101})$. In the case of infinite aspect ratio, they all lead to isotropic effective stress-free strain and no residual stress at the twin domain level, regardless of the number of herringbone structures included in the whole microstructure. Therefore, these specific microstructures do not store any elastic energy and should not develop any cracks during the SPT and cooling processes, assuming the microstructure transforms to the t phase all at once. For finite aspect ratio, the obtained residual stress is much smaller than those for all other microstructures. Interestingly, these microstructures are the ones expected in real t zirconia due to the twinning plane being compatible with the crystal structures.

We aim to expand the present study to residual stress fields resulting from the $t \to m$ SPT in pure zirconia. The m zirconia microstructures are slightly more difficult to construct due to a larger number of variants (namely 24) compared to the t case, but this should be very helpful to interpret recent experimental data acquired with electron microscopy [31] and with high resolution x-ray diffraction experiments at high temperatures obtained on synchrotron beamlines [14,40,39].

### CRediT authorship contribution statement

**Lukas Petrich:** Writing – review & editing, Methodology, Visualization, Formal analysis, Writing – original draft, Investigation, Software.
**Katell Derrien:** Writing – review & editing, Methodology, Writing – original draft, Investigation, Visualization, Software, Formal analysis.
**Volker Schmidt:** Supervision, Writing – review & editing, Funding acquisition.
**René Guinebretière:** Funding acquisition, Writing – review & editing, Conceptualization.
**Hervé Moulinec:** Investigation, Software, Formal analysis.
**Olivier Castelnau:** Writing – review & editing, Methodology, Writing – original draft, Funding acquisition, Visualization, Formal analysis, Investigation, Conceptualization.

### Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

### Acknowledgement

This work was carried out within the framework of the HoTMiX research program funded by the agencies ANR (France) and DFG (Germany) under grant No. ANR-09-CE09-0035 and SCHM 997/37-1, respectively, as well as by the NanOX-ML program funded by ANR and FWF (Austria) under grant No. ANR-23-CE09-0041.

### Appendix A. Exact results for a rank-2 thermo-elastic laminate

In this section, we recall the expression for the effective behavior of a rank-n thermo-elastic laminate [28] and we provide compact expressions for the residual stress field. Schematic representation of the microstructure of rank-1 and rank-2 laminates are shown in Fig. 8. A rank-1 laminate is constituted by a random stack of infinite, planar, and perfectly bonded layers exhibiting different thermo-elastic properties. A rank-2 laminate is constituted by random stack of rank-1 laminates. The laminate is submitted to uniform stress or strain boundary condition and uniform temperature change, see Milton [28] for more details.

In the rank-1 laminate, the unit vector normal to the layers is denoted **n**. These layers are denoted ‘mechanical phases’ in the following. Each layer has a different elastic stiffness **C** and stress-free strain $\epsilon^0$. The effective stiffness and effective stress-free strain of such a composite are denoted **C̄** and $\epsilon^0$ respectively.

In a rank-2 laminate, each layer is itself a rank-1 laminate, and **n** now denotes the unit normal to the planar boundary between each of these rank-1 laminates. The stiffness and stress-free strain of each layer of the rank-2 laminate are thus those computed for the constituting rank-1 laminates. Solving for a rank-n laminate can thus be done iteratively.

The mechanical phases constituting the laminate are supposed to exhibit a thermo-elastic behavior

$$\epsilon(\mathbf{x}) = \mathbf{S}^{(r)} : \sigma(\mathbf{x}) + \epsilon^{0(r)}, \quad \sigma(\mathbf{x}) = \mathbf{C}^{(r)} : (\epsilon(\mathbf{x}) - \epsilon^{0(r)}) \quad (A.1)$$

with $\mathbf{C}^{(r)}$ and $\mathbf{S}^{(r)}$ the elastic stiffness and compliance tensors, respectively, and $\epsilon^{0(r)}$ the stress-free (or thermal) strain, for the phase designed by its index $(r)$. Furthermore, $\sigma(\mathbf{x})$ is the stress tensor at position $\mathbf{x}$ within the laminate, while $\epsilon(\mathbf{x})$ is the total strain tensor at $\mathbf{x}$, i.e. $\epsilon(\mathbf{x})$ is the sum of the elastic and stress-free deformations. Here, we use the notation $\chi^{(r)}$ when the field $\chi(\mathbf{x})$ is uniform within phase $r$, i.e. when $\chi(\mathbf{x}) = \chi^{(r)}$ for all locations $\mathbf{x}$ in the phase $(r)$. When the stress-free strain is due to thermal dilation, it holds

$$\epsilon^{0(r)} = \alpha^{(r)} \Delta T \quad (A.2)$$

with $\alpha^{(r)}$ the local dilation modulus of phase $(r)$ and $\Delta T$ the temperature change, supposed homogeneous within the laminate. We consider hereafter the case of general anisotropy for both the elastic stiffness and thermal strain (or dilation tensor). The elastic strain field is given by
$$ \epsilon^e(x) = S^{(r)} : \sigma(x) \qquad (A.3) $$

and is due to both the elastic accommodation of the incompatible thermal strain and to elastic strain due to the applied macroscopic stress. The effective behavior is given by

$$ \bar{\epsilon} = \bar{S} : \bar{\sigma} + \bar{\epsilon}^0 \qquad (A.4) $$

with $\bar{S}$ the effective compliance and $\bar{\epsilon}^0$ the effective thermal strain. One defines for a purely elastic problem (i.e. for which $\epsilon^{(r)} = 0$ for all $r$), e.g. when $\Delta T = 0$) the stress concentration tensor $\mathbf{B}$ relating local and effective stresses

$$ \sigma(x) = \mathbf{B}(x) : \bar{\sigma} \qquad (A.5) $$

so that, for the thermo-elastic problem, the local stress can be expressed as

$$ \sigma(x) = \mathbf{B}(x) : \bar{\sigma} + \sigma^{res}(x) \qquad (A.6) $$

with $\sigma^{res}$ the field of residual stress. One can similarly define a strain localization tensor $\mathbf{A}$ such that, for a purely elastic problem,

$$ \epsilon(x) = \epsilon^e(x) = \mathbf{A}(x) : \bar{\epsilon} \qquad (A.7) $$

It is thus true that, for any microstructure,

$$ \bar{C} = \langle C : A \rangle, \quad \bar{S} = \langle S : B \rangle, \quad \bar{\epsilon}^0 = \langle \epsilon^0 : B \rangle, \qquad (A.8) $$

together with $\langle A \rangle = \langle B \rangle = \mathbf{I}$ where $\langle \cdot \rangle$ denotes the volume average over the whole laminate volume.

Accounting for their specific microstructure, the effective elastic stiffness $\bar{C}$ for rank-1 laminates is given by the following exact expression [28]

$$ [C_0(C_0\mathbf{I} - \bar{C})^{-1} - \mathbf{\Gamma}^1(\mathbf{n})]^{-1} = \langle [C_0(C_0\mathbf{I} - \mathbf{C}^{(r)})^{-1} - \mathbf{\Gamma}^1(\mathbf{n})]^{-1} \rangle \qquad (A.9) $$

with $C_0$ an arbitrary stiffness used to introduce the polarization field and $\mathbf{I}$ the fourth order identity tensor with components $I_{ijkl} = \frac{1}{2}(\delta_{ik}\delta_{jl} + \delta_{il}\delta_{jk})$. The projection tensor $\mathbf{\Gamma}^1$ only depends on the unit normal $\mathbf{n}$

$$ \mathbf{\Gamma}_{ijlm}^{1} = \frac{1}{2}(n_{i}\delta_{jl}n_{m} + n_{i}\delta_{jm}n_{l} + n_{j}\delta_{il}n_{m} + n_{j}\delta_{im}n_{l}) - n_{i}n_{j}n_{l}n_{m} \qquad (A.10) $$

and the complementary tensor $\mathbf{\Gamma}^2$ is defined as

$$ \mathbf{\Gamma}^1 + \mathbf{\Gamma}^2 = \mathbf{I} \qquad (A.11) $$

A well-known remarkable result is that, for homogeneous boundary conditions, the stress and strain in laminates are uniform within each layer, i.e.

$$ \sigma(x) = \sigma^{(r)}, \quad \epsilon(x) = \epsilon^{(r)} \qquad (A.12) $$

for all locations $x$ within the phase $(r)$ with some components being equal to their effective counterparts. This can be expressed as

$$ \Gamma_1 : \sigma(x) = \Gamma_1 : \bar{\sigma}, \quad \text{and} \quad \Gamma_2 : \epsilon(x) = \Gamma_2 : \bar{\epsilon} \qquad (A.13) $$

for all locations $x$. Considering the purely elastic problem, remarking that the local constitutive relation can be written $(\mathbf{\Gamma}^1 + \mathbf{\Gamma}^2) : \sigma = \mathbf{C} : (\mathbf{\Gamma}^1 + \mathbf{\Gamma}^2) : \epsilon$, and using the properties indicated above, it can be shown that the stress concentration tensor (which is uniform within phases) is given by the following expression

$$ \mathbf{B}^{(r)} = [\mathbf{\Gamma}^1 : \mathbf{S}^{(r)} - \mathbf{S}^{(r)} : \mathbf{\Gamma}^2]^{-1} : [\mathbf{S}^{(r)} : \mathbf{\Gamma}^1 - \mathbf{\Gamma}^2 : \bar{\mathbf{S}}] \qquad (A.14) $$

and the strain localization by

$$ A^{(r)} = [\mathbf{\Gamma}^2 : \mathbf{C}^{(r)} - \mathbf{C}^{(r)} : \mathbf{\Gamma}^1]^{-1} : [\mathbf{C}^{(r)} : \mathbf{\Gamma}^2 - \mathbf{\Gamma}^1 : \bar{\mathbf{C}}] \qquad (A.15) $$

To express the residual stress field, which turns out to be uniform in each phase, one has to remark that (A.13) combined with (A.11) leads to

$$ \sigma = \bar{\sigma} + \mathbf{\Gamma}^2 : (\sigma - \bar{\sigma}), \quad \epsilon = \bar{\epsilon} + \mathbf{\Gamma}^1 : (\epsilon - \bar{\epsilon}). \qquad (A.16) $$

These latter expressions can be plugged into the local behavior (A.1)$_2$ to give

$$ (\mathbf{S}^{(r)} : \mathbf{\Gamma}^2 - \mathbf{\Gamma}^1 : \mathbf{S}^{(r)}) : \sigma^{(r)} = \mathbf{\Gamma}^2 : (\bar{\epsilon} - \epsilon^{(r)}) - \mathbf{S}^{(r)} : \mathbf{\Gamma}^1 : \bar{\sigma}. \qquad (A.17) $$

The residual stress is defined as the field of stress remaining in the material when there is no applied stress, i.e. for $\bar{\sigma} = \mathbf{0}$ leading also to $\bar{\epsilon} = \bar{\epsilon}^0$. One thus gets

$$ \sigma^{\text{res}(r)} = (\mathbf{S}^{(r)} : \mathbf{\Gamma}^2 - \mathbf{\Gamma}^1 : \mathbf{S}^{(r)})^{-1} : \mathbf{\Gamma}^2 : (\bar{\epsilon}^0 - \epsilon^{(r)}) \qquad (A.18) $$

with $\bar{\epsilon}^0$ given by (A.8).

## Data availability

The data to reproduce these findings will be made available upon reasonable request.

## References

[1] K. Aizu, Possible species of ferromagnetic, ferroelectric and ferroelastic crystals, Phys. Rev. B 2 (1970) 754–772, https://doi.org/10.1103/physrevb.2.754.

[2] A. Authier, International Table for Crystallography, vol. D: Physical Properties of Crystals, Kluwer Academic Publishers, 2003.

[3] G.E. Backus, Long-wave elastic anisotropy produced by horizontal layering, J. Geophys. Res. 67 (1962) 4427–4440, https://doi.org/10.1029/jz067i011p04427.

[4] G. Bansal, A. Heuer, Precipitation in partially stabilized zirconia, J. Am. Ceram. Soc. 58 (1975) 235–238, https://doi.org/10.1111/j.1151-2916.1975.tb11453.x.

[5] G. Boittin, H. Moulinec, F. Silva, CRaFT: composite response and Fourier transform, https://lma-software-craft.cnrs.fr/, 2024.

[6] S.N. Chiu, D. Stoyan, W.S. Kendall, J. Mecke, Stochastic Geometry and Its Applications, 3 ed., J. Wiley & Sons, 2013.

[7] B. Das, P. Brodard, P. Bandyopadhyay, Raman spectroscopy assisted residual stress measurement of plasma sprayed and laser remelted zirconia splats and coatings, Surf. Coat. Technol. 378 (2019) 124920, https://doi.org/10.1016/j.surfcoat.2019.124920.

[8] G. Fadda, L. Truskinovsky, G. Zanzotto, Unified Landau description of the tetragonal, orthorhombic, and monoclinic phases of zirconia, Phys. Rev. B 66 (2002) 174107, https://doi.org/10.1103/PhysRevB.66.174107.

[9] G.A. Francfort, F. Murat, Homogenization and optimal bounds in linear elasticity, Arch. Ration. Mech. Anal. 94 (1986) 307–334, https://doi.org/10.1007/BF00280908.

[10] Y. Fu, A. Evans, W. Kriven, Microcrack nucleation in ceramics subjects to a phase transformation, J. Am. Ceram. Soc. 67 (1984) 626–630, https://doi.org/10.1111/j.1151-2916.1984.tb19607.x.

[11] O. Furat, L. Petrich, D. Finegan, D. Diercks, F. Usseglio-Viretta, K. Smith, V. Schmidt, Artificial generation of representative single Li-ion electrode particle architectures from microscopy data, npj Comput. Mater. 7 (2021) 105, https://doi.org/10.1038/s41524-021-00567-9.

[12] M. Groeber, S. Ghosh, M.D. Uchic, D.M. Dimiduk, A framework for automated analysis and simulation of 3D polycrystalline microstructures. Part 2: synthetic structure generation, Acta Mater. 56 (2008) 1274–1287, https://doi.org/10.1016/j.actamat.2007.11.040.

[13] T. Gu, O. Castelnau, S. Forest, E. Hervé-Luanco, F. Lecouturier, H. Proudhon, L. Thilly, Multiscale modeling of the elastic behavior of architected and nanostructured Cu-Nb composite wires, Int. J. Solids Struct. 121 (2017) 148–162, https://doi.org/10.1016/j.ijsolstr.2017.05.022.

[14] R. Guinebretière, T. Ors, V. Michel, E. Thune, M. Huger, S. Arnaud, N. Blanc, N. Boudet, O. Castelnau, Coupling between elastic strains and phase transition in dense pure zirconia polycrystals, Phys. Rev. Mater. 6 (2022) 013602, https://doi.org/10.1103/PhysRevMaterials.6.013602.

[15] M. Hayakawa, N. Kuntani, M. Oka, Structural study on the tetragonal to monoclinic transformation in arc-melted ZrO2-2 mol.%Y2O3 – I. experimental observations, Acta Metall. 37 (1989) 2223–2228, https://doi.org/10.1016/0001-6160(89)90148-X.

[16] M. Hayakawa, M. Tada, H. Okamoto, M. Oka, Structural study on the cubic to tetragonal transformation in arc-melted ZrO2-3 mol.%Y2O3, Trans. Jpn. Inst. Met. 27 (1986) 750–760, https://doi.org/10.2320/matertrans1960.27.750.

[17] M. Idiart, Nonlinear sequential laminates reproducing hollow sphere assemblages, C. R., Méc. 335 (2007) 363–368, https://doi.org/10.1016/j.crme.2007.04.003.

[18] P.M. Kelly, L.R. Francis Rose, The martensitic transformation in ceramics — its role in transformation toughening, Prog. Mater. Sci. 47 (2002) 463–557, https://doi.org/10.1016/S0079-6425(00)00005-0.

[19] E.H. Kisi, C.J. Howard, Elastic constants of tetragonal zirconia measured by a new powder diffraction technique, J. Am. Ceram. Soc. 81 (1998) 1682–1684, https://doi.org/10.1111/j.1151-2916.1998.tb02533.x.
[20] J.A. Krogstad, M. Lepple, Y. Gao, D.M. Lipkin, C.G. Levi, Effect of yttria content on the zirconia unit cell parameters, J. Am. Ceram. Soc. 94 (2011) 4548–4555, https://doi.org/10.1111/j.1551-2916.2011.04862.x.

[21] L.D. Landau, E.M. Lifshitz, Statistical Physics, 3 ed., Elsevier, 2013.

[22] S.M. Lang, Axial thermal expansion of tetragonal ZrO, between 1150° and 1700°C, J. Am. Ceram. Soc. 47 (1964) 641–644, https://doi.org/10.1111/j.1151-2916.1964.tb13125.x.

[23] R.A. Lebensohn, P. Ponte Castañeda, R. Brenner, O. Castelnau, Full-field vs. Homogenization methods to predict microstructure-property relations for polycrystalline materials, in: S. Ghosh, D. Dimiduk (Eds.), Computational Methods for Microstructure-Property Relationships, Springer, 2011, pp. 393–441.

[24] R.A. Lebensohn, A.D. Rollett, Spectral methods for full-field micromechanical modeling of polycrystalline materials, Comput. Mater. Sci. 173 (2020) 109336, https://doi.org/10.1016/j.commatsci.2019.109336.

[25] S. Lucarini, M.V. Upadhyay, J. Segurado, FFT based approaches in micromechanics: fundamentals, methods and applications, Model. Simul. Mater. Sci. Eng. 30 (2021) 023002, https://doi.org/10.1088/1361-651X/ac34e1.

[26] M. Mamivand, M. Asle Zaeem, H. El Kadiri, Phase field modeling of stress-induced tetragonal-to-monoclinic transformation in zirconia and its effect on transformation toughening, Acta Mater. 64 (2014) 208–219, https://doi.org/10.1016/j.actamat.2013.10.031.

[27] M. Mamivand, M. Asle Zaeem, H. El Kadiri, L.Q. Chen, Phase field modeling of the tetragonal-to-monoclinic phase transformation in zirconia, Acta Mater. 61 (2013) 5223–5235, https://doi.org/10.1016/j.actamat.2013.05.015.

[28] G.W. Milton, The Theory of Composites, Cambridge University Press, 2002.

[29] H. Moulinec, P. Suquet, A numerical method for computing the overall response of nonlinear composites with complex microstructure, Comput. Methods Appl. Mech. Eng. 157 (1998) 69–94, https://doi.org/10.1016/S0045-7825(97)00218-1.

[30] B.C. Muddle, R. Hannink, Crystallography of the tetragonal to monoclinic transformation in MgO partially stabilized zirconia, J. Am. Ceram. Soc. 69 (1986) 547–555, https://doi.org/10.1111/j.1151-2916.1986.tb04791.x.

[31] I. Ors, N. Gey, M. Huger, O. Castelnau, R. Guinebretière, Strains, stresses and microcracks in pure zirconia polycrystals, Mater. Charact. 220 (2025) 114667, https://doi.org/10.1016/j.matchar.2024.114667.

[32] T. Ors, F. Gouraud, V. Michel, M. Huger, N. Gey, J.S. Micha, O. Castelnau, R. Guinebretière, Huge local elastic strains in bulk nanostructured pure zirconia materials, Mater. Sci. Eng. A 806 (2021) 140817, https://doi.org/10.1016/j.msea.2021.140817.

[33] T. Ors, F. Gouraud, R. Guinebretière, M. Huger, V. Michel, O. Castelnau, Neutron diffraction measurements of residual stress distribution in large zirconia based refractory bricks produced by electro-fusion and casting, J. Eur. Ceram. Soc. 37 (2017) 2295–2302, https://doi.org/10.1016/j.jeurceramsoc.2017.01.009.

[34] C. Patapy, M. Huger, T. Chotard, R. Guinebretiere, N. Gey, A. Hazotte, M. Humbert, Solidification structure in pure zirconia liquid molten phase, J. Eur. Ceram. Soc. 33 (2013) 259–268, https://doi.org/10.1016/j.jeurceramsoc.2012.08.021.

[35] J. Paux, L. Morin, L. Gélébart, A.M. Amadou Sanoko, A discrete sine–cosine based method for the elasticity of heterogeneous materials with arbitrary boundary conditions, Comput. Methods Appl. Mech. Eng. 433 (2025), https://doi.org/10.1016/j.cma.2024.117488.

[36] Z. Pi, K. Wang, L. Yang, Y. Zhou, On the theoretical and phase field modeling of the stress state associated with ferroelastic twin nucleation and propagation near crack tip, Eng. Fract. Mech. 235 (2020) 107200, https://doi.org/10.1016/j.engfracmech.2020.107200.

[37] Z. Pi, F. Zhang, J. Chen, W. Zhu, L. Yang, Y. Zhou, Multiphase field theory for ferroelastic domain switching with an application to tetragonal zirconia, Comput. Mater. Sci. 170 (2019) 109165, https://doi.org/10.1016/j.commatsci.2019.109165.

[38] E. Plancher, V. Favier, C. Maurice, E. Bosso, N. Rupin, J. Stodolna, D. Loisnard, J.B. Marijon, J. Petit, J.S. Micha, O. Robach, O. Castelnau, Direct measurement of local constitutive relations, at the micrometre scale, in bulk metallic alloys, J. Appl. Crystallogr. 50 (2017) 940–948, https://doi.org/10.1107/S1600576717006185.

[39] R.R.P. Purushottam Raj Purohit, D. Fowan, S. Arnaud, N. Blanc, J.S. Micha, R. Guinebretière, O. Castelnau, Laue microdiffraction on polycrystalline samples above 1500 K achieved with the QMAX-μLaue furnace, J. Appl. Crystallogr. 57 (2024) 470–480, https://doi.org/10.1107/S1600576724001821.

[40] R.R.P. Purushottam Raj Purohit, D.P. Fowan, E. Thune, S. Arnaud, G. Chahine, N. Blanc, O. Castelnau, R. Guinebretière, Phase transition and twinning in polycrystals probed by in situ high temperature 3D reciprocal space mapping, Appl. Phys. Lett. 121 (2022) 181901, https://doi.org/10.1063/5.0109058.

[41] R.R.P. Purushottam Raj Purohit, T. Richeton, S. Berbenni, L. Germain, N. Gey, T. Connelly, O. Castelnau, Estimating single-crystal elastic constants of polycrystalline β metastable titanium alloy: a Bayesian inference analysis based on high energy x-ray diffraction and micromechanical modeling, Acta Mater. 208 (2021) 116762, https://doi.org/10.1016/j.actamat.2021.116762.

[42] R. Quey, L. Renversade, Optimal polyhedral description of 3D polycrystals: method and application to statistical and synchrotron x-ray diffraction data, Comput. Methods Appl. Mech. Eng. 330 (2018) 308–333, https://doi.org/10.1016/j.cma.2017.10.029.

[43] C. Redenbach, A. Liebscher, Random tessellations and their application to the modelling of cellular materials, in: V. Schmidt (Ed.), Stochastic Geometry, Spatial Statistics and Random Fields: Models and Algorithms, Springer, 2015, pp. 73–93.

[44] F. Seitl, L. Petrich, J. Staněk, C.E. Krill III, V. Schmidt, V. Beneš, Exploration of Gibbs-Laguerre tessellations for three-dimensional stochastic modeling, Methodol. Comput. Appl. Probab. 23 (2020) 669–693, https://doi.org/10.1007/s11009-019-09757-x.

[45] P.H. Serrao, S. Sandfeld, A. Prakash, OptiMic: a tool to generate optimized polycrystalline microstructures for materials simulations, SoftwareX 15 (2021) 100708, https://doi.org/10.1016/j.softx.2021.100708.

[46] N. Simha, Twin and habit plane microstructures due to the tetragonal to monoclinic transformation of zirconia, J. Mech. Phys. Solids 45 (1997) 261–292, https://doi.org/10.1016/S0022-5096(96)00074-9.

[47] M. Smirnov, A. Mirgorodsky, R. Guinebretière, Phenomenological theory of lattice dynamics and polymorphism of ZrO₂, Phys. Rev. B 68 (2003) 104106, https://doi.org/10.1103/PhysRevB.68.104106.

[48] P. Suquet, H. Moulinec, O. Castelnau, M. Montagnat, N. Lahellec, F. Grennerat, P. Duval, R. Brenner, Multiscale modeling of the mechanical behavior of polycrystalline ice under transient creep, in: O. Cazacu (Ed.), IUTAM Symp. Linking Scales in Computations: From Microstructure to Macro-Scale Properties, 2012, pp. 76–90.

[49] Y. Wang, H.Y. Wang, L.Q. Chen, A.G. Khachaturyan, Microstructural development of coherent tetragonal precipitates in magnesium-partially-stabilized zirconia: a computer simulation, J. Am. Ceram. Soc. 78 (1995) 657–661, https://doi.org/10.1111/j.1151-2916.1995.tb08228.x.

[50] M. Yashima, N. Ishizawa, M. Yoshimura, High-temperature x-ray study of the cubic-tetragonal diffusionless phase transition in the ZrO2-ErO1.5 system: I, phase change between two forms of a tetragonal phase, t'-ZrO2 and t"-ZrO2, in the compositionally homogeneous 14 mol% ErO1.5-ZrO2, J. Am. Ceram. Soc. 76 (1993) 641–648, https://doi.org/10.1111/j.1151-2916.1993.tb03654.x.

[51] X.S. Zhao, S.L. Shang, Z.K. Liu, J.Y. Shen, Elastic properties of cubic, tetragonal and monoclinic ZrO2 from first-principles calculations, J. Nucl. Mater. 415 (2011) 13–17, https://doi.org/10.1016/j.jnucmat.2011.05.016.

[52] J. Zhu, J. Luo, Study of transformation induced intergranular microcracking in tetragonal zirconia polycrystals with the phase field method, Mater. Sci. Eng. A 701 (2017) 69–84, https://doi.org/10.1016/j.msea.2017.06.060.