# Stress partitioning and effective behavior of N-phase laminates in anisotropic elasticity from a fast explicit method

Thiebaud Richeton

► To cite this version:

Thiebaud Richeton. Stress partitioning and effective behavior of N-phase laminates in anisotropic elasticity from a fast explicit method. Journal of Theoretical, Computational and Applied Mechanics, 2022, Volume 30 (2022), Issue 1, pp.1-14. 10.46298/jtcam.8506. hal-03347969v3

HAL Id: hal-03347969

https://hal.science/hal-03347969v3

Submitted on 31 Mar 2022

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
**Identifiers**

DOI 10.46298/jtcam.8506

OAI hal-03347969v3

**History**

Received Sep 20, 2021

Accepted Feb 14, 2022

Published Mar 30, 2022

**Associate Editor**

Laurence BRASSART

**Reviewers**

Anonymous
Reinaldo RODRÍGUEZ-RAMOS

**Open Review**

OAI hal-03625301

**Supplementary Material**

Software
SWHID [SoftwareHeritage]

**Licence**

CC BY 4.0
©The Authors

# Stress partitioning and effective behavior of N-phase laminates in anisotropic elasticity from a fast explicit method

Thiebaud RICHETON

Université de Lorraine, CNRS, Arts et Métiers ParisTech, LEM3, F-57000 Metz, France

In this work, a fast explicit method, easy to implement numerically, is proposed in order to compute the effective behavior and the distribution of stresses in a general N-phase laminate made of parallel, planar and perfectly bonded interfaces. The solutions are exact for a homogeneous far-field loading and work for an arbitrary number of phases, a general linear anisotropic elasticity, as well as different uniform thermal and plastic strains in the phases. A simple direct analytical formula is also derived to compute the stress in a given phase once the effective behavior of the laminate is known. Moreover, the correctness of the proposed method is checked by comparisons with finite element simulation results on a same boundary value problem, showing excellent agreements. An application of the method is performed for a near-β titanium alloy with elongated grains, by comparing the level of internal stresses for different elastic loadings within a N-phase laminate made of 100 000 orientations and a 2-phase laminate of equal volume fraction with maximal elastic contrast. Interestingly, the maximum von Mises stress of the 2-phase laminate is always the lowest, which is explained by a volume fraction effect. Finally, comparisons with elastic self-consistent models considering oblate spheroidal grains of different aspect ratios are performed.

**Keywords:** laminate, anisotropic elasticity, incompatibility stresses, effective behavior, composite structure

## 1 Introduction

At a material surface of discontinuity like a grain boundary or an interface between two dissimilar materials, incompatibility stresses may arise due to heterogeneous thermo-elastic properties and/or difference of plastic strains. These internal stresses can be computed numerically thanks to finite element methods (Roters et al. 2010) or spectral approaches (A. Lebensohn 2001). Analytical formulations also exist in the case of two-phase laminates (Stupkiewicz and Petryk 2002; Milton 2002; Franciosi and Berbenni 2007; Glüge and Kalisch 2014) which can also be viewed as infinite bi-crystals or bi-materials subjected to a homogeneous far-field loading (Gemperlová and Paidar 1985; Gemperlová et al. 1989; Richeton and Berbenni 2013; Richeton et al. 2015; Richeton 2017). Such analytical models provide instantaneous responses which are first-order estimations around interfaces subjected to more complex boundary conditions (Peralta et al. 1993; Tiba et al. 2015). Hence, it can be very useful to test, on a huge number of configurations, the influence of different parameters like material single crystal elastic constants, crystallographic misorientation, crystal volume fraction, interface inclination, amount of plastic deformation or loading type (Gemperlová et al. 1989; Richeton and Berbenni 2013; Richeton et al. 2015). In the same way, explicit analytical models for N-phase laminates are very useful to investigate the influence of the previously mentioned parameters in problems involving several parallel interfaces, like in thin films (e.g., Welzel et al. 2003; Abadias et al. 2018), in electronics packaging (e.g., Wong and Lim 2008), in nanolaminated materials (e.g., Mukhopadhyay et al. 2017; Wang et al. 2017) or in materials with welding joints which are classically modelled as a base metal, a heat affected zone and a weld metal (e.g., Du et al. 2020). Actually, laminated materials are widely used in the industry because of the relative ease to control their properties from the choice of the materials combination,
volume fraction and orientation (Glüge and Kalisch 2014). For the design of laminates, fast estimates of their properties is thus required.

The objective of the present paper is to derive an explicit method, which can be easily implemented in a numerical code, in order to compute the effective behavior and the distribution of stresses in a general N-phase laminate made of parallel, planar and perfectly bonded interfaces and which is assumed to be subjected to a far-field homogeneous loading and/or homogeneous temperature variation, in addition to the possibility of having different uniform plastic strains in each phase. The method should be exact for a small number of phase (N = 2, 3, 4, ...) as well as for an infinitely large number of phases. Such an objective may appear relatively standard considering previous works on laminates in the literature (e.g., Milton 2002; Franciosi and Berbenni 2007; El Omri et al. 2000; Glüge and Kalisch 2014; Glüge 2016). However, it is generally hard to find compact explicit formulas for direct use in engineering applications. Besides, at the difference of previous approaches where explicit formulas were restricted to a finite number of phases (Milton 2002; Franciosi and Berbenni 2007) or to isotropic elasticity (El Omri et al. 2000; Milton 2002; Glüge 2016) or to purely elastic behavior (Milton 2002; Glüge and Kalisch 2014), this contribution provides explicit compact formulas for an arbitrary number of phases, a completely general anisotropic elasticity as well as different thermal and plastic strains in the phases. In particular, a simple direct analytical formula is derived in order to compute the stress in a given phase once the effective behavior of the laminate has been determined. Furthermore, the code to compute the effective behavior and the stress distribution in an arbitrary N-phase laminate is made freely available (see Supplementary Material at the end of the document).

The paper is organized as follows. Section 2 reminds the expressions of incompatibility stresses in a two-phase laminate from which the effective compliance, effective plastic strain and effective thermal expansion tensors can be deduced. Then, Section 3 presents an explicit method to find the effective behavior and the stress distribution for a general N-phase laminate based on the solutions of the two-phase laminate. Section 4 is dedicated to the validation of the proposed method by comparisons with finite element simulation results on a same boundary value problem. In Section 5, an application of the method is performed for a near-β titanium alloy with elongated grains, by comparing the level of internal stresses for different elastic loadings within a N-phase laminate made of 100 000 orientations and a 2-phase laminate with maximal elastic contrast. Then, in Section 6, comparisons with elastic self-consistent models considering oblate spheroidal grains of different aspect ratios are performed. In particular, the numerical performance of the laminate model compared to the self-consistent models is discussed. Finally, concluding remarks are given in Section 7.

In the following, the Einstein summation convention over repeated indices and the contracted Voigt notation (Voigt 1966) (11 → 1, 22 → 2, 33 → 3, 23 → 4, 31 → 5, 12 → 6) are used, i.e. indices range from 1 to 6. For consistency, an engineering convention is considered for strain components, i.e., $ε_4 = 2ε_{23}$, $ε_5 = 2ε_{31}$, $ε_6 = 2ε_{12}$, while the components of the elastic compliance tensor $s_{ij}$ include the multiplying factors of 2 and 4 (Voigt 1966). Besides, the notation $\langle \rangle$ denotes a volume average and superscripts 1, 2 or m denote fields in materials 1, 2 or m, respectively.

# 2 Two-phase laminate

## 2.1 Averaging rules, constitutive behavior and continuity conditions

In this section, an infinite 2-phase laminate composed of two alternating materials or crystals perfectly bonded along planar interfaces is considered. This rank-1 laminate is supposed to have been deformed under the action of a homogeneous temperature variation $\Delta T$, as well as a macroscopic homogeneous and remotely applied stress $\Sigma_i$ satisfying the averaging rule

$$ \Sigma_i = \langle \sigma_i \rangle, \qquad (1) $$

where $\sigma_i$ is the local Cauchy stress. Linear thermo-elasticity is assumed, as well as the possibility of having different uniform plastic strains $\epsilon_i^p$ in each phase. Under these assumptions and through consideration of strains compatibility and balance of linear momentum without body force in a static small strain setting, it can be shown that $\sigma_i$ and the local strains $\epsilon_i$ are uniform in each
material (Gemperlová et al. 1989; El Omri et al. 2000; Milton 2002; Richeton and Berbenni 2013), so that

$$ \Sigma_i = f^1 \sigma_i^1 + f^2 \sigma_i^2, \qquad (2) $$

and

$$ E_i = \langle \epsilon_i \rangle = f^1 \epsilon_i^1 + f^2 \epsilon_i^2, \qquad (3) $$

where $f^1$ and $f^2 = 1 - f^1$ are the material volume fractions. $E_i$ is defined as the macroscopic strain. The strains in each material are expressed as

$$ \epsilon_i^m = s_{ij}^m \sigma_j^l + \epsilon_i^{p^m} + \alpha_i^{T^m} \Delta T, \qquad (4) $$

where $m = 1, 2$ while $s_{ij}$ and $\alpha_i^T$ are the components of the elastic compliance tensor and of the symmetrical thermal expansion tensor, respectively. Both $s_{ij}$ and $\alpha_i^T$ are supposed isotherm and uniform in each material. Besides, considering that the interface normals are parallel to the direction $e_2$ of a Cartesian frame ($e_1, e_2, e_3$), Equations (2) and (3) along with tractions continuity and strains compatibility impose the following continuity conditions (El Omri et al. 2000; Richeton and Berbenni 2013; Glüge and Kalisch 2014)

$$ \left\{ \begin{array}{l} \sigma_2^1 = \sigma_2^2 = \Sigma_2, \\ \sigma_4^1 = \sigma_4^2 = \Sigma_4, \\ \sigma_6^1 = \sigma_6^2 = \Sigma_6, \end{array} \right. \quad \text{and} \quad \left\{ \begin{array}{l} \epsilon_1^1 = \epsilon_1^2 = E_1, \\ \epsilon_3^1 = \epsilon_3^2 = E_3, \\ \epsilon_5^1 = \epsilon_5^2 = E_5. \end{array} \right. \qquad (5) $$

## 2.2 Stress partitioning

It was previously shown that the stresses in both materials can be written as (Richeton and Berbenni 2013; Richeton et al. 2015; Tiba et al. 2015; Richeton 2017)

$$ \begin{aligned} \sigma_i^1 &= \Sigma_i + f^2 G_{ik} [(s_{kj}^2 - s_{kj}^1)\Sigma_j + \epsilon_k^{p^2} - \epsilon_k^{p^1} + (\alpha_k^{T^2} - \alpha_k^{T^1})\Delta T], \\ \sigma_i^2 &= \Sigma_i - f^1 G_{ik} [(s_{kj}^2 - s_{kj}^1)\Sigma_j + \epsilon_k^{p^2} - \epsilon_k^{p^1} + (\alpha_k^{T^2} - \alpha_k^{T^1})\Delta T], \end{aligned} \qquad (6) $$

where the non-zero components of the symmetric tensor $G_{ij}$ are given by

$$ \begin{aligned} G_{11} &= (\tilde{s}_{35}^2 - \tilde{s}_{33}\tilde{s}_{55})/D, & G_{13} &= (\tilde{s}_{13}\tilde{s}_{55} - \tilde{s}_{15}\tilde{s}_{35})/D, \\ G_{33} &= (\tilde{s}_{15}^2 - \tilde{s}_{11}\tilde{s}_{55})/D, & G_{15} &= (\tilde{s}_{15}\tilde{s}_{33} - \tilde{s}_{13}\tilde{s}_{35})/D, \\ G_{55} &= (\tilde{s}_{13}^2 - \tilde{s}_{11}\tilde{s}_{33})/D, & G_{35} &= (\tilde{s}_{11}\tilde{s}_{35} - \tilde{s}_{13}\tilde{s}_{15})/D, \end{aligned} \qquad (7) $$

with $D = \tilde{s}_{11}\tilde{s}_{35}^2 + \tilde{s}_{33}\tilde{s}_{15}^2 + \tilde{s}_{55}\tilde{s}_{13}^2 - \tilde{s}_{11}\tilde{s}_{33}\tilde{s}_{55} - 2\tilde{s}_{13}\tilde{s}_{15}\tilde{s}_{35}$ and $\tilde{s}_{ij} = f^2 s_{ij}^1 + f^1 s_{ij}^2$.

It is noteworthy that these expressions can also be viewed as the solutions in an infinite bi-material (Gemperlová et al. 1989; Richeton and Berbenni 2013; Richeton et al. 2015) and that the effect grain boundary sliding could be considered as well (Richeton 2017). The differences $\sigma_i^1 - \Sigma_i$ and $\sigma_i^2 - \Sigma_i$ are called incompatibility stresses in materials 1 and 2, respectively. Besides, it must be noted that the tensor $G_{ij}$ is a function of $s_{ij}^1, s_{ij}^2$ and $f^1$ (or $f^2$) only so that we will use the notation $G_{ij}(s^1, s^2, f^1)$ in Section 3 where multiple tensors $G_{ij}$ are considered.

## 2.3 Effective behavior

By definition of the effective behavior, the macroscopic strain satisfies the following relation

$$ E_i = S_{ij}\Sigma_j + E_i^p + A_i^T \Delta T, \qquad (8) $$

where $S_{ij}$ is the effective compliance tensor, $E_i^p$ the effective plastic strain tensor and $A_i^T$ the effective thermal expansion tensor. From Equations (3), (4) and (6), we have also

$$ \begin{aligned} E_i &= \langle s_{ij} \sigma_j \rangle + \langle \epsilon_i^p \rangle + \langle \alpha_i^T \rangle \Delta T \\ &= \langle s_{ij} \rangle \Sigma_j - f^1 f^2 (s_{ij}^2 - s_{ij}^1) G_{jk} ((s_{kl}^2 - s_{kl}^1) \Sigma_l + \epsilon_k^{p^2} - \epsilon_k^{p^1} + (\alpha_k^{T^2} - \alpha_k^{T^1}) \Delta T) + \langle \epsilon_i^p \rangle + \langle \alpha_i^T \rangle \Delta T. \end{aligned} \qquad (9) $$
Then, by identification with Equation (8), the three effective tensors can be expressed as

$$S_{ij} = \langle s_{ij} \rangle - f^1 f^2 (s_{ik}^2 - s_{ik}^1) G_{kl} (s_{lj}^2 - s_{lj}^1) \quad (10)$$

$$E_i^p = \langle \epsilon_i^p \rangle - f^1 f^2 (s_{ik}^2 - s_{ik}^1) G_{kl} (\epsilon_l^{p2} - \epsilon_l^{p1}) \quad (11)$$

$$A_i^T = \langle \alpha_i^T \rangle - f^1 f^2 (s_{ik}^2 - s_{ik}^1) G_{kl} (\alpha_l^{T2} - \alpha_l^{T1}). \quad (12)$$

# 3 N-phase laminate

In the following, an infinite N-phase laminate composed of N alternating materials perfectly bonded along parallel planar interfaces is considered. All the other assumptions made for the 2-phase laminate in Section 2, including Equation (1) and Equation (4), are maintained so that stresses and strains are uniform within each phase (Glüge and Kalisch 2014).

## 3.1 Effective behavior

The solutions of the 2-phase laminate are first used to find the effective behavior of the N-phase laminate and then the stresses in each phase. The general methodology is described in Figure 1. It

<figure><img src="image_2.png" /><figcaption>**Figure 1** Description of the multiple steps procedure for determining the effective properties of a N-phase laminate and the stresses in each phase: $\ell$ denotes the step number and $K^\ell$ the number of phases at step $\ell$. In every phase, the properties are denoted $prop^{\ell,m} = (f^{\ell,m}, s^{\ell,m}, \epsilon^{\ell,m}, \alpha^{\ell,m})$.</figcaption></figure>

consists in a multiple steps procedure, where at each step $\ell$, the phases are grouped two by two and replaced by a new phase having the effective properties of a 2-phase laminate made of the corresponding two phases. The procedure is followed until only one phase is remaining, corresponding to

$$\ell = n = \lceil \log_2 N + 1 \rceil. \quad (13)$$

The number of phases or materials of the laminate at a step $\ell$ is denoted $K^\ell$. Hence, we have $K^\ell = N$, $K^n = 1$ and

$$K^{\ell+1} = \begin{cases} K^\ell/2 & \text{if } K^\ell \text{ even} \\ (K^\ell + 1)/2 & \text{otherwise.} \end{cases} \quad (14)$$

At step $\ell + 1$, the properties of a new phase $m$, i.e. the volume fractions $f^{\ell+1,m}$, elastic compliances $s_{ij}^{\ell+1,m}$, plastic strains $\epsilon_i^{\ell+1,m}$ and thermal expansions $\alpha_i^{\ell+1,m}$, are computed based on Equations (10)
to (12):

$$
\begin{align}
f^{\ell+1,m} &= f^{\ell,2m-1} + f^{\ell,2m} \nonumber \\
s_{ij}^{\ell+1,m} &= f^* s_{ij}^{\ell,2m-1} + (1-f^*) s_{ij}^{\ell,2m} - f^*(1-f^*)(s_{ik}^{\ell,2m} - s_{ik}^{\ell,2m-1}) G_{k\ell}(s^{2m-1}, s^{2m}, f^*) (s_{\ell j}^{\ell,2m} - s_{\ell j}^{\ell,2m-1}) \tag{15} \\
\epsilon_i^{p^{\ell+1,m}} &= f^* \epsilon_i^{p^{\ell,2m-1}} + (1-f^*) \epsilon_i^{p^{\ell,2m}} - f^*(1-f^*)(s_{ik}^{\ell,2m} - s_{ik}^{\ell,2m-1}) G_{k\ell}(s^{2m-1}, s^{2m}, f^*) (\epsilon_\ell^{p^{\ell,2m}} - \epsilon_\ell^{p^{\ell,2m-1}}) \nonumber \\
\alpha_i^{T^{\ell+1,m}} &= f^* \alpha_i^{T^{\ell,2m-1}} + (1-f^*) \alpha_i^{T^{\ell,2m}} - f^*(1-f^*)(s_{ik}^{\ell,2m} - s_{ik}^{\ell,2m-1}) G_{k\ell}(s^{2m-1}, s^{2m}, f^*) (\alpha_\ell^{T^{\ell,2m}} - \alpha_\ell^{T^{\ell,2m-1}}), \nonumber
\end{align}
$$

with $f^* = f^{\ell,2m-1}/(f^{\ell,2m-1} + f^{\ell,2m})$. By imposing the consistency condition

$$
f^* \sigma_i^{\ell,2m-1} + (1-f^*) \sigma_i^{\ell,2m} = \sigma_i^{\ell+1,m} \quad (16)
$$

between stresses at step $\ell$ and at step $\ell + 1$, the stresses at step $\ell$ can be related to the stresses at step $\ell + 1$ by application of Equation (6):

$$
\begin{align}
\sigma_i^{\ell,2m-1} &= \sigma_i^{\ell+1,m} + (1-f^*) G_{ik} (s^{2m-1}, s^{2m}, f^*) \beta \nonumber \\
\sigma_i^{\ell,2m-1} &= \sigma_i^{\ell+1,m} - f^* G_{ik} (s^{2m-1}, s^{2m}, f^*) \beta \tag{17}
\end{align}
$$

with $\beta = (s_{kj}^{\ell,2m} - s_{kj}^{\ell,2m-1})\sigma_j^{\ell+1,m} + \epsilon_k^{p^{\ell,2m}} - \epsilon_k^{p^{\ell,2m-1}} + (\alpha_k^{T^{\ell,2m}} - \alpha_k^{T^{\ell,2m-1}})\Delta T$. In case where $K^\ell$ is odd, the last phase of step $\ell + 1$ is simply set to be the same as the last phase of step $\ell$, i.e.

$$
f^{\ell+1,K^{\ell+1}} = f^{\ell,K^{\ell}}, \quad s_{ij}^{\ell+1,K^{\ell+1}} = s_{ij}^{\ell,K^{\ell}}, \quad \epsilon_i^{p^{\ell+1,K^{\ell}}} = \epsilon_i^{p^{\ell,K^{\ell}}}, \quad \alpha_i^{T^{\ell+1,K^{\ell+1}}} = \alpha_i^{T^{\ell,K^{\ell}}}, \quad \sigma_i^{\ell+1,K^{\ell+1}} = \sigma_i^{\ell,K^{\ell}}. \quad (18)
$$

From the procedure just described, it is noticeable that averaging rules and continuity conditions are automatically satisfied at each step, i.e.

$$
\forall \ell, \quad \sum_m f^{\ell,m} = 1, \quad \langle \sigma_i^\ell \rangle = \Sigma_i, \quad \langle \epsilon_i^\ell \rangle = E_i \tag{19}
$$

together with

$$
\forall \ell, m, \quad \sigma_2^{\ell,m} = \Sigma_2, \quad \sigma_4^{\ell,m} = \Sigma_4, \quad \sigma_6^{\ell,m} = \Sigma_6, \tag{20}
$$

and

$$
\forall \ell, m, \quad \epsilon_1^{\ell,m} = E_1, \quad \epsilon_3^{\ell,m} = E_3, \quad \epsilon_5^{\ell,m} = E_5. \tag{21}
$$

As a consequence, the remaining phase at step $n$ satisfies

$$
E_i = s_{ij}^{n,1} \Sigma_j + \epsilon_i^{p^{n,1}} + \alpha_i^{T^{n,1}} \Delta T \quad (22)
$$

which means, from Equation (8), that the effective tensors of the N-phase laminate are

$$
S_{ij} = s_{ij}^{n,1}, \quad E_i^p = \epsilon_i^{p^{n,1}}, \quad \text{and} \quad A_i^T = \alpha_i^{T^{n,1}}. \tag{23}
$$

In the above procedure which consists in grouping the phases in pair iteratively until one phase with the effective properties remains, the initial relative positions of the different phases have obviously no influence on the values obtained at the final step. This means that the effective properties of the laminate are not sensitive to the stacking order of the phases as it was already noted by Glüge and Kalisch (2014).

3.2 Stress partitioning

The stresses in each phase of the initial N-phase laminate can be determined from a top-down algorithm, starting from $\sigma_i^{n,1} = \Sigma_i$ at step $n$ and then going down to step 1 by application of Equation (17) step by step. However, a direct and much faster computation is possible by first noticing that for a 2-phase laminate

$$
\begin{align}
& \Sigma_i + (1 - f^1) G_{ik} (s^1, s^2, f^1) [(s_{kj}^2 - s_{kj}^1) \Sigma_j + \epsilon_k^{p2} - \epsilon_k^{p1} + (\alpha_k^{T2} - \alpha_k^{T1}) \Delta T] \nonumber \\
& = \Sigma_i + G_{ik} (s^1, S, 0) [(S_{kj} - s_{kj}^1) \Sigma_j + E_k^p - \epsilon_k^{p1} + (A_k^T - \alpha_k^{T1}) \Delta T], \tag{24}
\end{align}
$$
where $S_{ij}$, $E_i^p$, $A_i^T$ are, in this equation, the effective properties of the 2-phase laminate. Indeed, let us consider a fictitious cut of phase 1 parallel to the planar interface, such that one of two remaining parts, denoted hereafter 1*, has a negligible volume fraction ($f^{1*} \to 0$). Stresses in phase 1* are equal to those in the first phase due to uniformity of stresses within a phase (see Section 2). This corresponds to the first member of Equation (24) (cf. Equation (6)). Then, as a consequence of the preceding procedure, it is known that these stresses can also be computed from a 2-phase laminate composed of phase 1* and another phase of volume fraction equals to 1 which has thus necessarily the effective properties of the initial 2-phase laminate. This corresponds to the second member of Equation (24). The principle is illustrated in Figure 2.

<table><caption>**Figure 2** Description of the two configurations used to establish Equation (24).</caption><thead><tr><th></th><th></th></tr></thead><tbody><tr><td>$f^1, s^1, \epsilon^p, \alpha^T, \sigma^1$</td><td>$f \to 0, s^1, \epsilon^p, \alpha^T, \sigma^1$</td></tr><tr><td>$f^2, s^2, \epsilon^p, \alpha^T, \sigma^2$</td><td>$f = 1, s, E^p, A^T, \Sigma$</td></tr></tbody></table>

Such a thought experiment can of course also be made for the N-phase laminate, i.e. performing a fictitious cut in a given phase *m* such that one of the remaining part has a negligible volume fraction and then considering a 2-phase laminate made of the phase *m* with negligible volume fraction and another phase having the effective properties of the N-phase laminate. Accordingly, once the effective properties ($S_{ij}$, $E_i^p$, $A_i^T$) of the N-phase laminate are known, the stresses in each phase can be directly computed as

$$ \sigma_i^m = \Sigma_i + G_{ik} \overline{(s^m, S, 0)} [(S_{kj} - s_{kj}^m)\Sigma_j + E_k^p - \epsilon_k^{p^m} + (A_k^T - \alpha_k^{T^m})\Delta T]. \quad (25) $$

Equation (25) shows that the stresses in one phase can be expressed only with respect to the effective properties of the laminate and the properties of the phase. Since the effective properties of the laminate are independent of the relative positions of the phases (cf. Section 3.1), this means that the stresses in the phases are also not sensitive to the stacking order of the phases. Besides, the numerical implementation of the proposed method is very straightforward and is provided as a Matlab code (see Supplementary materials). It contains two functions, one that computes the tensor $G_{ij}$ (cf. Equation (7)) from two elastic compliance tensors and a volume fraction value and one that computes the effective properties of a 2-phase laminate according to Equation (15). The main code is just composed of a single iterative loop that calls the two aforementioned functions in order to compute the properties of the new phases step by step as described in Figure 1. At step *n*, the effective behavior of the N-phase laminate is known through Equation (23) and the stresses in any phase can then be directly computed from Equation (25).

Finally, it may be worth to indicate also the simple expressions of the stresses in the particular case of isotropic homogeneous thermo-elasticity. With $\mu$ the shear modulus and $\nu$ the Poisson's ratio, the plane components are indeed given by (Rey and Zaoui 1980; Richeton and Berbenni 2013)

$$ \begin{aligned} \sigma_1^m &= \Sigma_1 + \frac{2\mu}{1-\nu} [\langle \epsilon_1^p \rangle - \epsilon_1^{p^m} + \nu (\langle \epsilon_3^p \rangle - \epsilon_3^{p^m})], \\ \sigma_3^m &= \Sigma_3 + \frac{2\mu}{1-\nu} [\langle \epsilon_3^p \rangle - \epsilon_3^{p^m} + \nu (\langle \epsilon_1^p \rangle - \epsilon_1^{p^m})], \\ \sigma_5^m &= \Sigma_5 + \mu [\langle \epsilon_5^p \rangle - \epsilon_5^{p^m}], \end{aligned} \quad (26) $$

while the anti-plane components are given by Equation (20).

## 4 Validation by finite element simulations

Several Finite Element (FE) simulations were performed in order to check the correctness of the methodology described in Section 3. The present Section 4 describes the particular case
of an infinite 5-phase laminate with interface normals parallel to the direction **e**₂ and volume
fractions f¹ = 0.25, f² = 0.1, f³ = 0.3, f⁴ = 0.2 and f⁵ = 0.15. Cubic elastic constants of Cu at
room temperature were considered with c₁₁ = 170 GPa, c₁₂ = 124 GPa and c₄₄ = 75 GPa for the
five phases. Crystallographic orientations of the phases were randomly selected. For each phase,
a slip of magnitude 0.02 was assigned to the most stressed system within a pure elastic loading,
i.e. among the 12 FCC slip systems {111}⟨110⟩, the system with the highest absolute value of its
resolved shear stress |τ| was selected and affected γ = ±0.02 so that τγ > 0. No thermal strain
was considered.

The FE simulation was performed thanks to the use of the open source software suite freefem++ (Hecht 2012). A structured mesh made of 270 000 triangular quadratic Lagrangian elements was used for the meshing of a bar of length 100 (see Figure 3), which corresponds to 4 nodes along **e**₁ and **e**₃ and 5001 nodes along **e**₂. Periodic boundary conditions were imposed along the lateral sides. Arbitrary and fixed velocities were imposed on the bottom and upper faces of the bar (u₁ = -0.05, u₂ = 0, u₃ = -0.1 on the bottom face, u₁ = 0.05, u₂ = 0.1, u₃ = 0.1 on the top face). Phase 1 was affected to the region 75 < x₂ < 100, phase 2 to 65 < x₂ < 75, phase 3 to 35 < x₂ < 65, phase 4 to 15 < x₂ < 35 and phase 5 to 0 < x₂ < 15. A first simulation was performed in pure elasticity in order to determine the most stressed systems and then a second simulation was performed considering static plastic strains in each phase.

The distribution of stresses obtained by the FE simulation is shown in Figure 3. For comparison,

<figure><img src="image_3.png" /><figcaption>**Figure 3** Distribution of σ₁ (left), σ₃ (middle) and σ₅ (right) in MPa obtained from FE simulation.</figcaption></figure>

stresses were then computed from Equations (15), (23) and (25) considering the same elastic
constants, the same crystallographic orientations of the phases, the same plastic strains and
taking Σᵢ as the average stress vector given by the FE simulation. In Figure 4, it can be seen
that the matching between the stresses obtained numerically by a FE simulation and the ones
computed analytically from Equations (15), (23) and (25) is excellent. The tiny differences are due

<figure><img src="image_4.png" /><figcaption>**Figure 4** Comparisons of stresses: FE simulation (lines) and computed from Equations (15), (23) and (25) (circles).</figcaption></figure>
to the fact that the FE method is an approximate numerical method. In particular, it was checked
that the stress differences decrease with the number of nodes along **e**₂. These results, along with
the several other FE simulations performed, prove the correctness of the analytical formulas
established in Section 3.

5 Application to polycrystals with elongated grains

As an illustration of potential applications, the N-phase laminate model can be used as a quick
way to estimate the level of internal stresses in metals and alloys having very elongated grains
due to rolling or forging processes. For instance, near-β titanium alloys like Ti-1023 used in
forged pieces of aircraft for landing gears or rotor systems can achieve high specific strength
thanks to the building-up of complex α/β microstructures. However, these Ti-1023 alloys still
contain very large domains of close β orientations which corresponds to the prior-β grains
that were deformed during the forging steps above the β transus (Lhadi et al. 2018; Lhadi et al.
2020). These prior-β grains are millimeter size and are highly elongated along a same direction
corresponding to the axial axis of the billet. Due to the significant elastic anisotropy of the
cubic β phase (Purushottam raj purohit et al. 2021), strong mechanical contrasts exist between
these millimeter size regions which might be at the origin of the early crack initiation sometimes
observed in those materials (Lhadi et al. 2018; Lhadi et al. 2020). In order to get a very first
insight into such a mechanical issue, it can thus be interesting to simplify the microstructure
as an infinite laminate composed of different β orientations. By doing so, the focus is only on
the mechanical interactions between prior-β grains and on a same kind of interfaces, i.e. the
predominant long ones which are parallel to the axial axis of the billet.

Hence, a laminate composed of 100 000 different $\beta$ orientations of equal volume fraction was considered as an application of the methodology described in Section 3. The orientations were generated from the open-source software package Neper (Quey et al. 2018) so that their distribution was nearly uniform in presence of cubic symmetry. It is noteworthy that the effect of a crystallographic texture could have considered as well without difficulty. Elastic constants of Ti-1023, $c_{11} = 92.6$ GPa, $c_{12} = 82.5$ GPa and $c_{44} = 43.5$ GPa, were taken from Purushottam raj purohit et al. (2021) who provided an estimation based on a method coupling Bayesian inference analysis, high energy X-ray diffraction and elastic self-consistent modeling. These constant values correspond to a strong elastic anisotropy ratio of $A = 8.6$. The laminate interfaces normal was set parallel to the direction $\mathbf{e}_2$ of the Cartesian frame ($\mathbf{e}_1, \mathbf{e}_2, \mathbf{e}_3$). A uniaxial macroscopic stress of magnitude 300 MPa was applied without consideration of any plastic or thermal strains. The direction of the uniaxial stress was varied by rotation around $\mathbf{e}_3$. This rotation was described by an angle $\theta$, $\theta = 0^\circ$ corresponding to an uniaxial stress along $\mathbf{e}_1$ and $\theta = 90^\circ$ along $\mathbf{e}_2$. The stresses in each grain were computed from Equations (15), (23) and (25) and then the maximum von Mises stress among the 100 000 grains was recorded and plotted with respect to $\theta$, see Figure 5. For comparison, a 2-phase laminate made of grains of equal volume fraction was considered as well. The cristallographic orientations of the 2 grains were set so that the directions [111] and [1\bar{1}0] in one grain and [100] and [010] in the other grain were parallel to $\mathbf{e}_1$ and $\mathbf{e}_2$, respectively. $\langle 111 \rangle$ and $\langle 100 \rangle$ are actually the directions of maximum and minimum directional Young's modulus, respectively. Hence, the mechanical contrast should very strong for this 2-grain configuration.

In Figure 5, it is interesting to observe that the maximum von Mises stress of the 2-phase
laminate is always below the one with 100 000 orientations despite the strong elastic contrast
considered. However, it was shown in Section 3 that the stresses in a given phase of a N-phase
laminate can also be deduced from a 2-phase laminate where the given phase is in interaction
with the whole laminate (see Figure 2 and Equation (25)). Accordingly, the elastic contrast
between the phase where the maximum von Mises stress is located and the effective laminate
of 100 000 grains should normally not exceed the elastic contrast that exists between [111] $|| \mathbf{e}_1$
and [100] $|| \mathbf{e}_1$-oriented grains. The difference of maximum von Mises stress is actually rather
due to a volume fraction effect, incompatibility stresses in a phase of volume fraction $f$ scaling
indeed as $1-f$, see Equation (6). Yet, the grain volume fraction is $f = 0.5$ in the 2-phase laminate
whereas it is negligible ($f = 10^{-5}$) in the laminate with 100 000 orientations. Similar results could
<figure><img src="image_5.png" /><figcaption>**Figure 5** Comparisons of the maximum von Mises stresses obtained between a laminate made of 100 000 uniform orientations (dashed line) and a 2-phase laminate made of [111] $|| e_1 $ and [100] $|| e_1 $-oriented grains (solid line). $\theta$ represents a rotation around $e_3$ of the uniaxial stress direction. Cubic elastic constants of Ti-1023 were used (Purushottam raj purohit et al. 2021).</figcaption></figure>

also be obtained by considering strong mechanical contrasts resulting from difference of plastic or thermal strains. As a consequence, it can be inferred that considering a 2-phase laminate of equal volume fraction with extremal mechanical contrast may lead to underestimate the maximal level of internal stresses in polycrystals. Furthermore, the difference of positions of the local extrema between the 2-phase laminate and the one with 100 000 orientations can also be noticed in Figure 5.

## 6 Comparisons with elastic self-consistent models and discussion

It is noteworthy that the internal stress evaluation performed in Section 5 could also have been made from an elastic self-consistent model (Hershey 1954; Kröner 1958; Hill 1965a) with grains more realistic modeled as ellipsoidal inclusions. In a self-consistent model, the effective compliance tensor can be computed as

$$S_{ij} = \langle s_{ik}B_{k\ell} \rangle \langle B_{\ell j} \rangle^{-1}, \qquad (27)$$

where the stress concentration tensor $B_{ij}^m$ for any phase $m$ is given by

$$B_{ij}^m = (s_{ik}^m + \tilde{M}_{ik})^{-1} (S_{kj} + \tilde{M}_{kj}), \qquad (28)$$

and where $\tilde{M}_{ij}$ is the constraint tensor defined by

$$\tilde{M}_{ij} = (\delta_{ik} - S_{ik}^E)^{-1} S_{k\ell}^E S_{\ell j}. \qquad (29)$$

$\delta_{ij}$ is the Kronecker delta and $S_{ij}^E$ is the interior Eshelby tensor (Eshelby 1957), here expressed as a 6 × 6 matrix consistent with the used Voigt notation (Barnett and Cai 2018). The stress tensor in a given phase is then directly deduced from the stress concentration tensor

$$\sigma_i^m = B_{ij}^m \Sigma_j. \qquad (30)$$

For ease of direct comparisons with the results of the laminate model, this set of equations was also implemented in the Matlab software. The results of Figure 5 were then compared with the maximum von Mises stresses obtained from elastic self-consistent models considering oblate spheroidal grains ($a = b > c$) of different aspect ratios (Figure 6), the small axis of the oblate spheroids being parallel to the interfaces normal in the laminate model. A convergence towards the laminate solution is observed by increasing the aspect ratio $a/c$ of the oblate spheroids, both for the case with 100 000 orientations and for the specific case with two orientations only. For the case with 100 000 orientations, it is seen that, as long as elongated grains are considered,
<figure><img src="image_6.png" /><figcaption>**Figure 6** Comparisons of the maximum von Mises stress obtained with the laminate model and elastic self-consistent models (SC) considering oblate spheroidal grains ($a = b > c$) of different aspect ratios ($a/c$). Top: case of the two orientations considered in Figure 5. Bottom: case of 100 000 uniform orientations. $\theta$ represents a rotation around $e_3$ of the uniaxial stress direction, see Section 5. Cubic elastic constants of Ti-1023 were used (Purushottam raj purohit et al. 2021).</figcaption></figure>

the laminate model provides fairly good qualitative trends for the variation of the maximal von Mises stress with a rotation of the uniaxial macroscopic stress. However, the maximal von Mises stresses are always overestimated compared to the values provided by self-consistent models with more realistic grain aspect ratios. These over estimations represent about 20 % of the uniaxial macroscopic stress for $a/c = 10$ and about 30 % for $a/c = 5$.

Besides, it must be underlined that Equation (27), from which is computed the effective compliance tensor $S_{ij}$ in the self-consistent model, is an implicit equation since both the stress concentration tensor and the Eshelby tensor depend on $S_{ij}$, see Equations (28) and (29). In the present work, this implicit equation was solved through a fixed-point iteration method (Walpole 1969). This might however cause convergence issues for strong mechanical contrasts. By opposition, the laminate methodology of Section 3 is fully explicit and thus, is expected to be more robust and faster. Moreover, the Eshelby tensor $S_{ij}^E$ in an anisotropic medium cannot be expressed analytically and must be computed numerically. For this work, the Eshelby tensor was computed thanks to the Matlab code provided by Cai 2018 where $S_{ij}^E$ is obtained from the Hill P tensor (Hill 1965b) which is itself computed from numerical integrals that are performed using the quadv function in Matlab. On the contrary, the expression of the tensor $G_{ij}$ in Equation (7) is fully analytical and hence very simple to code. Actually, it can be shown that the Eshelby stress solutions can be directly obtained from the tensor $G_{ij}$ by appropriate integrations (Richeton and Berbenni 2014).

In order to illustrate the numerical differences between the two types of model, the effect of the number of phases on the computation of the effective stiffness tensor, as well as the effect of the absolute tolerance value set in the quadv function in Matlab to evaluate the Eshelby tensor was analyzed. First, Figure 7 shows the relative difference of the norm of the effective stiffness tensor of 100 000 uniform orientations computed by the laminate model and an elastic self-consistent model considering oblate spheroids with $a = b = 100c$ as a function of the absolute tolerance set in the quadv function. It is known that, the more the grains are elongated, the more
<figure><img src="image_7.png" /></figure>

**Figure 7** Black '+' (left axis): relative difference of the norm of the effective stiffness tensor of 100 000 uniform orientations computed by the laminate model and an elastic self-consistent model considering oblate spheroidal grains with $a = b = 100c$ as a function of the absolute tolerance set in the quadv function in Matlab to compute the Eshelby tensor. Red '×' (right axis): corresponding ratios between the CPU time needed by the elastic self-consistent to achieve the calculation over the one needed by the laminate model. Simulations were repeated 10 times to provide average CPU time values. The Frobenius norm of a matrix was considered. The cubic elastic constants of Ti-1023 were used (Purushottam raj purohit et al. 2021).

accuracy is needed on the numerical evaluation of the Eshelby tensor (Gavazzi and Lagoudas 1990). For the present extreme case with $a = b = 100c$, Figure 7 shows that an absolute tolerance lower than $10^{-12}$ should be set in order to obtain a convergence towards an accurate result. At the same time, Figure 7 exhibits also the average ratio between the CPU time needed by the elastic self-consistent to achieve the calculation over the one needed by the laminate model. This ratio is about 7 showing, as expected, that the laminate is much faster. Moreover, this ratio increases exponentially for absolute tolerance smaller than $10^{-13}$. This means that, when using a self-consistent model, a preliminary numerical study is needed to fix the value of the absolute tolerance (or the number of Gauss points if another numerical method is used as in (Gavazzi and Lagoudas 1990)) in order to both ensure a sufficient accuracy on the evaluation of the Eshelby tensor and also to maintain a reasonable CPU time. Then, Figure 8 shows that this relative difference of the norm of the effective stiffness tensor decreases with the number of phases considered. However, it must be underlined that the relative difference is already pretty small (i.e.

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

**Figure 8** Relative difference of the norm of the effective stiffness tensor computed by the laminate model and an elastic self-consistent model considering oblate spheroidal grains with $a = b = 100c$ as a function of the number of phases considered. Phases correspond to crystallographic orientations randomly picked within a list of 100 000 uniform orientations. Simulations were repeated 1000 times to provide average values. The Frobenius norm of a matrix was considered. The absolute tolerance in the quadv function in Matlab was set to $10^{-12}$ to evaluate the Eshelby tensor in the self-consistent model. The cubic elastic constants of Ti-1023 were used (Purushottam raj purohit et al. 2021).
smaller than $4 \times 10^{-5}$) when considering only two phases. This small difference agrees with the convergence towards the laminate solution displayed in Figure 6 when the aspect ratio was increased in the case considering two orientations only. This result can actually be understood from Equation (24). This equation shows indeed that the stresses within a phase of a two-phase laminate are the same as those in a phase of negligible volume fraction and same properties belonging to a two-phase laminate where the other phase has the effective properties of the first two-phase laminate. In the same way, the self-consistent model computes the stresses as a result of the mechanical interactions between one grain and the infinite homogeneous equivalent medium made of the two grains.

# 7 Conclusion

This work provides a fast explicit method to compute the effective behavior and the distribution of stresses in a general N-phase laminate made of parallel, planar and perfectly bonded interfaces. The formulas are based on the solutions of the two-phase laminate. They are very compact which is convenient for numerical implementation and direct use in engineering applications. In contrast with previous approaches, this contribution considers, at the same time, an arbitrary number of phases, a completely general anisotropic elasticity as well as different thermal and plastic strains in the phases. In addition, a simple direct analytical formula is derived in order to compute the stress in a given phase once the effective behavior of the laminate has been determined. The developed method proves that both the effective properties of the laminate and the stresses in the phases are insensitive to the stacking order of the phases. The code to compute the effective behavior and the stress distribution in an arbitrary N-phase laminate is made freely available for download (link in the Supplementary materials section).

Besides, the correctness of the proposed method is checked by comparisons with finite element simulation results on a same boundary value problem, showing excellent agreements. An application of the method is also performed for a near-$\beta$ titanium alloy with elongated grains, by comparing, for different elastic loadings, the level of internal stresses within a N-phase laminate made of 100 000 orientations and a 2-phase laminate with maximal elastic contrast. Interestingly, the maximum von Mises stress of the 2-phase laminate is always the lowest despite the strong elastic contrast considered, which is explained by a volume fraction effect. Similar results could be obtained by considering strong mechanical contrasts resulting from difference of plastic or thermal strains.

Finally, comparisons with elastic self-consistent models considering oblate spheroidal grains of different aspect ratios are performed. A convergence towards the laminate solution is observed by increasing the aspect ratio of the oblate spheroids. However, it is noteworthy that the laminate model is based only on fully analytical and explicit equations, which is more robust from a numerical point of view. It is also shown that the execution of the laminate model is much faster than that of self-consistent models.

As a perspective, it may be worth underlying that the developed method might be extended to obtain the effective strain energy potential of a N-phase viscoplastic laminate, i.e. grouping the N phases in pairs and applying the existing compact expressions for the effective strain energy potential of a two-phase laminate (Debotton 2005; Idiart 2008) for each pair iteratively until one phase with the effective strain energy potential is obtained. The strains in each phase could then be deduced from a top-down algorithm.

# References

A. Lebensohn, R. (2001). N-site modeling of a 3D viscoplastic polycrystal using Fast Fourier Transform. *Acta Materialia* 49(14):2723-2737. [DOI].

Abadias, G., E. Chason, J. Keckes, M. Sebastiani, G. B. Thompson, E. Barthel, G. L. Doll, C. E. Murray, C. H. Stoessel, and L. Martinu (2018). Review Article: Stress in thin films and coatings: Current status, challenges, and prospects. *Journal of Vacuum Science & Technology A: Vacuum, Surfaces, and Films* 36(2):020801. [DOI], [OA].
Barnett, D. and W. Cai (2018). Properties of the Eshelby tensor and existence of the equivalent ellipsoidal inclusion solution. *Journal of the Mechanics and Physics of Solids* 121:71–80. [DOI].

Cai, W. (2018). Matlab codes that implement the calculation of the Eshelby tensors in isotropic and anisotropic elastic media. [URL].

Debotton, G (2005). Transversely isotropic sequentially laminated composites in finite elasticity. *Journal of the Mechanics and Physics of Solids* 53(6):1334–1361. [DOI].

Du, C., Q. Pan, S. Chen, and S. Tian (2020). Effect of rolling process on fatigue performance of stir zone of AA6061-T6 double-side friction stir welding joint. *Fatigue & Fracture of Engineering Materials & Structures* 44(3):748–761. [DOI].

El Omri, A., A. Fennan, F. Sidoroff, and A. Hihi (2000). Elastic-plastic homogenization for layered composites. *European Journal of Mechanics - A/Solids* 19(4):585–601. [DOI], [HAL].

Eshelby, J. (1957). The determination of the elastic field of an ellipsoidal inclusion and related problems. *Proceedings of the Royal Society of London. Series A. Mathematical and Physical Sciences* 241(1226):376–396. [DOI], [HAL].

Franciosi, P. and S. Berbenni (2007). Heterogeneous crystal and poly-crystal plasticity modeling from a transformation field analysis within a regularized Schmid law. *Journal of the Mechanics and Physics of Solids* 55(11):2265–2299. [DOI], [HAL].

Gavazzi, A. C. and D. C. Lagoudas (1990). On the numerical evaluation of Eshelby’s tensor and its application to elastoplastic fibrous composites. *Computational Mechanics* 7(1):13–19. [DOI].

Gemperlová, J. and V. Paidar (1985). Stresses in deformed cubic bicrystals with a tilt (001) grain boundary. *Czechoslovak Journal of Physics* B 35(3):351–354. [DOI], [HAL].

Gemperlová, J., V. Paidar, and F. Kroupa (1989). Compatibility stresses in deformed bicrystals. *Czechoslovak Journal of Physics* B 39(4):427–446. [DOI], [HAL].

Glüge, R. (2016). Effective plastic properties of laminates made of isotropic elastic plastic materials. *Composite Structures* 149:434–443. [DOI], [HAL].

Glüge, R. and J. Kalisch (2014). The effective stiffness and stress concentrations of a multi-layer laminate. *Composite Structures* 111:580–586. [DOI].

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

Hershey, A. V. (1954). The elasticity of an isotropic aggregate of anisotropic cubic crystals. *Journal of Applied Mechanics* 21(3):236–240. [DOI], [HAL].

Hill, R. (1965a). A self-consistent mechanics of composite materials. *Journal of the Mechanics and Physics of Solids* 13(4):213–222. [DOI], [HAL].

Hill, R. (1965b). Continuum micro-mechanics of elastoplastic polycrystals. *Journal of the Mechanics and Physics of Solids* 13(2):89–101. [DOI], [HAL].

Idiart, M. I. (2008). Modeling the macroscopic behavior of two-phase nonlinear composites by infinite-rank laminates. *Journal of the Mechanics and Physics of Solids* 56(8):2599–2617. [DOI].

Kröner, E. (1958). Berechnung der elastischen Konstanten des Vielkristalls aus den Konstanten des Einkristalls. *Zeitschrift für Physik* 151(4):504–518. [DOI].

Lhadi, S., S. Berbenni, N. Gey, T. Richeton, and L. Germain (2018). Micromechanical modeling of the effect of elastic and plastic anisotropies on the mechanical behavior of β-Ti alloys. *International Journal of Plasticity* 109:88–107. [DOI], [HAL].

Lhadi, S., R. r. p. Purushottam raj purohit, T. Richeton, N. Gey, S. Berbenni, O. Perroud, and L. Germain (2020). Elasto-viscoplastic tensile behavior of as-forged Ti-1023 alloy: Experiments and micromechanical modeling. *Materials Science and Engineering: A* 787:139491. [DOI], [OA].

Milton, G. (2002). *The Theory of Composites*. Cambridge University Press. [DOI].

Mukhopadhyay, T., A. Mahata, S. Adhikari, and M. A. Zaeem (2017). Effective mechanical properties of multilayer nano-heterostructures. *Scientific Reports* 7(15818):1–13. [DOI], [OA].

Peralta, P., A. Schober, and C. Laird (1993). Elastic stresses in anisotropic bicrystals. *Materials Science and Engineering: A* 169(1-2):43–51. [DOI].

Purushottam raj purohit, R. r. p., T. Richeton, S. Berbenni, L. Germain, N. Gey, T. Connolley, and O. Castelnau (2021). 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 Materialia* 208:116762 (1–21). [DOI], [HAL].

Quey, R., A. Villani, and C. Maurice (2018). Nearly uniform sampling of crystal orientations.
*Journal of Applied Crystallography* 51(4):1162-1173. [DOI], [HAL].

Rey, C. and A. Zaoui (1980). Slip heterogeneities in deformed aluminum bicrystals. *Acta Metallurgica* 28(6):687-697. [DOI], [HAL].

Richeton, T. and S. Berbenni (2013). Effects of heterogeneous elasticity coupled to plasticity on stresses and lattice rotations in bicrystals: a Field Dislocation Mechanics viewpoint. *European Journal of Mechanics - A/Solids* 37:231-247. [DOI], [HAL].

Richeton, T. and S. Berbenni (2014). From bicrystals to spherical inclusions: A superposition method to derive analytical expressions of stress fields in presence of plastic strain gradients. *International Journal of Solids and Structures* 51(3-4):794-807. [DOI], [OA].

Richeton, T., I. Tiba, S. Berbenni, and O. Bouaziz (2015). Analytical expressions of incompatibility stresses at Σ3⟨111⟩ twin boundaries and consequences on single-slip promotion parallel to twin plane. *Philosophical Magazine* 95:12-31. [DOI], [HAL].

Richeton, T. (2017). Incompatibility stresses and lattice rotations due to grain boundary sliding in heterogeneous anisotropic elasticity. *Crystals* 7 (203):1-14. [DOI], [OA].

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

Stupkiewicz, S. and H. Petryk (2002). Modelling of laminated microstructures in stress-induced martensitic transformations. *Journal of the Mechanics and Physics of Solids* 50(11):2303-2331. [DOI].

Tiba, I., T. Richeton, C. Motz, H. Vehoff, and S. Berbenni (2015). Incompatibility stresses at grain boundaries in Ni bicrystalline micropillars analyzed by an anisotropic model and slip activity. *Acta Materialia* 83:227-238. [DOI], [HAL].

Voigt, W. (1966). *Lehrbuch der Kristallphysik*. Springer. [DOI].

Walpole, L. (1969). On the overall elastic moduli of composite materials. *Journal of the Mechanics and Physics of Solids* 17(4):235-251. [DOI].

Wang, J., Q. Zhou, S. Shao, and A. Misra (2017). Strength and plasticity of nanolaminated materials. *Materials Research Letters* 5:1-19. [DOI], [OA].

Welzel, U., M. Leoni, and E. J. Mittemeijer (2003). The determination of stresses in thin films; modelling elastic grain interaction. *Philosophical Magazine* 83(5):603-630. [DOI], [HAL].

Wong, E. H. and T. B. Lim (2008). A More Comprehensive Solution for Tri-Material Layers Subjected to Thermal Stress. *IEEE Transactions on Components and Packaging Technologies* 31(1):54-64. [DOI], [HAL].

**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the authors—the copyright holder. To view a copy of this license, visit creativecommons.org/licenses/by/4.0.

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

**Authors' contributions** The author confirms sole responsibility for carrying the study, performing the numerical simulations and writing the manuscript.

**Supplementary Material** MATLAB codes to compute the effective behavior and the stress fields in an N-phase laminate, application to an alloy with elongated grains and comparisons with self-consistent models. Freefem++ scripts to check the finite element validity of the laminate solutions. Permalink: swh:1:dir:d7373f7f498776bffeaecof0507a3a637cf39e1e; origin=https://github.com/AniPlas/Laminate;visit=swh:1:snp:07d3ff0677884a21549531f80eb9b316cefc3a97;anchor=swh:1:rev:6bfcf8de1452757d7b04ee17e24a5d3ef5667f97.

**Acknowledgements** This work was supported by the French State through the program "Investment in the future" operated by the National Research Agency (ANR) and referenced by ANR-11-LABX-0008-01 (LabEx DAMAS).

**Ethics approval and consent to participate** Not applicable.

**Consent for publication** Not applicable.

**Competing interests** The author declares that they have no competing interests.

**Journal's Note** JTCAM remains neutral with regard to the content of the publication and institutional affiliations.