Maximum-Entropy Based Estimates of Stress and Strain in Thermoelastic Random Heterogeneous Materials

Abstract

Mean-field methods are a common procedure for characterizing random heterogeneous materials. However, they typically provide only mean stresses and strains, which do not always allow predictions of failure in the phases since exact localization of these stresses and strains requires exact microscopic knowledge of the microstructures involved, which is generally not available. In this work, the maximum entropy method pioneered by Kreher and Pompe (Internal Stresses in Heterogeneous Solids, Physical Research, vol. 9, 1989) is used for estimating one-point probability distributions of local stresses and strains for various classes of materials without requiring microstructural information beyond the volume fractions. This approach yields analytical formulae for mean values and variances of stresses or strains of general heterogeneous linear thermoelastic materials as well as various special cases of this material class. Of these, the formulae for discrete-phase materials and the formulae for polycrystals in terms of their orientation distribution functions are novel. To illustrate the theory, a parametric study based on Al-Al₂O₃ composites is performed. Polycrystalline copper is considered as an additional example. Through comparison with full-field simulations, the method is found to be particularly suited for polycrystals and materials with elastic contrasts of up to 5. We see that, for increasing contrast, the dependence of our estimates on the particular microstructures is increasing, as well.

1 Introduction

Heterogeneous materials, owing to their fabrication process, generally possess random microstructures, allowing for the application of statistical continuum theories to mechanical problems, as described by, e.g., Beran [3]. By using this framework to project the random heterogeneous material properties onto homogeneous effective properties, the “mean field” problem is obtained out of the more complex “full field” problem. This projection is achieved via homogenization methods, which can be numerical in nature, such as those used in FE [9], FFT [17] or NTFA [8] approaches. These numerical methods depend on full-field calculations of representative volumes as an intermediate step, which requires detailed microstructure knowledge. Full-field calculations are generally more numerically expensive than analytical solutions, which offer explicit or implicit solutions in terms of straightforward formulas. Exact (i.e., precise to arbitrary tolerances) analytical homogenizations are known for special microstructures [2], while analytical estimates such as self-consistent estimates, the Mori-Tanaka method [16] and differential schemes [19] suffice in many common applications. For an overview of these and other homogenization methods, see Nemat-Nasser [18] or Torquato [24] and the references therein. These methods generally require only the one-point probability function of material properties, along with some microstructural assumptions.

Solving the homogenized problem yields the mean or effective fields. Many homogenization methods additionally allow calculation of phase mean fields from the overall mean fields by providing phase averages of stress or strain localization tensors; in the two-phase case, they can even be calculated exactly. However, considering only the overall or phase mean fields is insufficient for problems in which local fluctuations determine the result, such as material failure, for which no homogenization theory is known to the authors (as discussed by, e.g., Kröner [15] and Dyskin [6]). For these problems, it is necessary to extract information about local fluctuations from the effective field. While this is again possible by performing numerical simulations of fully known microstructures, as is generally done alongside numerical homogenization methods in two-scale simulations such as FE², the required microstructural data is not easily acquired.

Kreher and Pompe [14] detail an approach to produce analytical approximations of first- and second-order statistical moments of stresses and strains in heterogeneous thermoelastic materials under consideration of effective material values obtained by homogenization methods or experiments. This approach employs the maximum entropy methods commonly seen in statistical thermodynamics and physics (see Jaynes [11], [12]), and is therefore often abbreviated as MEM in the following.

Other approaches have been used to extract information in terms of higher statistical moments of stress and strain distributions. Kreher and Pompe [13] mention an exact analytical solution obtainable by differentiating the effective strain energy, if its dependence on the phase stiffnesses is known. This approach was used most prominently in Ponte Castaneda’s works on second-order homogenization methods [20].

The statistical nature of this work is reminiscent of, but only tangentially related to, Uncertainty Quantification (UQ), the statistical study of models operating on uncertain parameters [5]. In the language of UQ, the method considered here calculates the aleatoric uncertainty of the stress and strain values of a single randomly chosen point of the microstructure, if the effective properties of the continuum, the phase volume fractions, phase properties, and the effective load are known exactly, but the microstructure is entirely uncertain. As far as the authors are aware, there are no UQ works considering that particular problem.

In the present work, the approach of Kreher and Pompe [14] is recapitulated and updated. Some improvements are made concerning notation and mathematical rigorousness. Furthermore, analytical solutions for many material classes are presented, expanding on the depiction of Kreher and Pompe [14]. These range from the isotropic two-phase linear elastic case to the general linear thermoelastic anisotropic case with arbitrary distributions of material properties. Among the cases of intermediate complexity, various combinations of global and local isotropy are discussed. The discussion also touches on polycrystals, yielding a solution in terms of orientation distribution functions (ODFs). Example calculations are compared with FFT full-field simulations. A parametric study serves to illustrate the dependencies of the MEM solution on its input. This paper is intended to offer some didactic value and stimulate interest in maximum entropy methods. The authors hope to provide an introductory text and reference work suited for practical application. For that reason compact overviews are given for the special cases discussed. The explicit formulae depicted in boxes are meant to be applied by engineers who wish to directly estimate the stresses and strains in the materials they are considering.

The manuscript is organized as follows. In Sect. 2, the method of maximal information entropy is briefly recapitulated, followed by a short explanation of its application to continuum mechanics. In Sect. 3, a selection of formulae for various material classes is presented. In Sect. 4, Al-Al₂O₃ composites of various microstructures and polycrystalline copper serve to illustrate the theory. In Sect. 5, the results are summarized.

Notation

A symbolic tensor notation is preferred throughout the text. Scalars are denoted by light-face type characters, e.g., $a,b,\alpha ,\beta ,W$. First-order tensors are denoted by bold-face type lower case Latin characters, e.g., $\boldsymbol{x}$, $\boldsymbol{y}$. Second-order tensors are denoted by bold upper case Latin characters and Greek characters, e.g., $\boldsymbol{A}$, $\boldsymbol{\sigma}$, $\boldsymbol{\varepsilon}$. Fourth-order tensors are denoted by blackboard bold upper case Latin characters, e.g., ${\mathbb{ C}}, {\mathbb{ S}}$. The scalar product between first-, second-, and fourth-order tensors is denoted as $\boldsymbol{a}\cdot \boldsymbol{b}$, $\boldsymbol{A}\cdot \boldsymbol{B}$, and ${\mathbb{ A}}\cdot {\mathbb{ B}}$, respectively. The tensor product is denoted by ⊗. The Rayleigh product is denoted by $\boldsymbol{Q}\star \mathbb{ A}$; it denotes a rotation of the tensor ${\mathbb{ A}}$, if $\boldsymbol{Q}$ is a proper orthogonal tensor ($\boldsymbol{Q}\in \mathit{Orth}^{+}$). The components of any tensor with respect to its orthonormal basis are denoted by $a_{i}, A_{ij}, C_{\mathit{ijkl}}$ respectively for first-, second-, and fourth-order tensors with the corresponding bases $\{\boldsymbol{e}_{i}\}$, $\{\boldsymbol{e}_{i}\otimes \boldsymbol{e}_{j} \}$ and $\{\boldsymbol{e}_{i}\otimes \boldsymbol{e}_{j} \otimes \boldsymbol{e}_{k}\otimes \boldsymbol{e}_{l}\}$. Einstein’s summation convention is applied wherever this index notation is used. The identity on symmetric second-order tensors is denoted by ${\mathbb{ I}}^{\mathrm{{S}}}$ and has the components $I^{\mathrm{{S}}}_{\mathit{ijkl}} = (\delta _{ik}\delta _{jl} + \delta _{il} \delta _{jk})/2$. It can be decomposed into the projectors ${\mathbb{ P}}_{1} = \boldsymbol{1}\otimes \boldsymbol{1}/ 3$ and ${\mathbb{ P}}_{2} = {\mathbb{ I}}^{\mathrm{{S}}} - {\mathbb{ P}}_{1}$. A stiffness tensor ℂ is a positive definite fourth-order tensor with minor and major symmetries, i.e. $C_{\mathit{ijkl}} = C_{\mathit{jikl}} = C_{\mathit{klij}}$. The linear map of a second-order tensor over a fourth-order tensor is denoted by ${\mathbb{ C}}[\boldsymbol{\varepsilon}]$. A bilinear form $\boldsymbol{\varepsilon}\cdot {\mathbb{ C}}[ \boldsymbol{\varepsilon}]$ is computed as $\varepsilon _{ij} C_{\mathit{ijkl}} \varepsilon _{kl}$.The inverse of a fourth-order tensor with respect to ${\mathbb{ I}}^{\mathrm{{S}}}$ is denoted by ${{\mathbb{ C}}}^{\mathrm{-1}}$. The expectation value of a quantity $\psi$ is denoted as $\langle \psi \rangle$. Effective quantities are denoted by a bar, e.g., $\bar{\psi }$.

2 Maximum Entropy Methods in Continuum Mechanics

2.1 General Principle

The problem of localization in continuum mechanics of random heterogeneous materials can be formulated as follows: Given a material of known components but unknown exact structure, as well as known effective material behavior including effective stresses and strains, what is known about the local stresses and strains? In the maximum entropy approach, the unknown properties are modeled as random properties. As the material’s components are generally known, it is straightforward to find the probability function of material properties, but difficult to derive from those the probability distributions of local stresses and strains, even though the mean stresses and strains are known.

Consider a fully unknown probability distribution $p(\boldsymbol{x})$, such as the distribution of strains $p(\boldsymbol{\varepsilon})$, of which some macroscopic values such as means or variances are known. Generally, there are an infinite number of possible $p$. One can calculate the information entropy or information content of a probability distribution according to Jaynes [11] as follows:

$$ S = - \int p(\boldsymbol{x}) {\ln \left (p( \boldsymbol{x})\right )}\!\mathrm{d} \boldsymbol{x}. $$

(1)

The principle of maximal information entropy distinguishes those $p$ which contain a maximum of information. A maximum of random information corresponds to a minimum of determined information. Any such determined information could be formulated as further assumptions about the probability distribution. Therefore, the state of maximal information entropy corresponds to a state of minimal unfounded assumptions. The approach may be formally stated as:

$$\underset{p}{max}S,F_{\beta }(p)=0,$$

(2)

where $F_{\beta }$ are functionals representing constraints any $p$ must satisfy.

A solution of maximal entropy is generally an approximation. It does not correspond to the real solution because not all real constraints were considered. For example, when considering heterogeneous materials, one does not model every production parameter which might impact the final structure.

2.2 Application to Linear Thermoelastic Materials

Linear thermoelastic materials can be described by a fourth-order elasticity tensor ℂ, the second-order tensor of thermal expansion coefficients $\boldsymbol{\alpha}$ and the scalar temperature difference $\Delta \theta$ with respect to a reference temperature. In this framework, the constitutive law takes the form:

$$ \boldsymbol{\sigma}= {\mathbb{ C}}[ \boldsymbol{\varepsilon}- \boldsymbol{\alpha}\Delta \theta ]. $$

(3)

By assuming the temperature difference to be homogeneous, it is possible to consider $\boldsymbol{\varepsilon}^{\uptheta }:= \boldsymbol{\alpha}\Delta \theta$ as the relevant material constant instead. This “stress-free strain” can also be used to model various other conditions of internal stress, such as those caused by phase transitions in crystals.

To allow for easier calculations in the following steps, it is expedient to additively divide the strain $\boldsymbol{\varepsilon}$ into two parts, one associated with an external stress and one with the stress-free strain. By the principle of superposition, the thus defined load cases $\mathrm{{I}}$ and $\mathrm{{I\!I}}$ can model the original combined load, as discussed by Kreher and Pompe [14]. The constitutive laws are as follows:

$${\mathit{\sigma }}^{\mathrm{I}}=\mathbb{C}[{\mathit{\epsilon }}^{\mathrm{I}}],$$

(4)

$${\mathit{\sigma }}^{\mathrm{II}}=\mathbb{C}[{\mathit{\epsilon }}^{\mathrm{II}}-{\mathit{\epsilon }}^{\mathrm{\Theta }}],〈{\mathit{\sigma }}^{\mathrm{II}}〉=0.$$

(5)

By the stochastic Hill-Mandel condition as defined by Kreher and Pompe [14] (cf. the more commonly used ergodic Hill-Mandel condition [10]), effective homogenized stresses and strains can be defined. For the load cases $\mathrm{{I}}$ and $\mathrm{{I\!I}}$, they are:

$$〈{\mathit{\sigma }}^{\mathrm{I}}〉={\overline{\mathit{\sigma }}}^{\mathrm{I}}=\overline{\mathbb{C}}[{\overline{\mathit{\epsilon }}}^{\mathrm{I}}],{\overline{\mathit{\epsilon }}}^{\mathrm{I}}=〈{\mathit{\epsilon }}^{\mathrm{I}}〉,$$

(6)

$$〈{\mathit{\sigma }}^{\mathrm{II}}〉=\mathbf{0}={\overline{\mathit{\sigma }}}^{\mathrm{II}},{\overline{\mathit{\epsilon }}}^{\mathrm{II}}=〈{\mathit{\epsilon }}^{\mathrm{II}}〉.$$

(7)

The effective elasticity tensor $\bar{{\mathbb{ C}}}$ cannot generally be calculated from the phase elasticity tensors and must therefore be prescribed independently. There is another similarly independent quantity, the effective strain energy density. The strain energy density is defined as

$$ w = \int _{\boldsymbol{\varepsilon}^{\uptheta }}^{ \boldsymbol{\varepsilon}} \boldsymbol{\sigma}\cdot \mathop{}\!\mathrm{d}\boldsymbol{\varepsilon}= \frac{1}{2} \boldsymbol{\sigma}\cdot (\boldsymbol{\varepsilon}- \boldsymbol{\varepsilon}^{\uptheta }). $$

(8)

In load case $\mathrm{{I}}$, the effective strain energy density is given by the stochastic Hill-Mandel condition as $\bar{w}^{\mathrm{{I}}}= \langle \boldsymbol{\sigma} \rangle \cdot \langle \boldsymbol{\varepsilon}\rangle$. However, $\bar{w}^{\mathrm{{I\!I}}}$ is more complicated:

$$ \bar{w}^{\mathrm{{I\!I}}}= \frac{1}{2}\langle \boldsymbol{\sigma} \cdot (\boldsymbol{\varepsilon}^{\mathrm{{I\!I}}}- \boldsymbol{\varepsilon}^{\uptheta })\rangle = \frac{1}{2} \bar{\boldsymbol{\sigma}}^{\mathrm{{I\!I}}}\cdot \bar{\boldsymbol{\varepsilon}}^{\mathrm{{I\!I}}}- \frac{1}{2} \langle \boldsymbol{\sigma}^{\mathrm{{I\!I}}}\cdot \boldsymbol{\varepsilon}^{\uptheta }\rangle = \mbox{$ 0 $} - \frac{1}{2}\langle \boldsymbol{\sigma}^{\mathrm{{I\!I}}}\cdot \boldsymbol{\varepsilon}^{\uptheta }\rangle . $$

(9)

This quantity cannot be described macroscopically by effective stresses, strains, or $\bar{{\mathbb{ C}}}$ and is therefore prescribed independently.

The maximum entropy method for linear thermoelastic localization problems can now be stated as follows: With a given random distribution of material properties $p_{1}^{\mathrm{{C}}}(\boldsymbol{\varepsilon}^{\uptheta }, { \mathbb{ C}})$, given effective material properties $\bar{{\mathbb{ C}}}$ and $\bar{u}^{\mathrm{{I\!I}}}$, and the effective homogenized strains $\bar{\boldsymbol{\varepsilon}}^{\mathrm{{I}}}$ and $\bar{\boldsymbol{\varepsilon}}^{\mathrm{{I\!I}}}$, calculate the joint probability distribution $p_{1}({\mathbb{ C}}, \boldsymbol{\varepsilon}^{\uptheta }, \boldsymbol{\varepsilon}^{\mathrm{{I}}}, \boldsymbol{\varepsilon}^{\mathrm{{I\!I}}})$ by maximizing the information entropy of $p_{1}$ while honoring various constraints derived from mechanical laws. These constraints are the prescribed effective values, yielding five equations:

$$\begin{aligned} \bar{\boldsymbol{\varepsilon}}^{\mathrm{{I}}}&= \langle \boldsymbol{\varepsilon}^{\mathrm{{I}}}\rangle , \end{aligned}$$

(10)

$$\begin{aligned} \bar{\boldsymbol{\varepsilon}}^{\mathrm{{I\!I}}}&= \langle \boldsymbol{\varepsilon}^{\mathrm{{I\!I}}}\rangle , \end{aligned}$$

(11)

$$\begin{aligned} \bar{{\mathbb{ C}}}[\bar{\boldsymbol{\varepsilon}}^{\mathrm{{I}}}] &= \langle {\mathbb{ C}}[\boldsymbol{\varepsilon}^{\mathrm{{I}}}] \rangle , \end{aligned}$$

(12)

$$\begin{aligned} \langle {\mathbb{ C}}[\boldsymbol{\varepsilon}^{\mathrm{{I\!I}}}- \boldsymbol{\varepsilon}^{\uptheta }]\rangle &= \bar{{\mathbb{ C}}}[\bar{\boldsymbol{\varepsilon}}^{\mathrm{{I \!I}}}-\bar{\boldsymbol{\varepsilon}}^{\uptheta }] = \mbox{$ 0 $}, \end{aligned}$$

(13)

$$\begin{aligned} \langle \boldsymbol{\varepsilon}^{\uptheta }\cdot {\mathbb{ C}}[ \boldsymbol{\varepsilon}^{\mathrm{{I\!I}}}- \boldsymbol{\varepsilon}^{\uptheta }]\rangle &= -2\bar{w}^{\mathrm{{I \!I}}}, \end{aligned}$$

(14)

the stochastic Hill-Mandel condition, which states that $\langle \boldsymbol{\varepsilon}\cdot \boldsymbol{\sigma}\rangle = \langle \boldsymbol{\varepsilon}\rangle \cdot \langle \boldsymbol{\sigma}\rangle$ for any compatible $\boldsymbol{\varepsilon}$ and divergence-free $\boldsymbol{\sigma}$, yielding four equations:

$$\begin{aligned} \langle \boldsymbol{\varepsilon}^{\mathrm{{I}}}\cdot {\mathbb{ C}}[ \boldsymbol{\varepsilon}^{\mathrm{{I}}}]\rangle &= \bar{\boldsymbol{\varepsilon}}^{\mathrm{{I}}}\cdot \bar{{\mathbb{ C}}}[\bar{\boldsymbol{\varepsilon}}^{\mathrm{{I}}}], \end{aligned}$$

(15)

$$\begin{aligned} \langle \boldsymbol{\varepsilon}^{\mathrm{{I\!I}}}\cdot { \mathbb{ C}}[\boldsymbol{\varepsilon}^{\mathrm{{I}}}]\rangle &= \bar{\boldsymbol{\varepsilon}}^{\mathrm{{I\!I}}}\cdot \bar{{\mathbb{ C}}}[\bar{\boldsymbol{\varepsilon}}^{\mathrm{{I}}}], \end{aligned}$$

(16)

$$\begin{aligned} \langle \boldsymbol{\varepsilon}^{\mathrm{{I}}}\cdot {\mathbb{ C}}[ \boldsymbol{\varepsilon}^{\mathrm{{I\!I}}}- \boldsymbol{\varepsilon}^{\uptheta }]\rangle &= \bar{\boldsymbol{\varepsilon}}^{\mathrm{{I}}}\cdot \bar{{\mathbb{ C}}}[\bar{\boldsymbol{\varepsilon}}^{\mathrm{{I \!I}}}-\bar{\boldsymbol{\varepsilon}}^{\uptheta }] = 0, \end{aligned}$$

(17)

$$\begin{aligned} \langle \boldsymbol{\varepsilon}^{\mathrm{{I\!I}}}\cdot { \mathbb{ C}}[\boldsymbol{\varepsilon}^{\mathrm{{I\!I}}}- \boldsymbol{\varepsilon}^{\uptheta }]\rangle &= \bar{\boldsymbol{\varepsilon}}^{\mathrm{{I\!I}}}\cdot \bar{{\mathbb{ C}}}[\bar{\boldsymbol{\varepsilon}}^{\mathrm{{I \!I}}}-\bar{\boldsymbol{\varepsilon}}^{\uptheta }] = 0, \end{aligned}$$

(18)

and finally, the prescribed distribution of material properties, which can be calculated from $p_{1}$ by marginalization:

$$ \iint p_{1}(\boldsymbol{\varepsilon}^{\mathrm{{I}}}, \boldsymbol{\varepsilon}^{\mathrm{{I\!I}}}, {\mathbb{ C}}, \boldsymbol{\varepsilon}^{\uptheta })\mathop{}\!\mathrm{d} \boldsymbol{\varepsilon}^{\mathrm{{I}}}\mathop{}\!\mathrm{d} \boldsymbol{\varepsilon}^{\mathrm{{I\!I}}}= p_{1}^{\mathrm{{C}}}({ \mathbb{ C}}, \boldsymbol{\varepsilon}^{\uptheta }). $$

(19)

These constraints are implemented into the maximization of the information entropy by the method of Lagrange multipliers (see, e.g., Bertsekas [4]). A detailed explanation of the calculation of the multipliers used here has been supplied by Kreher and Pompe [14]. As a technical improvement on that account, the multiplier for the final constraint is formulated as a function-valued constraint from the space of square-integrable functions $\mathcal{L}^{2}$. Multiplying the constraint with its Lagrange multiplier $\lambda _{10}({\mathbb{ C}}, \boldsymbol{\varepsilon}^{\uptheta })$ via the dot product in $\mathcal{L}^{2}$ yields

$$\begin{aligned} &\iint \lambda _{10}({\mathbb{ C}}, \boldsymbol{\varepsilon}^{\uptheta }) p_{1}^{\mathrm{{C}}}({\mathbb{ C}}, \boldsymbol{\varepsilon}^{\uptheta }) \mathop{}\!\mathrm{d}{ \mathbb{ C}}\mathop{}\!\mathrm{d}\boldsymbol{\varepsilon}^{\uptheta }= \\ &\qquad \iiiint \lambda _{10}({\mathbb{ C}}, \boldsymbol{\varepsilon}^{\uptheta }) p_{1}( \boldsymbol{\varepsilon}^{\mathrm{{I}}}, \boldsymbol{\varepsilon}^{\mathrm{{I\!I}}}, {\mathbb{ C}}, \boldsymbol{\varepsilon}^{\uptheta })\mathop{}\!\mathrm{d} \boldsymbol{\varepsilon}^{\mathrm{{I}}}\mathop{}\!\mathrm{d} \boldsymbol{\varepsilon}^{\mathrm{{I\!I}}}\mathop{}\!\mathrm{d}{ \mathbb{ C}}\mathop{}\!\mathrm{d}\boldsymbol{\varepsilon}^{\uptheta }, \end{aligned}$$

(20)

after which further calculations follow the same path as those in Kreher and Pompe [14]. While this approach arguably stands on firmer mathematical foundations, it imparts the restriction that $p_{1}^{\mathrm{{C}}} \in L \subset \mathcal{L}^{2}$, where $L$ is finite-dimensional. A formulation of Lagrange multipliers which accounts for infinite-dimensional constraints would resolve this restriction, but no such formulation is known to the authors. In practical application, this restriction is generally no hindrance, as all material property distributions can be approximated through finite-dimensional functions.

After the Lagrange multipliers have been calculated, the method arrives at an analytical solution for the probability distribution $p_{1}$. The Lagrange multipliers are then substituted with new variables to allow for a more compact depiction, which will be shown in the following section along with several special cases.

3 Solutions for Various Material Classes

This section is organized as a collection of formulae intended for practical application. As such, we begin with a short definition of various descriptors:

• Discrete phases. For materials with discrete phases, the material property distribution function can be described as a finite sum of Dirac delta distributions $p_{1}^{\mathrm{{C}}}({\mathbb{ C}}, \boldsymbol{\varepsilon}^{\uptheta }) = \sum _{\alpha = 1}^{N} c_{\alpha }\delta ({\mathbb{ C}}- {\mathbb{ C}}_{\alpha }) \delta ( \boldsymbol{\varepsilon}^{\uptheta }- \boldsymbol{\varepsilon}^{\uptheta }_{\alpha })$, where $N$ is the number of phases and $c_{\alpha }$ their respective incidence probability, generally modeled through the volume fraction. This is not to be confused with the common term “piecewise constant material” which denotes a material whose spatial distribution of material properties ${\mathbb{ C}}(\boldsymbol{x})$, $\boldsymbol{\varepsilon}^{\uptheta }(\boldsymbol{x})$ is piecewise constant. Note that in this framework every orientation of an anisotropic material constitutes a formally different phase. As such, some common piecewise constant materials, such as isotropically distributed uniform polycrystals, have an infinite number of phases and are thus non-discrete. Materials with discrete phases are however always piecewise constant, excepting unrealistic mathematical curiosities such as a material whose phase changes from material point to material point.

• Global and local isotropy. Locally isotropic materials consist of isotropic phases. Globally isotropic materials are characterized by isotropic effective properties, as well as isotropic first-order bounds on these properties. This concept of global isotropy is a restriction on the one-point distribution of material properties, not higher-point correlation functions; furthermore, it does not take into account tensors of higher order than four. Therefore, it is a strictly weaker assumption than statistical isotropy as defined by e.g. Torquato [24]. As far as the formulae of the maximum entropy method are concerned, both concepts are equivalent and will be used interchangeably.

• Polycrystals. The material properties ℂ and $\boldsymbol{\varepsilon}^{\uptheta }$ of a multi-phase or non-uniform polycrystal can be described as a finite sum of reference properties ${\mathbb{ C}}_{\alpha }$ and $\boldsymbol{\varepsilon}^{\uptheta }_{\alpha }$ transformed by their orientations $\boldsymbol{Q}_{\alpha }\in \mathit{Orth}^{+}$. After applying this substitution, the material property distribution function is no longer a function $p_{1}^{\mathrm{{C}}}({\mathbb{ C}}, \boldsymbol{\varepsilon}^{\uptheta })$, but rather $p_{1}^{\mathrm{{C}}}(\boldsymbol{Q}) = \sum _{\alpha = 1}^{N} c_{\alpha }p(\boldsymbol{Q}_{\alpha })$. In this formulation, $p(\boldsymbol{Q}_{\alpha })$ is the orientation distribution function (ODF) of phase $\alpha$.

Using these definitions as well as the distinction between linear thermoelastic materials and the special case of linear elastic materials, where $\boldsymbol{\varepsilon}^{\uptheta }$ is constant across phases, yields the material classification depicted in Fig. 1. In Table 1, this classification is related to a table of contents of this section. The sections marked with asterisks were first derived by Kreher and Pompe [14], though they did not formulate all of these solutions explicitly. Note that the cases marked with asterisks result in singularities when using the more general forms, while many of the non-marked solutions were only derived for convenience of use.

Fig. 1

Set diagram of material classes discussed in this paper

Table 1 Table of contents for Sect. 3 by material class. Asterisks denote formulae adapted from Kreher and Pompe [14]

The following solutions are all formulated in terms of strains. The stress fields $\boldsymbol{\sigma}^{\mathrm{{I}}}$ and $\boldsymbol{\sigma}^{\mathrm{{I\!I}}}$ can be obtained through transformation via the constitutive laws of equations (4) and (5):

$$\begin{aligned} \begin{pmatrix} \boldsymbol{\sigma}^{\mathrm{{I}}} \\ \boldsymbol{\sigma}^{\mathrm{{I\!I}}} \end{pmatrix} &= \begin{pmatrix} {\mathbb{ C}}& 0 \\ 0 & {\mathbb{ C}} \end{pmatrix} \begin{bmatrix} \boldsymbol{\varepsilon}^{\mathrm{{I}}} \\ \boldsymbol{\varepsilon}^{\mathrm{{I\!I}}}- \boldsymbol{\varepsilon}^{\uptheta }\end{bmatrix}. \end{aligned}$$

(21)

The phase stress distribution can be obtained through the same transformation applied to the phase strain distribution, yielding simple forms for the mean stress and stress covariance matrix:

$$\begin{aligned} \begin{pmatrix} \boldsymbol{\gamma}^{\mathrm{{I}}}_{\sigma }\\ \boldsymbol{\gamma}^{\mathrm{{I\!I}}}_{\sigma }\end{pmatrix} &= \begin{pmatrix} {\mathbb{ C}}[\boldsymbol{\gamma}^{\mathrm{{I}}}] \\ {\mathbb{ C}}[\boldsymbol{\gamma}^{\mathrm{{I}}}- \boldsymbol{\varepsilon}^{\uptheta }] \end{pmatrix}, \end{aligned}$$

(22)

$$\begin{aligned} {\mathbb{ K}}_{\sigma }&= \begin{pmatrix} {\mathbb{ C}}& 0 \\ 0 & {\mathbb{ C}} \end{pmatrix} {\mathbb{ K}} \begin{pmatrix} {\mathbb{ C}}& 0 \\ 0 & {\mathbb{ C}} \end{pmatrix}. \end{aligned}$$

(23)

Note that the symbol ${\mathbb{ K}}$ here denotes a 2x2 array of fourth-order tensors.

3.1 Abbreviations

As the following equations may become very large, it is expedient to define a few general abbreviations for particularly unwieldy terms:

$$ w_{0} = \langle \boldsymbol{\varepsilon}^{\uptheta }\cdot { \mathbb{ C}}[\boldsymbol{\varepsilon}^{\uptheta }]\rangle - \boldsymbol{\varepsilon}_{\mathrm{{s}}} \cdot {\mathbb{ C}}_{+} \boldsymbol{\varepsilon}_{\mathrm{{s}}} - ( \boldsymbol{\varepsilon}_{\mathrm{{s}}} - \boldsymbol{\varepsilon}_{\mathrm{{e}}}) \cdot ({\mathbb{ C}}_{-}^{-1} - {\mathbb{ C}}_{+}^{-1})^{-1}[\boldsymbol{\varepsilon}_{\mathrm{{s}}} - \boldsymbol{\varepsilon}_{\mathrm{{e}}}], $$

(24)

$${\mathbb{C}}_{+}=〈\mathbb{C}〉,{\mathbb{C}}_{-}={〈{\mathbb{C}}^{-1}〉}^{-1},$$

(25)

$${\mathit{\epsilon }}_{\mathrm{e}}=〈{\mathit{\epsilon }}^{\mathrm{\Theta }}〉,{\mathit{\epsilon }}_{\mathrm{s}}={〈\mathbb{C}〉}^{-1}[〈\mathbb{C}[{\mathit{\epsilon }}^{\mathrm{\Theta }}]〉].$$

(26)

3.1.1 Global Isotropy

In the global isotropic case, global material parameters can be simplified using the fourth-order projectors:

$$\begin{aligned} \bar{{\mathbb{ C}}}&= 3 \bar{K}{\mathbb{ P}}_{1} + 2 \bar{G}{ \mathbb{ P}}_{2}, \end{aligned}$$

(27)

$${\mathbb{C}}_{+}=3K_{+}{\mathbb{P}}_1+2G_{+}{\mathbb{P}}_2,3K_{+}={\mathbb{P}}_1\cdot 〈\mathbb{C}〉,2G_{+}=\frac{1}{\sqrt{5}}{\mathbb{P}}_2\cdot 〈\mathbb{C}〉,$$

(28)

$${\mathbb{C}}_{-}=3K_{-}{\mathbb{P}}_1+2G_{-}{\mathbb{P}}_2,3K_{-}={\mathbb{P}}_1\cdot {〈{\mathbb{C}}^{-1}〉}^{-1},2G_{-}=\frac{1}{\sqrt{5}}{\mathbb{P}}_2\cdot {〈{\mathbb{C}}^{-1}〉}^{-1}.$$

(29)

Furthermore, the stress-free strains are spherical:

$$\begin{aligned} \bar{\boldsymbol{\varepsilon}}^{\mathrm{{I\!I}}}&= \bar{\varepsilon }^{\mathrm{{I\!I}}}\boldsymbol{1}, \end{aligned}$$

(30)

$${\mathit{\epsilon }}_{\mathrm{e}}={\epsilon }_e\mathbf{1},{\epsilon }_e=\frac{1}{3}〈{\mathit{\epsilon }}^{\mathrm{\Theta }}〉\cdot \mathbf{1},$$

(31)

$${\mathit{\epsilon }}_{\mathrm{s}}={\epsilon }_s\mathbf{1},{\epsilon }_s=\frac{1}{9K_{+}}〈\mathbb{C}[{\mathit{\epsilon }}^{\mathrm{\Theta }}]〉\cdot \mathbf{1}.$$

(32)

3.1.2 Local Isotropy

Local isotropy implies global isotropy of the various expectation values. Thus it follows that:

$${\mathbb{C}}_{+}=3K_{+}{\mathbb{P}}^1+2G_{+}{\mathbb{P}}^2,K_{+}=〈K〉,G_{+}=〈G〉,$$

(33)

$${\mathbb{C}}_{-}=3K_{-}{\mathbb{P}}^1+2G_{-}{\mathbb{P}}^2,K_{-}={〈K^{-1}〉}^{-1},G_{-}={〈G^{-1}〉}^{-1}.$$

(34)

Again, the stress-free strains are spherical:

$${\mathit{\epsilon }}_{\mathrm{e}}={\epsilon }_e\mathbf{1},{\epsilon }_e=〈{\epsilon }^{\mathrm{\Theta }}〉,$$

(35)

$${\mathit{\epsilon }}_{\mathrm{s}}={\epsilon }_e\mathbf{1},{\epsilon }_s=\frac{〈K{\epsilon }^{\mathrm{\Theta }}〉}{K_{+}}.$$

(36)

It should be noted that the isotropy does not extend to effective elasticity tensors and strains in this case.

3.2 General Solution for Linear Thermoelastic Materials

The general solution is depicted in box 1.

Box 1

MEM solution for general linear thermoelastic materials

3.3 General Solution for Linear Elastic Materials

Linear elastic materials have constant $\boldsymbol{\varepsilon}^{\uptheta }$ and are thus particularly simple. As $\boldsymbol{\varepsilon}^{\uptheta }$ shows no statistical variation, load case $\mathrm{{I\!I}}$ can be completely disregarded in the maximum entropy method. The solution follows as depicted in box 2.

Box 2

MEM solution for general linear elastic materials

3.4 Linear Thermoelastic Materials with Discrete Phases

Materials with discrete phases have material probability distributions $p_{1}^{\mathrm{{C}}}$ that can be described using finite sums of Dirac delta distributions. This means that expectation values containing only material parameters can be expressed as sums:

$$ \langle f({\mathbb{ C}}, \boldsymbol{\varepsilon}^{\uptheta }) \rangle = \sum _{\alpha =0}^{N} c_{\alpha }f({\mathbb{ C}}_{\alpha }, \boldsymbol{\varepsilon}^{\uptheta }_{\alpha }), $$

(37)

where $N$ is the number of phases and $c_{\alpha }$ the phase incidence probability, which is equal to the phase volume fraction in the ergodic case.

This simplification impacts the abbreviations listed in Sect. 3.1 and the form of $p_{1}$. The solution takes the form depicted in box 3.

Box 3

MEM solution for linear thermoelastic materials with discrete phases

3.4.1 Local Anisotropy, Global Isotropy

Applying the discrete simplification to the abbreviations from Sect. 3.1.1 yields globally isotropic discrete abbreviations. Using these, the solution takes the form depicted in box 4.

Box 4

MEM solution for a locally anisotropic, globally isotropic material with discrete phases

3.4.2 Local Isotropy, Global Anisotropy

Applying the discrete simplification to the abbreviation from Sect. 3.1.1 yields globally isotropic discrete abbreviations. Substituting these into the solution of Sect. 3.4 yields no significant simplifications because $\bar{{\mathbb{ C}}}$ is still anisotropic.

However, $\mathrm{det }( {\mathbb{ K}}_{\alpha })$ can be simplified similarly to Sect. 3.6.2, leading to:

$$\begin{aligned} {\mathbb{ K}}_{\alpha }&= \begin{pmatrix} d^{\mathrm{{I}}}{\mathbb{ C}}^{-1}_{\alpha }& k {\mathbb{ C}}^{-1}_{\alpha }\\ k {\mathbb{ C}}^{-1}_{\alpha }& d^{\mathrm{{I\!I}}}{\mathbb{ C}}^{-1}_{\alpha }\end{pmatrix}, \end{aligned}$$

(38)

$$\begin{aligned} \mathrm{det }( {\mathbb{ K}}_{\alpha }) &= (d^{\mathrm{{I}}}d^{\mathrm{{I\!I}}}- k^{2})^{6} {\mathrm{det }}( {\mathbb{ C}}^{-1}_{\alpha })^{2} \\ &= (d^{\mathrm{{I}}}d^{\mathrm{{I\!I}}}- k^{2})^{6} \left ( \frac{(1+\nu _{i})^{3} (1-3\nu ^{2}-2\nu ^{3})}{E^{6}}\right )^{2}. \end{aligned}$$

(39)

3.4.3 Local and Global Isotropy

Combining local and global isotropy yields the formulae in box 5, which are similar to the globally isotropic but locally anisotropic case with an additional simplification of $\boldsymbol{\gamma}$.

Box 5

MEM solution for globally and locally isotropic materials with discrete phases

3.5 Two-Phase Linear Thermoelastic Materials

In the case of only two phases, the phasewise mean strains can be calculated analytically from the effective elasticity tensor $\bar{{\mathbb{ C}}}$. As a result, $\bar{w}^{\mathrm{{I\!I}}}$ and $\bar{\boldsymbol{\varepsilon}}^{\mathrm{{I\!I}}}$ can be computed analytically as well, as demonstrated by Rosen and Hashin [21]:

$$\begin{aligned} \bar{\boldsymbol{\varepsilon}}^{\mathrm{{I\!I}}}&= \langle \boldsymbol{\varepsilon}^{\uptheta }\rangle + (\langle { \mathbb{ C}}^{-1}\rangle - \bar{{\mathbb{ C}}}^{-1})({\mathbb{ C}}_{1}^{-1} - {\mathbb{ C}}_{2}^{-1})^{-1} [\boldsymbol{\varepsilon}^{\uptheta }_{2} - \boldsymbol{\varepsilon}^{\uptheta }_{1}], \end{aligned}$$

(40)

$$\begin{aligned} \bar{w}^{\mathrm{{I\!I}}}&= \frac{1}{2} (\boldsymbol{\varepsilon}^{\uptheta }_{2} - \boldsymbol{\varepsilon}^{\uptheta }_{1})({ \mathbb{ C}}_{1}^{-1} - {\mathbb{ C}}_{2}^{-1})^{-1}[ \bar{\boldsymbol{\varepsilon}}^{\mathrm{{I\!I}}}- \langle \boldsymbol{\varepsilon}^{\uptheta }\rangle ] \\ &= \frac{1}{2} (\boldsymbol{\varepsilon}^{\uptheta }_{2} - \boldsymbol{\varepsilon}^{\uptheta }_{1})({\mathbb{ C}}_{1}^{-1} - {\mathbb{ C}}_{2}^{-1})^{-1}(\langle {\mathbb{ C}}^{-1}\rangle - \bar{{\mathbb{ C}}}^{-1})({\mathbb{ C}}_{1}^{-1} - {\mathbb{ C}}_{2}^{-1})^{-1}[ \boldsymbol{\varepsilon}^{\uptheta }_{2} - \boldsymbol{\varepsilon}^{\uptheta }_{1}]. \end{aligned}$$

(41)

This causes a major simplification of the equations, which are depicted in box 6.

Box 6

MEM solution for two-phase linear thermoelastic materials

3.5.1 Local Anisotropy, Global Isotropy

Using the abbreviations from Sect. 3.1.1 leads to the solution in box 7.

Box 7

MEM solution for locally anisotropic, globally isotropic two-phase linear thermoelastic materials

3.5.2 Local Isotropy, Global Anisotropy

This case is most easily stated using the two-phase anisotropic case in box 6, where the abbreviations from Sect. 3.1.1 should be used. Additionally, $\mathrm{det }( {\mathbb{ K}}_{\alpha })$ can be simplified as described in Sect. 3.4.2, leading to:

$$\begin{aligned} {\mathrm{det }}( {\mathbb{ K}}_{\alpha }) &= (d^{\mathrm{{I}}}d^{\mathrm{{I\!I}}}- k^{2})^{6} \left (\frac{(1+\nu _{i})^{3} (1-3\nu ^{2}-2\nu ^{3})}{E^{6}}\right )^{2}. \end{aligned}$$

(42)

3.5.3 Local and Global Isotropy

Combining global and local isotropies yields a depiction similar to the one for global isotropies with a further simplified formula for $\boldsymbol{\gamma}$. The solution is shown in box 8.

Box 8

MEM solution for locally and globally isotropic two-phase thermoelastic materials

3.6 Linear Elastic Materials with Discrete Phases

Similarly to the general linear elastic case in Sect. 3.3, load case $\mathrm{{I\!I}}$ can be completely disregarded. Contrary to linear thermoelastic materials, there are no special formulae for two-phase linear elastic materials. The linear elastic solution for materials with discrete phases is shown in box 9.

Box 9

MEM solution for linear elastic materials with discrete phases

3.6.1 Local Anisotropy, Global Isotropy

The solution can be obtained analogously to the linear thermoelastic equivalent, leading to box 10.

Box 10

MEM solution for locally anisotropic, globally isotropic linear elastic materials with discrete phases

3.6.2 Local Isotropy, Global Anisotropy

If the local elasticity tensor ℂ is isotropic, so is the local covariance matrix ${\mathbb{ K}}_{i}$. With Poisson’s ratio $\nu$ and Young’s modulus $E$ follows:

$$ \begin{aligned} {\mathbb{ K}}_{\alpha }&= \frac{1}{2\lambda _{Q\mathrm{{I}}}} { \mathbb{ C}}^{-1}, \\ \mathrm{det }( {\mathbb{ K}}_{\alpha }) &= \frac{1}{(2\lambda _{Q\mathrm{{I}}})^{6}} {\mathrm{det }}( {\mathbb{ C}}^{-1} ) \\ & = \frac{1}{(2\lambda _{Q\mathrm{{I}}}E_{\alpha })^{6}} \begin{vmatrix} 1 & -\nu _{\alpha }& -\nu _{\alpha }& 0 & 0 & 0 \\ -\nu _{\alpha }& 1 & -\nu _{\alpha }& 0 & 0 & 0 \\ -\nu _{\alpha }& -\nu _{\alpha }& 1 & 0 & 0 & 0 \\ 0 & 0 & 0 & 1 + \nu _{\alpha }& 0 & 0 \\ 0 & 0 & 0 & 0 & 1 + \nu _{\alpha }& 0 \\ 0 & 0 & 0 & 0 & 0 & 1 + \nu _{\alpha }\end{vmatrix} \\ &= \frac{(1+\nu _{\alpha })^{3} (1-3\nu ^{2}-2\nu ^{3})}{(2\lambda _{Q\mathrm{{I}}}E_{\alpha })^{6}}. \end{aligned} $$

(43)

The rest of the equations are most easily stated using the more general solution in box 9, where the abbreviations from Sect. 3.1.2 should be used.

3.6.3 Local and Global Isotropy

For fully isotropic linear elastic materials, the solution is given in box 11.

Box 11

MEM solution for globally and locally isotropic linear elastic materials with discrete phases

3.7 Linear Thermoelastic Polycrystals

3.7.1 Globally Anisotropic Single-Phase Polycrystals

In single-phase polycrystals, the material properties ℂ and $\boldsymbol{\varepsilon}^{\uptheta }$ can be described by reference material properties ${\mathbb{ C}}_{0}$ and $\boldsymbol{\varepsilon}^{\uptheta }_{0}$ and the orientation tensor of second order $\boldsymbol{Q}\in \mathit{Orth}^{+}$:

$$\begin{aligned} {\mathbb{ C}}&= \boldsymbol{Q}\star {\mathbb{ C}}_{0}, \end{aligned}$$

(44)

$$\begin{aligned} \boldsymbol{\varepsilon}^{\uptheta }&= \boldsymbol{Q} \star \boldsymbol{\varepsilon}^{\uptheta }_{0}, \end{aligned}$$

(45)

and in component notation:

$$\begin{aligned} C_{\mathit{ijkl}} &= Q_{im} Q_{jn} Q_{ko}Q_{lp} C_{0\mathit{mnop}}, \end{aligned}$$

(46)

$$\begin{aligned} \varepsilon ^{\uptheta }_{ij} &= Q_{ik} Q_{jl} \boldsymbol{\varepsilon}^{\uptheta }_{0kl}. \end{aligned}$$

(47)

Expectation values of material properties can therefore be reduced to expectation values of orientation tensors. For a property $A$ with tensorial order $n$, this can be written as

$$ \langle A\rangle = \langle \boldsymbol{Q}\star A_{0}\rangle = \langle \boldsymbol{Q}^{\star n}\rangle [A_{0}], $$

(48)

where $\langle \boldsymbol{Q}^{\star n}\rangle$ takes the component form:

$$\begin{aligned} \langle \boldsymbol{Q}^{\star 0}\rangle [A_{0}] &= A_{0}, \end{aligned}$$

(49)

$$\begin{aligned} \langle \boldsymbol{Q}^{\star 2}\rangle [ \boldsymbol{A}_{0}] &= \int p(\boldsymbol{Q}) Q_{ik} Q_{jl} \mathop{}\!\mathrm{d}\boldsymbol{Q}A_{0kl}, \end{aligned}$$

(50)

$$\begin{aligned} \langle \boldsymbol{Q}^{\star 4}\rangle [{\mathbb{ A}}_{0}] &= \int p(\boldsymbol{Q}) Q_{im} Q_{jn} Q_{ko}Q_{lp} \mathop{}\! \mathrm{d}\boldsymbol{Q}A_{0\mathit{mnop}}. \end{aligned}$$

(51)

Determining these orientation integrals is therefore sufficient to calculate all expectation values. Thus it is possible to simplify the abbreviations of Sect. 3.1 for single-phase polycrystals by replacing all expectation values:

$$ w_{0} = \boldsymbol{\varepsilon}^{\uptheta }\cdot {\mathbb{ C}}[ \boldsymbol{\varepsilon}^{\uptheta }] - \boldsymbol{\varepsilon}_{\mathrm{{s}}} \cdot {\mathbb{ C}}_{+} \boldsymbol{\varepsilon}_{\mathrm{{s}}} - ( \boldsymbol{\varepsilon}_{\mathrm{{s}}} - \boldsymbol{\varepsilon}_{\mathrm{{e}}}) \cdot ({\mathbb{ C}}_{-}^{-1} - {\mathbb{ C}}_{+}^{-1})^{-1}[\boldsymbol{\varepsilon}_{\mathrm{{s}}} - \boldsymbol{\varepsilon}_{\mathrm{{e}}}], $$

(52)

$${\mathbb{C}}_{+}=〈\mathbb{C}〉=〈{\mathit{Q}}^{\star 4}〉[{\mathbb{C}}_0],{\mathbb{C}}_{-}={〈{\mathbb{C}}^{-1}〉}^{-1}={(〈{\mathit{Q}}^{\star 4}〉[{\mathbb{C}}_0^{-1}])}^{-1},$$

(53)

$${\mathit{\epsilon }}_{\mathrm{e}}=〈{\mathit{\epsilon }}^{\mathrm{\Theta }}〉=〈{\mathit{Q}}^{\star 2}〉[{\mathit{\epsilon }}_0^{\mathrm{\Theta }}],{\mathit{\epsilon }}_{\mathrm{s}}={〈\mathbb{C}〉}^{-1}[〈\mathbb{C}[{\mathit{\epsilon }}^{\mathrm{\Theta }}]〉]={(〈{\mathit{Q}}^{\star 4}〉[{\mathbb{C}}_0])}^{-1}[〈{\mathit{Q}}^{\star 2}〉.{\mathbb{C}}_0[{\mathit{\epsilon }}_0^{\mathrm{\Theta }}]]$$

(54)

After this simplification, the rest of the solution follows similarly to the general linear thermoelastic case, as shown in box 12.

Box 12

MEM solution for globally anisotropic single-phase polycrystals

3.7.2 Globally Isotropic Single-Phase Polycrystals

Globally isotropic polycrystals are characterized by isotropic effective material properties, which occur alongside particular orientation integrals:

$$\begin{aligned} \langle {\mathbb{ Q}}^{\star 2k +1}\rangle &= 0, \quad k \in \mathbb{N}, \end{aligned}$$

(55)

$$\begin{aligned} \langle {\mathbb{ Q}}^{\star 0}\rangle &= 1, \end{aligned}$$

(56)

$$\begin{aligned} \langle {\mathbb{ Q}}^{\star 2}\rangle &= {\mathbb{ P}}_{1}, \end{aligned}$$

(57)

$$\begin{aligned} \langle {\mathbb{ Q}}^{\star 4}\rangle &= \frac{{\mathbb{ P}}_{1} \otimes {\mathbb{ P}}_{1}}{{\mathbb{ P}}_{1}\cdot {\mathbb{ P}}_{1}} + \frac{{\mathbb{ P}}_{2} \otimes {\mathbb{ P}}_{2}}{{\mathbb{ P}}_{2}\cdot {\mathbb{ P}}_{2}} + \frac{{\mathbb{ P}}_{3} \otimes {\mathbb{ P}}_{3}}{{\mathbb{ P}}_{3}\cdot {\mathbb{ P}}_{3}} = {\mathbb{ P}}_{1} \otimes {\mathbb{ P}}_{1} + \frac{1}{5} { \mathbb{ P}}_{2} \otimes {\mathbb{ P}}_{2} + \frac{1}{3} { \mathbb{ P}}_{3} \otimes {\mathbb{ P}}_{3}. \end{aligned}$$

(58)

The abbreviations from Sect. 3.1 can be further simplified:

$$\begin{aligned} {\mathbb{ C}}_{+}&= \langle \boldsymbol{Q}^{\star 4}\rangle [{ \mathbb{ C}}_{0}] = ({\mathbb{ P}}_{1} \cdot {\mathbb{ C}}_{0}) { \mathbb{ P}}_{1} + \frac{{\mathbb{ P}}_{2} \cdot {\mathbb{ C}}_{0}}{5} {\mathbb{ P}}_{2}, \end{aligned}$$

(59)

$$\begin{aligned} {\mathbb{ C}}_{-}&= (\langle \boldsymbol{Q}^{\star 4}\rangle [{ \mathbb{ C}}_{0}^{-1}])^{-1} = \frac{1}{{\mathbb{ P}}_{1} \cdot {\mathbb{ C}}_{0}^{-1}} { \mathbb{ P}}_{1} + \frac{5}{{\mathbb{ P}}_{2} \cdot {\mathbb{ C}}_{0}^{-1}} { \mathbb{ P}}_{2}, \end{aligned}$$

(60)

$$\begin{aligned} \boldsymbol{\varepsilon}_{\mathrm{{e}}} &= \langle \boldsymbol{Q}^{\star 2}\rangle [ \boldsymbol{\varepsilon}^{\uptheta }_{0}] = {\mathbb{ P}}_{1}[ \boldsymbol{\varepsilon}^{\uptheta }_{0}] = \mathrm{sph}( \boldsymbol{\varepsilon}^{\uptheta }_{0} ), \end{aligned}$$

(61)

$$\begin{aligned} \boldsymbol{\varepsilon}_{\mathrm{{s}}} &= (\langle \boldsymbol{Q}^{\star 4}\rangle [{\mathbb{ C}}_{0}])^{-1} [ \langle \boldsymbol{Q}^{\star 2}\rangle {\mathbb{ C}}_{0}[ \boldsymbol{\varepsilon}^{\uptheta }_{0}]] \\ &= \left (\frac{1}{{\mathbb{ P}}_{1} \cdot {\mathbb{ C}}_{0}} { \mathbb{ P}}_{1} + \frac{5}{{\mathbb{ P}}_{2} \cdot {\mathbb{ C}}_{0}} {\mathbb{ P}}_{2} \right ) [{\mathbb{ P}}_{1} {\mathbb{ C}}_{0}[ \boldsymbol{\varepsilon}^{\uptheta }_{0}]] \\ &= \frac{1}{{\mathbb{ P}}_{1} \cdot {\mathbb{ C}}_{0}} { \mathbb{ P}}_{1} {\mathbb{ C}}_{0}[\boldsymbol{\varepsilon}^{\uptheta }_{0}]. \end{aligned}$$

(62)

All other equations are equivalent to the general polycrystalline case.

3.7.3 Multi-Phase Polycrystals

Multi-phase polycrystals can be modeled by replacing the reference properties ${\mathbb{ C}}_{0}$ and $\boldsymbol{\varepsilon}^{\uptheta }_{0}$ with phase reference properties ${\mathbb{ C}}_{\alpha }$ and $\boldsymbol{\varepsilon}^{\uptheta }_{\alpha }$. Similarly, the orientation distribution function becomes:

$$ p_{1}^{\mathrm{{C}}}(\boldsymbol{Q}) = \sum _{\alpha =1}^{N} c_{\alpha }p_{1}^{\mathrm{{C}}}(\boldsymbol{Q}_{\alpha }), $$

(63)

which leads to the abbreviations becoming sums as well:

$$\begin{aligned} {\mathbb{ C}}_{+}&= \langle {\mathbb{ C}}\rangle = \sum _{\alpha =1}^{N} c_{\alpha }\langle \boldsymbol{Q}^{\star 4}\rangle _{\alpha }[{ \mathbb{ C}}_{\alpha }], \end{aligned}$$

(64)

$$\begin{aligned} {\mathbb{ C}}_{-}&= \langle {\mathbb{ C}}^{-1}\rangle ^{-1} = \left (\sum _{\alpha =1}^{N} c_{\alpha }\langle \boldsymbol{Q}^{ \star 4}\rangle _{\alpha }[{\mathbb{ C}}_{\alpha }^{-1}]\right )^{-1}, \end{aligned}$$

(65)

$$\begin{aligned} \boldsymbol{\varepsilon}_{\mathrm{{e}}} &= \langle \boldsymbol{\varepsilon}^{\uptheta }\rangle = \sum _{\alpha =1}^{N} c_{\alpha }\langle \boldsymbol{Q}^{\star 2}\rangle _{\alpha }[ \boldsymbol{\varepsilon}^{\uptheta }_{\alpha }], \end{aligned}$$

(66)

$$\begin{aligned} \boldsymbol{\varepsilon}_{\mathrm{{s}}} &= \langle {\mathbb{ C}} \rangle ^{-1}\langle {\mathbb{ C}}[\boldsymbol{\varepsilon}^{\uptheta }]\rangle = \left (\sum _{\alpha =1}^{N} c_{\alpha }\langle \boldsymbol{Q}^{\star 4}\rangle _{\alpha }[{\mathbb{ C}}_{ \alpha }]\right )^{-1} \left [\sum _{\beta =1}^{N} c_{\beta }\langle \boldsymbol{Q}^{\star 2}\rangle _{\beta }{\mathbb{ C}}_{\beta }[ \boldsymbol{\varepsilon}^{\uptheta }_{\beta }] \right ]. \end{aligned}$$

(67)

The solution follows from the single-phase polycrystalline case in box 12 by substituting Equation 63, leading to the changed equations shown in box 13, with all other equations remaining the same.

Box 13

Modified formulae for multi-phase polycrystals

4 Examples

4.1 Computational Setup

As a comparison to the MEM results, full-field simulations were performed with solvers based on the fast Fourier transform (FFT). The methodology is based on the approach of Moulinec and Suquet [17]. The conjugate gradient method proposed by Zeman et al. [27] is used as it is particularly suited for linear problems. By relying upon the staggered grid discretization [22], so-called ringing artifacts or Gibbs’ oscillations which might impact the statistics of the results are minimized. All calculations were carried out in Python 3 using the FFTW library [7]. Further details on the FFT-based solver used can be found in Wicht et al. [25], [26].

The FFT simulations were performed on cubic representative volume elements with a resolution of 256 voxels on each side.

4.2 Isotropic Al-Al₂O₃

In this section, the metal matrix composite Al-Al₂O₃ is investigated using the MEM and FFT simulations for comparison. Aluminum and ceramic are modeled as discrete isotropic elastic phases with parameters according to Table 2.

Table 2 Material parameters for the Al-Al₂O₃ composite according to Agrawal and Sun [1]

As shown in the formulae of Sect. 3.5, in the case of only two phases, strains caused by effective stresses and strains caused by eigenstrains are independent of each other. First, thermal eigenstrains are neglected and an effective mean stress of $\bar{\boldsymbol{\sigma }}^{\mathrm{I}}= 100~\mbox{MPa}$$\boldsymbol{e}_{1} \otimes \boldsymbol{e}_{1}$ is applied to the composite. The value and direction of this loading are chosen arbitrarily, yet representative for all unidirectional stresses because of the MEM’s linearity and the isotropy of the material.

Three different microstructures with an Al₂O₃ volume fraction of 30% as shown in Fig. 2 are used in the FFT simulations. Shown in the middle is a polycrystal microstructure generated via Voronoi tessellation. 8000 grains are each randomly chosen to contain either of the two phases, yielding a co-continuous structure. Two additional microstructures are metal matrix inclusion structures, on the left with 300 spherical inclusions of equal radius and on the right with cylindrical inclusions with a ratio of length to diameter of 10, in the following referred to as fibers. The numerically computed effective stiffness tensor of all structures are isotropic and similar to each other, as seen in Table 3. The anisotropic error is defined as the Frobenius norm of the difference between the effective stiffness and its projection into the space of isotropic fourth-order tensors; more details may be found in Sect. 4.4.

Fig. 2

Al₂O₃ phase of microstructures used for the FFT reference solutions. From left to right: spherical inclusions, two-phase polycrystal, fiber inclusions

Table 3 Numerical effective stiffnesses and anisotropic error for three different Al-Al₂O₃ microstructures shown in Fig. 2

The distributions of phase stresses along the $\boldsymbol{e}_{1} \otimes \boldsymbol{e}_{1}$-direction for the two phases according to the MEM are shown in Fig. 3. Two Gaussian distributions with means as given by the exact analytical solution for two phases can be seen. Statistical moments with orders exceeding two vanish. For comparison, the distributions according to the FFT simulations are shown in Fig. 4. For the polycrystal and the fiber inclusions, higher statistical moments apparently do not vanish. In the case of the polycrystal and the fiber inclusions, the variance is high compared to the MEM. The variance of the spherical inclusions is comparatively low.

Fig. 3

Phase-wise MEM $\sigma _{11}$ distributions for Al-Al₂O₃ under $\bar{\sigma }_{11}$-load for three different microstructures

Fig. 4

Phase-wise MEM $\sigma _{11}$ distributions for Al-Al₂O₃ under $\bar{\sigma }_{11}$-load for polycrystal microstructure; comparison with polycrystal, sphere and fiber inclusion FFT simulation

The previous diagrams show the marginalized distribution of only one component of the stress. Even in unidirectional mean stress states, however, the variances are multidirectional. To depict the correlation between different components, it is expedient to use a more expressive type of diagram derived from the covariance ellipses known from multivariate statistics. An annotated example is given in Fig. 5. The main feature of the diagram is the covariance ellipse, which is an isosurface of the marginalized probability density $p(x_{1}, x_{2})$, chosen such that the ellipse encloses a certain cumulative probability. In this paper, the ellipses are chosen such that there is a 95% confidence that any stress state falls into the ellipse; put differently, the 95% most probable stress states are found within the ellipse. By the size of the ellipse, one can estimate how wide the stresses are spread. Its tilt reflects the correlation between the components, with the ellipse of an uncorrelated distribution being aligned with the coordinate system. As the multivariate Gaussian distribution is symmetric, the center of the ellipse corresponds to the mean values of the distribution. However, this ellipse illustrates only the marginal distribution, discarding information about correlations with components which are not shown. To ameliorate this shortcoming, the principal axes of the covariance ellipsoid of the full distribution are projected into the diagram. If certain dimensions are uncorrelated, the projected axes degenerate into points in the center; if all correlations concerning the components shown in the diagram are accounted for, the projected axes overlap with the principal axes of the shown ellipse.

Fig. 5

Annotated example covariance ellipse of a marginal distribution $p(x_{1}, x_{2})$ of a full distribution $p(\boldsymbol{x})$, showing the information conveyed in this diagram type

In Fig. 6, the covariance ellipses for $\sigma _{11}$ and $\sigma _{22}$ are shown for the MEM solution and the three FFT simulations. Higher statistical moments of the FFT simulations are neglected in this diagram. The difference in effective stiffnesses between the various structures is evidenced by the distance between the midpoints of the ellipses. Furthermore, it becomes clear that the apparent higher variance of the polycrystal and fiber inclusion simulations compared to the MEM is present only in the direction of loading. As the stress covariance matrix of the MEM solution is proportional to the stiffness tensor in Mandel notation, the shape of the covariance ellipse is not dependent on the loading direction, whereas such a dependence is evidenced by the simulations. The principal axes of the MEM covariance ellipsoid are depicted, showing that one covariance is unaccounted for. Due to the materials’ isotropy, this covariance can only be that of the shown directions with the 33-direction.

Fig. 6

Phase-wise MEM covariance ellipses of $\sigma _{11}$ and $\sigma _{22}$ for Al-Al₂O₃ under $\bar{\sigma }_{11}$-load; comparison with polycrystal, sphere and fiber inclusion FFT simulation

In lieu of a six-dimensional diagram, the representation of the covariance hyperellipsoid is separated into two projections: one which captures the covariances between the normal stresses, whereas the other depicts the shear stresses. As the covariances of the MEM solution are proportional to the phase stiffnesses, this separation always captures the full MEM information for a stiffness matrix without normal-shear correlations, i.e. for isotropic, cubic, orthotropic, trigonal, tetragonal, and transversely isotropic phases in the coordinate systems where their symmetry is apparent. This approach leads to Fig. 7. For the multivariate normal distributions considered here, this style of diagram has higher information content than the two-dimensional probability density plots. In this particular case, the FFT covariance is anisotropic not only in the normal, but also in the shear directions. A loading in 11-direction causes not only higher normal stress covariances in that direction, but also increases the stress covariances in the 12- and 13-direction. Otherwise, the shear covariances are similar for both methods considered.

Fig. 7

MEM covariance ellipsoid projections into normal stress space (left) and shear stress space (right) along with marginal distributions for Al-Al₂O₃ under $\bar{\sigma }_{11}$-loading; comparison with polycrystal FFT simulation

In the Al-Al₂O₃ composite, differing thermal expansion coefficients cause spherical eigenstrains of .355% and .099%, respectively, at a process temperature of $\delta \theta = 150\,{}^{\circ }\mbox{C}$ [1]. The resulting one-dimensional distributions are shown in Fig. 8. The MEM solution shows marked differences to each of the FFT solutions. In particular, higher statistical moments are not negligible in the FFT solutions and the variances of the MEM solution are considerably larger.

Fig. 8

Phase-wise MEM $\sigma _{11}$ distributions for Al-Al₂O₃ under thermal strain

The covariance ellipsoid projections are depicted in Fig. 9. As the purely thermal load case is isotropic, the covariance ellipsoids have the same width in each of the three directions. However, the MEM overestimates the degree of correlation between different stress directions, as the difference between the shapes of the ellipsoids shows.

Fig. 9

MEM covariance ellipsoid projections into normal stress space (left) and shear stress space (right) along with marginal distributions for Al-Al₂O₃ under thermal strain; comparison with polycrystal FFT simulation

4.3 Parametric Study

In the following, the sensitivity of the MEM with regard to its input parameters will be investigated. All simulations are performed on the polycrystal structure described in the previous section, as the polycrystal is a structure for which arbitrary phase volume fractions may be realized.

Both variances and mean values of the stress distribution are affected by variations in the Al₂O₃ volume fraction $c_{2}$. In the degenerate cases $c_{2} = 0$ and $c_{2} = 1$, the material consists of a single phase; at $c_{2} = 0.5$, a maximum of stress variances is reached. This is illustrated in Fig. 10 on the left. The diagram depicts the mean of the probability distribution, adding and subtracting a standard deviation to convey a sense of the variance. This may be understood as a one-dimensional cut through a covariance ellipsoid for each value of $c_{2}$. On the right, the same diagram is shown for the 22-component of stress. The variances in the 11- and 22-diagrams are proportional to each other. The MEM prediction that the shape of the covariance ellipsoids is the same for any volume fraction holds true for the FFT.

Fig. 10

Phase-wise $\sigma _{11}$ (left) and $\sigma _{22}$ (right) mean (solid line) and spread of one standard deviation (area) for polycrystalline Al-Al₂O₃ structure of varying Al₂O₃ content under $\bar{\sigma }_{11}$ load, calculated with MEM and FFT

In general, the MEM becomes less reliable as the contrast in phase stiffnesses increases. This is illustrated in Fig. 11 for $\sigma _{11}$ and $\sigma _{22}$, where the phase contrast parameter $\alpha = E_{2}/E_{1}$ is varied from 1 to 11. According to [14], the simple MEM is usable for phase contrasts up to about 5; an advanced method incorporating some structure assumptions is detailed in [13]. While there is no visible qualitative change at that point, the divergence in the variances does increase with rising $\alpha$. As in the study of $c_{2}$, the shape of the covariance ellipsoids remains roughly the same. The effect of $\alpha$ on the covariance, which is markedly different for the two phases, is predicted approximately by the MEM.

Fig. 11

Phase-wise $\sigma _{11}$ (left) and $\sigma _{22}$ (right) mean (solid line) and spread of one standard deviation (area) for polycrystalline material of varying elastic contrast under $\bar{\sigma }_{11}$-load, calculated with MEM and FFT

Finally, we consider the influence of material structure. While the maximum entropy method does not incorporate any explicit structural information beyond the one-point probability function of material properties $p_{1}^{\mathrm{{C}}}$, some amount of higher-order structural information is contained in the effective properties. Instead of using the numerical effective stiffness, the effective stiffness $\bar{{\mathbb{ C}}}$ is parameterized as a weighted average of the Voigt and Reuss bounds $\bar{{\mathbb{ C}}}= \delta {\mathbb{ C}}_{+} + (1-\delta ){ \mathbb{ C}}_{-}$, $\delta \in (0, 1)$, which are strict bounds for the elastic behavior. As Fig. 12 shows, the method is highly sensitive to variations of the effective stiffness tensor. On reaching the Reuss or Voigt bound, the variances of both phases vanish.

Fig. 12

MEM phase-wise $\sigma _{11}$ (left) and $\sigma _{22}$ (right) mean (solid line) and spread of one standard deviation (area) for Al-Al₂O₃ with varying $\bar{{\mathbb{ C}}}$ under $\bar{\sigma }_{11}$-load

4.4 Uniform Globally Isotropic Polycrystal

As an example for polycrystals, we consider a cubic polycrystal with a uniform ODF, yielding an isotropic effective behaviour. The polycrystal is composed of copper single crystals of different orientation. Copper has been chosen for its significant degree of anisotropy. Its elastic constants are given as $C_{1111} = 170.2~\mbox{GPa}$, $C_{1122} = 114.9~\mbox{GPa}$, and $C_{1212} = 61.0~\mbox{GPa}$ [23], which are the components of the stiffness tensor in the crystal lattice coordinate system. As in the previous example, an unidirectional stress of $\bar{\boldsymbol{\sigma}}^{\mathrm{{I}}}= 100~\mbox{MPa}\ \boldsymbol{e}_{1} \otimes \boldsymbol{e}_{1}$ is applied to the material.

An FFT simulation is used as a reference. The microstructure is generated via Voronoi tessellation, with 8000 crystals distributed over a cube 256 voxels wide in each direction. The orientations were sampled from a uniform distribution, yielding a numerical effective stiffness which is only slightly anisotropic. The relative anisotropic error was found to be negligible:

$$\frac{\parallel {\overline{\mathbb{C}}}_{\mathrm{iso}}-\overline{\mathbb{C}}\parallel }{\parallel \overline{\mathbb{C}}\parallel }<0.05\mathrm{\%},{\overline{\mathbb{C}}}_{\mathrm{iso}}=(\overline{\mathbb{C}}\cdot {\mathbb{P}}_1){\mathbb{P}}_1+\frac{1}{5}(\overline{\mathbb{C}}\cdot {\mathbb{P}}_2){\mathbb{P}}_2.$$

(68)

One should note that a similar definition of the anisotropic error using the compliance tensor yields slightly different values. The numerical stiffness is used as input of the MEM to ensure comparability between the methods. However, where the FFT result is calculated with 8000 grain orientations drawn from a uniform distribution, the MEM results can be calculated using the uniform distribution itself.

In Fig. 13, the distribution of $\sigma _{11}$ is shown for both methods. For the FFT, the diagram was compiled from the available full-field information, while in the case of the MEM, Monte Carlo integration was used to calculate this marginalized distribution. The overall variance of the distribution is on the order of 10 MPa in both cases, which is not negligible, yet both methods yield quite similar distributions.

Fig. 13

Probability distribution of $\sigma _{11}$ for isotropic polycrystalline copper, calculated with MEM and FFT

For an arbitrarily chosen orientation, the normal and shear projections of the covariance ellipsoid of stress distribution are depicted in Fig. 14. Note that the coordinate system is rotated to coincide with the local crystal coordinate system. As a result, the covariance ellipsoid of the MEM exhibits no shear-normal correlations. The shear-normal correlations of the FFT are not represented. The variances of the single grain are on the same order of magnitude as those depicted in Fig. 13, and their local shape is similar in both MEM and FFT. A comparison with the similarly rotated effective stress $\bar{\sigma } = {(5.6, 26.5, 67.9, -60.0, -27.6, 17.2)}^{ \mathsf{T}}$ MPa reveals that the single grain mean stresses are far from the effective stress of the polycrystal, yet FFT results are predicted satisfactorily by the MEM. Interpreting the differences in covariances is not straightforward. While the MEM suffers from some approximations, the FFT results are taken from only one single grain. For the purpose of calculating statistical moments of a given orientation, they are suboptimal because the single grain interacts only with a few grains in the realization considered. By contrast, the MEM result is not dependent on a single realization and might be more representative.

Fig. 14

MEM covariance ellipsoid projections into normal stress space (left) and shear stress space (right) along with marginal distributions for a single copper grain under $\bar{\sigma }_{11}$-load, rotated into the crystal coordinate system; comparison with polycrystal FFT simulation

5 Conclusions

In the present work, a framework for a probabilistic estimate for stress and strain fields of random heterogeneous linear thermoelastic materials was considered. Special cases for various material classes were derived, including:

• solutions for linear thermoelastic and linear elastic materials with discrete phases,

• of these, specialized solutions for all combinations of global and local isotropy and anisotropy,

• and finally polycrystal solutions in terms of ODFs, in particular multi-phase polycrystals, uniform polycrystals, and globally isotropic uniform polycrystals.

While discrete phases seem to be a restriction of the method’s capabilities, they are the most common, the most easily numerically implemented, and can model all other cases approximately.

For an Al-Al₂O₃ composite, a parameter study thereof and a copper polycrystal, MEM results and FFT simulations were compared. The main conclusions of these studies were:

• The method’s applicability is, particularly for high elastic phase contrasts, highly dependent on the microstructure of the material considered.

• For polycrystals, the method is well suited, as it was shown capable of predicting both the overall stress distributions and the grain mean stresses of copper.

• As the MEM covariances are proportional to the local stiffness, their shape cannot be directly dependent on the load direction. The general magnitude of the covariances is nonetheless often comparable with FFT results.

An expressive type of stress state diagram analogous to the covariance ellipses known from multivariate statistics was used to visualize these results.

As a future development, the MEM may be refined by selectively including structural information, which would broaden its applicability to more microstructures. Furthermore, an in-depth study of the covariance shapes and the possibility of incorporating dependencies on loading directions may be warranted. The method presented in this paper may provide impetus for nonlinear homogenization methods that are typically based on mean values of quantities involves, and may profit from accurate estimates of strain and stress field variances.