A statistical framework for generating microstructures of two-phase random materials: application to fatigue analysis

TL;DR

This paper presents a flexible statistical microstructure model for two-phase materials using Gaussian random fields, enabling improved fatigue analysis and microstructure design through efficient Bayesian parameter estimation and Monte Carlo simulations.

Contribution

It introduces a novel surrogate microstructure model based on Gaussian random fields with controllable parameters, enhancing realism and applicability in fatigue analysis.

Findings

Model accurately captures complex microstructural features.

Efficient Bayesian estimation from real images.

Enables uncertainty quantification in fatigue predictions.

Abstract

Random microstructures of heterogeneous materials play a crucial role in the material macroscopic behavior and in predictions of its effective properties. A common approach to modeling random multiphase materials is to develop so-called surrogate models approximating statistical features of the material. However, the surrogate models used in fatigue analysis usually employ simple microstructure, consisting of ideal geometries such as ellipsoidal inclusions, which generally does not capture complex geometries. In this paper, we introduce a simple but flexible surrogate microstructure model for two-phase materials through a level-cut of a Gaussian random field with covariance of Mat\'ern class. Such parametrization of the covariance function allows for the representation of a few key design parameters while representing the geometry of inclusions in a more general setting for a large class of random heterogeneous two-phase media. In addition to the traditional morphology descriptors such as porosity, size and aspect ratio, it provides control of the regularity of the inclusions interface and sphericity. These parameters are estimated from a small number of real material images using Bayesian inversion. An efficient process of evaluating the samples, based on the Fast Fourier Transform, makes possible the use of Monte-Carlo methods to estimate statistical properties for the quantities of interest in a given material class. We demonstrate the overall framework of the use of the surrogate material model in application to the uncertainty quantification in fatigue analysis, its feasibility and efficiency, and its role in the microstructure design.

Equations

$$\chi (\mathbf{x};\omega )=\{\begin{array}{ll}1,& \mathbf{x}\text{ in inclusion},\\ 0,& \mathbf{x}\text{ in matrix},\end{array}\mathbf{x}\in D\subset {\mathbb{R}}^d,\omega \in \Omega ,$$

$$\kappa (\mathbf{x};\omega )={\kappa }_I\cdot \chi (\mathbf{x};\omega )+{\kappa }_M\cdot (1-\chi (\mathbf{x};\omega )),\mathbf{x}\in D,\omega \in \Omega ,$$

$$\chi (\mathbf{x};\omega )=\{\begin{array}{ll}1,& \text{if }|m(\mathbf{x};\omega )|\ge \tau ,\text{inclusions},\\ 0,& \text{if }|m(\mathbf{x};\omega )|<\tau ,\text{matrix},\end{array}$$

$$\kappa (\mathbf{x})=\frac{{\kappa }_I+{\kappa }_M}{2}+\frac{{\kappa }_I-{\kappa }_M}{2}\cdot \mathrm{sign}\left(|m(\mathbf{x})|-\tau \right).$$

$$S_1={\varphi }_0$$

$$S_2(\mathbf{x},\mathbf{y})$$

$$=2{\varphi }_0-4T\left(\tau ,\sqrt{\frac{1-\mathcal{C}(\mathbf{x},\mathbf{y})}{1+\mathcal{C}(\mathbf{x},\mathbf{y})}}\right)-4T\left(\tau ,\sqrt{\frac{1+\mathcal{C}(\mathbf{x},\mathbf{y})}{1-\mathcal{C}(\mathbf{x},\mathbf{y})}}\right),$$

$$S_2(\mathbf{x},\mathbf{x})={\varphi }_0,\underset{\Vert \mathbf{x}-\mathbf{y}\Vert \to \infty }{\lim }{S}_2(\mathbf{x},\mathbf{y})={\varphi }_0^2.$$

$$\tau ({\varphi }_0)=\sqrt{2}{\mathrm{erf}}^{-1}(1-{\varphi }_0),$$

$$S_2(\mathbf{x},\mathbf{y})=S_2(\Vert \mathbf{x}-\mathbf{y}\Vert ),\mathcal{C}(\mathbf{x},\mathbf{y})=\mathcal{C}(\Vert \mathbf{x}-\mathbf{y}\Vert ).$$

$$\mathcal{C}(\boldsymbol{x},\boldsymbol{y})=\sigma^{2}\,\mathcal{M}_{\nu}\left(\sqrt{2\nu}\cdot r\right),\qquad r=\sqrt{(\boldsymbol{x}-\boldsymbol{y})\boldsymbol{\Theta}^{\scalebox{0.7}[1.0]{$-$}1}(\boldsymbol{x}-\boldsymbol{y})},$$

$${\mathcal{M}}_{\nu }(x)=\frac{1}{2^{\nu -1}\Gamma (\nu )}{x}^{\nu }{\mathcal{K}}_{\nu }(x)$$

$$\mathbf{\Theta }={\mathbf{R}}^t\left[\begin{array}{ccc}{\ell }_1^2& & \\ & \ddots & \\ & & {\ell }_d^2\end{array}\right]\mathbf{R},$$

$${\left(\mathrm{Id}-\frac{1}{2\nu }\nabla \cdot (\mathbf{\Theta }\nabla )\right)}^{\frac{1}{2}(\nu +\frac{d}{2})}m(\mathbf{x};\omega )=\eta \mathcal{W}(\mathbf{x};\omega ),$$

$${\eta }^2=\frac{(2\pi {)}^{\frac{d}{2}}\sqrt{\det \mathbf{\Theta }}\Gamma (\nu +d/2)}{{\nu }^{\frac{d}{2}}\Gamma (\nu )}.$$

$$m(\mathbf{x};\omega )={\mathcal{F}}^{-1}\left\{\hat{G}(\mathbf{\xi })\cdot \hat{\mathcal{W}}(\mathbf{\xi };\omega )\right\}(\mathbf{x}),$$

$$|\hat{G}(\mathbf{\xi }){|}^2=\frac{(2\pi {)}^{\frac{d}{2}}\sqrt{\det \mathbf{\Theta }}\Gamma (\nu +d/2)}{{\nu }^{\frac{d}{2}}\Gamma (\nu )}\cdot {\left(1+\frac{1}{2\nu }\mathbf{\xi }\cdot \mathbf{\Theta }\mathbf{\xi }\right)}^{-(\nu +\frac{d}{2})}.$$

$$\mathrm{Sphericity}=\frac{{\pi }^{1/3}(6V_{pore}{)}^{2/3}}{0.833\cdot S_{pore}},r={\left(\frac{V_{pore}}{V_{total}}\right)}^{\frac{1}{3}},$$

$$S_2^{data}({\mathbf{r}}_n)=\frac{1}{N_n}\sum _k\chi ({\mathbf{x}}_k)\chi ({\mathbf{x}}_k+{\mathbf{r}}_n),$$

$$f\left(\frac{1}{N}\sum _{k=1}^N{z}_k\right)\approx \frac{1}{\sqrt{2\pi N^{-1}p\cdot (1-p)}}\mathrm{Exp}\left(-\frac{{\left(\frac{1}{N}{\sum }_{k=1}^N{z}_k-p\right)}^2}{2N^{-1}p\cdot (1-p)}\right),\forall {\mathbf{r}}_n.$$

$$\mathrm{Likelihood}=\prod _n\frac{{\mathrm{Constant}}_n}{\sqrt{S_2(r_n)\cdot (1-S_2(r_n))}}\mathrm{Exp}\left(-\frac{{\left(S_2^{data}({\mathbf{r}}_n)-S_2(r_n)\right)}^2}{2N_n^{-1}{S}_2(r_n)\cdot (1-S_2(r_n))}\right),$$

$$S_2^{model}(r;\mathbf{\xi })=2{\varphi }_0-4T\left(\tau ,\sqrt{\frac{1-\mathcal{C}(r;\mathbf{\xi })}{1+\mathcal{C}(r;\mathbf{\xi })}}\right)-4T\left(\tau ,\sqrt{\frac{1+\mathcal{C}(r;\mathbf{\xi })}{1-\mathcal{C}(r;\mathbf{\xi })}}\right)$$

$${\mathbf{\xi }}_{MLE}=\arg \underset{\mathbf{\xi }\in \Xi }{\max }\mathrm{Likelihood}(\mathbf{\xi }),$$

$$\boldsymbol{\Sigma}=-\left.\left(\frac{\operatorname{d}\!^{2}\log\operatorname{Likelihood}(\boldsymbol{\xi})}{\operatorname{d}\!\,(\log\boldsymbol{\xi})^{2}}\right)^{\scalebox{0.7}[1.0]{$-$}1}\right|_{\boldsymbol{\xi}=\boldsymbol{\xi}_{MLE}}.$$

$$\mathbf{u}(\mathbf{x})=\mathbf{E}\cdot \mathbf{x}+\tilde{\mathbf{u}}(\mathbf{x}),$$

$$\mathbf{\epsilon }(\mathbf{x})=\frac{1}{2}\left(\nabla \mathbf{u}(\mathbf{x})+\nabla \mathbf{u}(\mathbf{x}{)}^t\right),$$

$$\mathbf{\epsilon }(\mathbf{x})=\mathbf{E}+\tilde{\mathbf{\epsilon }}(\mathbf{x}),\tilde{\mathbf{\epsilon }}(\mathbf{x})=\frac{1}{2}\left(\nabla \tilde{\mathbf{u}}(\mathbf{x})+\nabla \tilde{\mathbf{u}}(\mathbf{x}{)}^t\right),$$

$$\mathbf{\sigma }(\mathbf{x})=\mathbf{C}(\mathbf{x})\cdot \mathbf{\epsilon }(\mathbf{x})=\left(K(\mathbf{x})-\frac{2}{d}G(\mathbf{x})\right)\mathrm{tr}(\mathbf{\epsilon }(\mathbf{x}))\mathbf{Id}+2G(\mathbf{x})\mathbf{\epsilon }(\mathbf{x}),$$

$$K(\mathbf{x})=K_I\chi (\mathbf{x})+K_M(1-\chi (\mathbf{x})),G(\mathbf{x})=G_I\chi (\mathbf{x})+G_M(1-\chi (\mathbf{x})),$$

$$\nabla \cdot \mathbf{\sigma }(\mathbf{x})=0\text{in }D,\underset{D}{\int }\mathbf{\sigma }(\mathbf{x})d\mathbf{x}=\mathbf{\Sigma },$$

$$\displaystyle\nabla\cdot\bigl{(}\boldsymbol{C}(\boldsymbol{x})\cdot\tilde{\boldsymbol{\varepsilon}}(\boldsymbol{x})\bigr{)}=-\nabla\cdot\bigl{(}\boldsymbol{C}(\boldsymbol{x})\cdot\boldsymbol{E}\bigr{)}\quad\text{in }D,$$

Full text

{AMS}

82D30, 60G60, 62M40, 65C60, 62F15

1 Introduction

Real materials, employed in engineering applications or that occur in nature, possess random heterogeneous microstructures, which play a crucial role in their macroscopic behavior and in determining their effective properties [36, 37, 15, 16, 35]. The determination of how particular microstructures lead to effective macroscale properties has been a major goal of material science for decades. The particular features of the microstructure of a material, as for example porosities, cracks, polycrystalline shape and texture, etc., have an important impact on the lifetime of an industrial structure and highly depend on the choice of manufacturing process. The standard paradigm is to split the design analysis of structures in two steps. First, a global mechanical analysis is performed on the complete structure using the constitutive law of the material, and second, a local fatigue analysis based on the global stabilized stress-strain response is applied to a local elementary volume element of material. Standard fatigue analysis is generally based on a phenomenological formulae and only characterizes the average response of the material, neglecting the microstructure variability at the material. A micromechanical model representing the microstructure as a plastic inclusion has been initially proposed by Orowan [57] and interpreted and developed later for polycrysalline solids in [24]. Recently, methods have been proposed to analyze the statistical output of the microstructure and to estimate not only the average statistical lifetime, but also the spread of the lifetime due to the variability of the microstructure [32].

A common approach to modeling random multiphase heterogeneous materials is to develop so-called surrogate models that ideally mimic key macroscale features of the material. These models approximate the geometric distribution of phases, and effective properties are calculated using homogenization, often implemented using Monte-Carlo methods. Such surrogates for two-phase materials have been used to study and to characterize properties of the porous media, e.g., permeability [11] or effective constitutive models [33, 30, 52, 23, 22, 47], in fatigue and damage analysis [28, 49], in microstructure design [7]. However, these models usually employ simple microstructure, consisting only of ideal geometries such as ellipsoidal inclusions. Understandably, these idealized surrogate models have limited use as they generally do not capture complex geometries and distributions of phases such as inclusions.

Alternatively, realistic microstructures can be obtained experimentally by extracting and analyzing samples from real materials, e.g., Computed Tomography (CT) or Scanning Electron Microscopy (SEM). These methods are time-consuming and expensive, possibly destructive, and cannot be deployed on a large database. In order to complete the experimental databases, one can use algorithms to generate artificial microstructures, based on the statistical information from only a small number of original samples. Such stochastic reconstruction methods provide effective tools for efficient reproduction of heterogeneous microstructures [62, 73, 46, 56]. Recently developed methods based on convolutional neural networks [53, 14] fall into this category. Although these methods can lead to an accurate microstructure reconstruction, their mathematical and statistical structures are often based on heuristic arguments.

In this paper, we first propose a simple but flexible surrogate microstructure model for two-phase materials. Key design parameters of the microstructure are provided, while the geometry of inclusions or pores is represented in a general setting, which permits a close representation of the microstructure of materials. In a second stage, we use the surrogate microstructure as the representative volume element (RVE) for a fatigue analysis of the Dang Van type [24, 59, 20, 18]. More precisely, we focus on the high cycle fatigue regime where both structure and microstructure are in an elastic shakedown regime under an external cyclic loading. The novelty is two fold: on one hand, we introduce realistic microstructures in the analysis, and on the other hand, we propose tools for a statistical analysis of the results at the local scale of the elementary volume element.

We define our surrogate material model through a level-cut of a Gaussian random field [61, 5, 67, 43, 64], called the intensity field. Such a model is entirely controlled by the mean and the covariance function of this field. The choice of the covariance defines the morphology of the microstructure. Thus, considering the covariances from a given class constrains our design parameter space. In our model we consider covariances of Matérn class [71, 48, 34]. In addition to the traditional morphology descriptors such as porosity, size and aspect ratio, this covariance class provides control of the smoothness of the inclusion interface. This regularity parameter is also related to the sphericity of pores or inclusions.

While defined only by a few ”design” parameters, such a surrogate model is able to mimic a large class of random heterogeneous two-phase media. These parameters can be estimated from a small number of real material images by approximating the image statistical descriptors by the surrogate ones. Here, we use Bayesian inversion to find the probability distributions of the design parameters, from which they are then drawn for each material sample.

The intensity field is constructed by convolving the Matérn-type covariance with the white noise. In the case of statistically homogeneous media, this can be done by Fast Fourier Transform. Moreover, the Fourier Transform of the Matérn covariance is given in a closed analytical form. This leads to a very efficient sampling process, which makes possible the use of Monte-Carlo methods to evaluate statistical properties for the quantities of interest in a given material class.

An important goal of this work is to demonstrate the overall framework of the use of the surrogate material model in applications to homogenization and fatigue analysis, its feasibility and efficiency, and its role in the microstructure design.

The paper is structured as follows. First, we describe in detail the properties and assumptions underlying the construction of a class of surrogate material models in Section 2. We present the level-set of Gaussian field model (Subsection 2.1), the Matérn covariance (Subsection 2.2) and the microstructure reconstruction technique from real material images (Subsection 2.3). In Section 3, we give the formulation of the associated linear elasticity problem together with the description of quantities of mechanical interest used in fatigue analysis. The simulation results are taken up in Section 4, including analysis and prediction of damage evolution and failure (Subsection 4.1), and a discussion of the use of the resulting models for the design of microstructure employing the process parameters (Subsection 4.2). We present final discussions and conclusions in Section 5.

2 Surrogate material model

Let us consider a random two-phase composite material in a bounded domain $D\subset\mathbb{R}^{d}$, $d=2,3$, the phases representing a matrix with inclusions (e.g., pores, precipitates, etc.). For a particular microstructure realization, the material phases are defined by the characteristic function

[TABLE]

where $\Omega$ is the space of the material samples (as in [68]), such that each sample point $\omega\in\Omega$ corresponds to a realization of a spatial random field $\chi(\boldsymbol{x};\omega)$, $\boldsymbol{x}\in D$. Thus, for a particular realization $\omega$, the spatial distribution of a material property $\kappa$ over $D$ can be written as

[TABLE]

where $\kappa_{I}$ and $\kappa_{M}$ are the corresponding properties of the inclusions and the matrix, respectively.

2.1 Level-cut of a Gaussian field

The characteristic function of the material phase can be expressed as a level-cut of an intensity field $m(\boldsymbol{x};\omega)$:

[TABLE]

$\omega\in\Omega$, where the level $\tau\geq 0$ controls the volume fractions of the inclusions. We define the intensity $m(\boldsymbol{x})$ as a zero-mean Gaussian random field [2] with covariance $\mathcal{C}(\boldsymbol{x},\boldsymbol{y})$. Suggested in earlier works of [13, 61, 5], the level-set model is widely used in porous media reconstruction [67, 43, 39] and in geometric inverse problems for interfaces [64, 26, 17].

Given (3), the formula (2) can be rewritten as

[TABLE]

Let us denote the standard deviation of $m(\boldsymbol{x};\omega)$ by $\sigma=\sqrt{\mathcal{C}(\boldsymbol{x},\boldsymbol{x})}$. Then, since $\operatorname*{sign}\!\left(|m(\boldsymbol{x})|\!-\!\tau\right)$ $\!=\!\operatorname*{sign}\left(\frac{|m(\boldsymbol{x})|}{\sigma}-\frac{\tau}{\sigma}\right)$, we can consider the Gaussian field \sigma^{\scalebox{0.7}[1.0]{-}1}m(\boldsymbol{x}) instead of $m(\boldsymbol{x})$ with the level parameter $\frac{\tau}{\sigma}$. Thus, without loss of generality, we can assume $\mathcal{C}(\boldsymbol{x},\boldsymbol{x})=1$ and thus $\sigma=1$.

The first two moments of the random field $\chi(\boldsymbol{x},\omega)$ are given in closed form, similar to [6, 67], by the following lemma (see A for the proof).

Lemma 2.1**.**

Let $S_{1}=\operatorname{\mathbb{E}}\left[\chi(\boldsymbol{x})\right]$ and $S_{2}(\boldsymbol{x},\boldsymbol{y})=\operatorname{\mathbb{E}}\left[\chi(\boldsymbol{x})\,\chi(\boldsymbol{y})\right]$ be respectively the one- and two-point correlation functions [68] of the random field $\chi(\boldsymbol{x},\omega)$ defined by (3). Then, they can be written as

[TABLE]

where $T(\tau,x)=\frac{1}{2\pi}\int\limits_{0}^{x}e^{-\frac{\tau^{2}}{2}(t^{2}+1)}\frac{\operatorname{d}\!t}{t^{2}+1}$ is Owen’s T function [58, 60], and $\phi_{0}$ denotes the expected volume fraction of the inclusions (or porosity). It also holds

[TABLE]

From (5), it follows that the level-set parameter $\tau$ is uniquely defined through the given expected volume fraction $\phi_{0}$,

[TABLE]

where $\operatorname{erf}(x)$ denotes the Gauss error function.

{rmk}

Defining in (3) the phase through the absolute value of $m(\boldsymbol{x};\omega)$, see also [5, 6], allows us to distinguish the matrix and the inclusions for all values of the inclusions volume fraction $\phi_{0}$. This departs from the model in [61, 67, 39], where the phase is defined simply by the sign of $m(\boldsymbol{x};\omega)-\tau$, $\tau\in\mathbb{R}$, and the ”matrix” becomes ”inclusions” when $\phi_{0}>0.5$ (see Figure 1).

In this work, we primarily focus on statistically isotropic materials [68]. The isotropy constraint on $\mathcal{C}(\boldsymbol{x},\boldsymbol{y})$ is taken for simplicity, but it does not limit the generality of the approach. In this case, the covariance $\mathcal{C}(\boldsymbol{x},\boldsymbol{y})$ and, thus, the two-point correlation $S_{2}(\boldsymbol{x},\boldsymbol{y})$ are rotationally invariant (stationary) and depend only on the distance between the points:

[TABLE]

In what follows, we will write a stationary function of two arguments $\boldsymbol{x}$ and $\boldsymbol{y}$ as a function of the distance $r=\lVert\boldsymbol{x}-\boldsymbol{y}\rVert$.

2.2 Matérn covariance

The microstructure morphology of the above level-cut model is defined by the covariance function $\mathcal{C}(\boldsymbol{x},\boldsymbol{y})$ of the intensity field $m(\boldsymbol{x};\omega)$. Parametrization of this covariance provides the so-called ”design” parameters of the microstructure. In general, this covariance can have an infinite dimensional parameter space. In order to reduce the model and to have a limited number of design parameters, one can consider a projection of the material onto a reduced material class, that is, one may consider a parameterized class of covariance kernels. Here we use the Matérn covariance class [71, 48, 34], that is of the form

[TABLE]

with unit standard deviation $\sigma=1$ and Matérn kernel

[TABLE]

where $\Gamma$ and $\mathcal{K}_{\nu}$ denote the Euler Gamma function and the modified Bessel function of the second kind [1, 4, 70], respectively. The scalar parameter $\nu>0$ defines the differentiability (smoothness) of the field, while the second order tensor $\boldsymbol{\Theta}$ defines the shape of inclusions. In particular, $\boldsymbol{\Theta}=\ell^{2}\,\boldsymbol{\operatorname*{Id}}$ corresponds to an isotropic covariance with correlation length $\ell$. In general, $\boldsymbol{\Theta}$ can be written through a rotation and scaling as

[TABLE]

where the matrix $\boldsymbol{R}$ defines a rotation of axes in $\mathbb{R}^{d}$, and $\ell_{k}$, $k=1,\ldots,d$, are the correlation lengths in each axis direction.

The Matérn covariance model (11) is widely used in statistics [65, 29], geostatistics [51] and machine learning [72], and represents a large class of popular covariance kernels. In particular, when $\nu$ varies from $0.5$ to infinity, it presents a continuous family of kernels from the exponential function to the square exponential (Gaussian), respectively.

It is shown in [71] that the Matérn covariance function is related to the Green’s function of a stochastic PDE. Namely, the Gaussian random field $m(\boldsymbol{x};\omega)$ with Matérn covariance (11) is the solution of the following linear stochastic PDE in $\mathbb{R}^{d}$, $d\in\mathbb{N}$, [45, 63, 21]

[TABLE]

where $\mathcal{W}\sim\mathcal{N}(0,\operatorname*{Id})$ is a spatial white noise in $\mathbb{R}^{d}$, and the normalization parameter $\eta$ is given by

[TABLE]

In practice, equation (14) is solved on a bounded domain with arbitrary boundary conditions [25, 31, 41].

In the general case, the shape operator can be spatially varying: $\boldsymbol{\Theta}=\boldsymbol{\Theta}(\boldsymbol{x})$ (see, for example, [63]), which leads to a non-stationary covariance. As mentioned above, we consider in this work the stationary covariance case and thus the spatially invariant shape operator $\boldsymbol{\Theta}$. This allows the solution of the equation (14) using the Fourier Transform:

[TABLE]

where $\mathcal{F}^{-1}$ denotes the operator of the inverse Fourier Transform, $\hat{\mathcal{W}}$ is the Fourier Transform of the white noise $\mathcal{W}$, and $\hat{G}$ is the Fourier Transform of the Green’s operator for the operator in (14), which is given by (see [4, Vol.II, section $8.13$, formula $(3)$])

[TABLE]

Use of the Fast Fourier Transform algorithm in convolution (16), along with periodic boundaries, makes possible a fast generation of the synthetic samples.

Thus, given design parameters $\phi_{0}$ (or $\tau$), $\nu$ and $\boldsymbol{\Theta}$, the process of generating a microstructure realization for the surrogate material model described above can be summarized as follows:

Draw a realization of the spatial white noise $\mathcal{W}(\boldsymbol{x})$. 2. 2.

Compute the corresponding intensity field $m(\boldsymbol{x})$ by solving (14). For the stationary case, it is given by convolution (16), computed using FFT. 3. 3.

The characteristic function $\chi(\boldsymbol{x})$ of the inclusions phase is given by the level-cut (3) of the intensity $m(\boldsymbol{x})$.

This surrogate material model is completely determined by only a few design parameters $\phi_{0}$, $\nu$ and $\boldsymbol{\Theta}$. However, by varying these parameters, the model can produce a large number of classes of random heterogeneous media (see examples in Figure 2). Besides, the parameters have clear intuitive meanings:

• •

The parameter $\phi_{0}$ is the expected volume fraction of the inclusions (or porosity).

• •

The regularity parameter $\nu$ defines the differentiability (smoothness) of the inclusions interface.

• •

The shape operator $\boldsymbol{\Theta}$ defines the metric to measure the distance in the covariance function and thus controls the form of the inclusions: their size, anisotropy level (aspect ratio) and orientation.

{rmk}

In [8], it is shown that, given the white noise is almost sure in $\mathcal{H}^{-d/2-\epsilon}(\mathbb{R}^{d})$, $\epsilon>0$, the solution of (14) is in $\mathcal{H}^{\nu-\epsilon}(\mathbb{R}^{d})$. Then, from the Sobolev embedding theorem, the intensity field $m(\boldsymbol{x})$ is continuous for $\nu>d/2$. And for $\nu\leq d/2$, we deal with a discontinuous intensity field $m(\boldsymbol{x})$, which produces a disjoint particles microstructure (see Figure 2(a)).

The regularity parameter $\nu$ is also related to the pore sphericity, one of the important parameters describing a microstructure [12, 42, 74]. We use the sphericity definition, proposed in [69], and define the relative pore size $r$:

[TABLE]

where $V_{pore}$ and $S_{pore}$ denote the volume of the pore and the surface of the pore interface respectively, and $V_{total}$ is the total volume of the sample. A correction factor of $0.833$ is taken into account due to the digitized structure of the surface [74]. Figure 3 shows the distribution of the relative pore size and the pore sphericity, computed for 3D surrogate samples of $2^{9\cdot 3}\approx 10^{8}$ voxels. We observe that with growing size the sphericity depends on $\nu$: the larger $\nu$ is, the more spherical are the pores. With the correlation length fixed, higher porosity $\phi_{0}$ provides larger maximum pore size.

We use the Matérn covariance class because of its flexibility, i.e. it covers a large variety of morphologies (microstructural shapes) with a relatively simple parametrization. Nevertheless, other convenient covariance kernel can be used instead (see, for example, [48, 65, 44, 67]). In the general case, in the absence of an explicit formula for the power spectrum, like (17), it can always be directly computed using FFT.

2.3 Bayesian inference of the model parameters

Let us suppose that several binary images (realizations $\omega$) $\chi(\boldsymbol{x};\omega)$ of a real material are given. Under the stationarity assumption, we can approximate the two-point probability function $S_{2}(r)$ of this material with the spatial average

[TABLE]

where $\boldsymbol{x}_{k}$ runs through all the voxels in all the provided images, $\chi(\boldsymbol{x})$ being zero outside the images (zero padding), and $N_{n}$ is the number of points $\boldsymbol{x}_{k}+\boldsymbol{r}_{n}$ belonging to the images. We compute the cross-correlation (19) using FFT with zero padding.

The two-point probability function $S_{2}(r)=\operatorname{\mathbb{E}}\left[\chi(\boldsymbol{x};\omega)\chi(\boldsymbol{x}+\boldsymbol{r};\omega)\right]$, $r=\lVert\boldsymbol{r}\rVert$, is the probability of the event when two points separated with distance $r$ are in the same phase. Then, the random variable $z_{k}(\boldsymbol{r};\omega)=\chi(\boldsymbol{x}_{k};\omega)\chi(\boldsymbol{x}_{k}+\boldsymbol{r};\omega)$ has the Bernoulli distribution with parameter $p=S_{2}(r)$. In (19) we compute the sum of ”trials” of $z(r;\omega,\boldsymbol{x})$, which is, thus, from the Binomial distribution. Since this tends with a large number of trials to the normal distribution $\mathcal{N}\left(p,\frac{p(1-p)}{N}\right)$ with the probability density function, we have

[TABLE]

Thus, the product of the densities for all $\boldsymbol{r}_{n}$ provides the likelihood for our Bayesian inverse problem:

[TABLE]

where we can substitute the real two-point correlation $S_{2}(r)$ within our surrogate model. The surrogate two-point probability function is given by formula (7):

[TABLE]

where $\mathcal{C}(r;\boldsymbol{\xi})$ is the intensity field covariance kernel, $\phi_{0}$ is the average volume fraction of the pores (porosity), $\tau$ is given by (9), and $\boldsymbol{\xi}$ denotes the vector of design parameters of the model. In our case, $\mathcal{C}(r;\boldsymbol{\xi})$ is of Matérn type in the form (11) with unit variance $\sigma=1$, and $\boldsymbol{\xi}=(\phi_{0},\nu,\boldsymbol{\Theta})$.

Assuming a uniform prior, the posterior distribution is simply determined by the likelihood. We use the Laplace approximation of the posterior: we approximate it with a multivariate lognormal distribution $\log\mathcal{N}(\log\boldsymbol{\xi}_{MLE},\boldsymbol{\Sigma})$, where the maximal likelihood estimator $\boldsymbol{\xi}_{MLE}$ is given by

[TABLE]

and the covariance matrix is approximated with

[TABLE]

In the isotropic case, we have $\boldsymbol{\Theta}=\ell^{2}\,\boldsymbol{\operatorname*{Id}}$. Then $\ell$ simply plays the role of the characteristic length scale and does not affect the microstructure morphology. So, we simply fix $\ell\equiv\ell_{MLE}$ as a deterministic parameter.

3 Fatigue Analysis: mechanical formulation

In what follows, we will perform a mechanical analysis on the surrogate microstructure with periodic boundary conditions under a given mean cyclic loading. The output of the computation will be analyzed in terms of quantities of interest for a standard high cycle fatigue analysis. The goal is to estimate the uncertainty of these quantities with respect to the variability of the underlying surrogate microstructures. From the point of view of a fatigue design method, the surrogate microstructure represents the material RVE and carries the statistical features of the real material as discussed in the previous sections. Moreover, the RVE represents a material point of a large structure and the given cyclic loading is the local shakedown response of the structure under service loading. In order to simplify the discussion and the subsequent computations, we will next only consider an elastic material behavior for the RVE. This assumption does however not restrain the generality of the present method and can be applied to materials with plasticity, viscosity, damage, etc., and can equally cover multiples material phases.

3.1 Linear elasticity problem

Let $D=(0,1)^{3}$ be a representative volume element (RVE) of a heterogeneous two-phase material. The phases are defined by a stochastic characteristic function $\chi(\boldsymbol{x};\omega)$ as in (1), periodic on $D$. The material is submitted to the macroscopic stress $\boldsymbol{\Sigma}$ to be imposed on average on the RVE, which results in a macroscopic strain $\boldsymbol{E}$. Thus, for all $\boldsymbol{x}\in\mathbb{R}^{d}$, the total displacement field $\boldsymbol{u}(\boldsymbol{x})$ is of the form [66]

[TABLE]

where the fluctuation $\tilde{\boldsymbol{u}}(\boldsymbol{x})$ is periodic on $D$.

Under the assumption of small strains, the local microscopic strain tensor field, observed on the stochastic RVE after loading, is

[TABLE]

and then condition (25) can be written in terms of strains:

[TABLE]

where $\tilde{\boldsymbol{\varepsilon}}(\boldsymbol{x})$ has zero mean and is periodic on $D$, and $\boldsymbol{E}$ denotes the macrostrain.

The local stress tensor field $\boldsymbol{\sigma}(\boldsymbol{x})$ is

[TABLE]

where $\boldsymbol{C}(\boldsymbol{x})$ is the fourth order stiffness tensor, $K(\boldsymbol{x})$ and $G(\boldsymbol{x})$ are the bulk and shear moduli, respectively, which are defined by the phase $\chi(\boldsymbol{x})$:

[TABLE]

where $K_{I},G_{I}$ are the bulk and shear moduli of the inclusions, and $K_{M},G_{M}$ – of the matrix, respectively. The equilibrium is given by

[TABLE]

where $\boldsymbol{\Sigma}$ denotes the macrostress.

Thus, for a given macrostress $\boldsymbol{\Sigma}$, we want to find the local strain field $\boldsymbol{\varepsilon}(\boldsymbol{x})=\boldsymbol{E}+\tilde{\boldsymbol{\varepsilon}}(\boldsymbol{x})$, such that the macrostrain $\boldsymbol{E}$ and the fluctuation $\tilde{\boldsymbol{\varepsilon}}(\boldsymbol{x})=\frac{1}{2}\left(\nabla\tilde{\boldsymbol{u}}(\boldsymbol{x})+\nabla\tilde{\boldsymbol{u}}(\boldsymbol{x})^{t}\right)$ are the solution of the system

[TABLE]

3.2 Mechanical Quantity of Interest for multidimensional loading

The high cycle fatigue analysis of metals is based on an elastic shakedown assumption at the structural scale and on material observations at the microscopic scale. Moreover, plasticity is driven by the deviatoric part of stresses, and cracks and porosities are opened by the positive spherical part of the stresses, which provides a natural split in the stress space. Therefore, a characteristic loading path can be characterized by a combination of bulk and shear loading as follows:

[TABLE]

where the angle $\theta$ varies from [math] to $\pi/2$ and characterizes the share of hydrostatic and deviatoric loading respectively. Lode coordinates of this loading are $z=\cos\theta$, $r=\sin\theta$ and $\theta_{s}=0$. In the stress space, the loading is a unit vector in the meridional plane $(z,r)$ and varies continuously from an isotropic traction with eigenstress $(\frac{1}{\sqrt{3}},\frac{1}{\sqrt{3}},\frac{1}{\sqrt{3}})$ obtained at $\theta=0$ to a pure shear with eigenstress $(-\frac{1}{\sqrt{2}},0,\frac{1}{\sqrt{2}})$ obtained at $\theta=\pi/2$.

We can compute for a given microstructure $\chi(\boldsymbol{x};\omega)$, $\boldsymbol{x}\in D$, the elastic solutions $\boldsymbol{\sigma}_{0}(\boldsymbol{x};\omega)$ and $\boldsymbol{\sigma}_{\pi/2}(\boldsymbol{x};\omega)$ corresponding to the loads $\boldsymbol{\Sigma}_{0}$ and $\boldsymbol{\Sigma}_{\pi/2}$, respectively. We decompose the local stress field $\boldsymbol{\sigma}(\boldsymbol{x})$ into a deviatoric part $\boldsymbol{\boldsymbol{s}}=\boldsymbol{\sigma}-\frac{1}{3}\operatorname{tr}(\boldsymbol{\sigma})\boldsymbol{\operatorname*{Id}}$ and a spherical part $\frac{1}{3}\operatorname{tr}(\boldsymbol{\sigma})\boldsymbol{\operatorname*{Id}}$, and introduce the local ”traction” $\operatorname{tr}_{+}(\boldsymbol{\sigma})=\max(\operatorname{tr}(\boldsymbol{\sigma}),0)$. Then, for each angle $\theta$, we can obtain the associated quantities

[TABLE]

and

[TABLE]

and deduce the associated combined quantity of interest over the volume $D$

[TABLE]

where

[TABLE]

is the microscopic damage parameter, and $\lVert\boldsymbol{s}\rVert=\sqrt{\boldsymbol{s}:\boldsymbol{s}}$. This weighted combination of deviatoric norm and of positive trace is often used in practical fatigue criteria [38]. The slope $0.3$ is an empirical value. The subdomain of integration $O_{t}(\omega,\theta)$ is here a ball of radius $\ell$ centered at the point $\boldsymbol{x}$ which maximizes $d(\boldsymbol{x};\omega,\theta)$, the pores being excluded.

We are interested in the expected value $\operatorname{\mathbb{E}}\left[Q(\theta)\right]$ and the variance $\operatorname{Var}\left[Q(\theta)\right]$, which are functions of $\theta$, and the ”worst” loading combination will be the one leading to the maximum expected value $\max\limits_{0\leq\theta\leq\frac{\pi}{2}}\operatorname{\mathbb{E}}\left[Q(\theta)\right]$.

The quantity $Q$ presents the mean squared damage parameter around the maximum point. Since, in general, the maximum point can be a singularity, the integral (47) of the damage parameter is preferred as quantity of interest. Having the estimation of the damage parameter in hand, one can proceed with the analysis of the number of loading cycles to failure [12, 38] or the failure probability [19].

4 Fatigue Analysis: results and discussions

In this section, we apply the surrogate microstructure model, discussed in Section 2, to the linear elasticity problem, formulated in Section 3, in order to perform a risk analysis of a heterogeneous material with respect to fatigue. Access to the fast microstructure sampling algorithm allows the approximation of the probability distribution of the fatigue criteria using Monte-Carlo methods. Furthermore, we can construct in this way a mapping from the microstructure design parameters to quantities of interest, which can serve for analysis of the influence of the microstructure morphology properties on the quantities of interest. Thus, in this section, we discuss two types of results:

• •

First, using the Laplace approximation for Bayesian inference of the microstructure design parameters (see Section 2.3), a surrogate material is reconstructed from a few CT images of a real material. The statistical properties of the simplified fatigue criteria (47) for the surrogate material are analyzed using the Monte-Carlo method, employing a fast Fourier-based sampling process and linear elasticity solver.

• •

Second, we perform the sensitivity analysis of the fatigue criteria and the homogenized elastic moduli with respect to the design parameters: porosity $\phi_{0}$ and pore regularity $\nu$. This makes possible the optimization of the microstructure design with respect to fatigue.

All simulations have been performed using an Intel Xeon E5 processor ($48$x $3.00$GHz) and $504$ GB RAM, computing Monte-Carlo samples in parallel. Each Monte-Carlo iteration includes:

Sample surrogate microstructure using Matérn covariance kernel (Section 2). 2. 2.

Solution of the linear elasticity problem (31) for the loadings $\theta=0$ and $\theta=\pi/2$ to obtain the local strains and stresses for the current microstructure sample (Section 3). Given the problem linearity, the results for intermediate loads $\theta\in(0,\pi/2)$ are computed as linear combination of the cases $\theta=0$ and $\theta=\pi/2$. 3. 3.

Computation of quantities of interest for the current sample (Subsection 3.2).

Runtime of one such iteration is about $1.5\,s$ for a 2D sample of $2^{8\cdot 2}\approx 6.5\cdot 10^{4}$ voxels, or about $70\,s$ for a 3D sample of $2^{7\cdot 3}\approx{2\cdot 10^{6}}$ voxels.

The linear elasticity problem (31) is solved in terms of local strains, using a Fourier-based Krylov solver, proposed by [9], which is a modification of the original method of [54], providing the convergence in case of the infinite contrast. See also [55, 50]. The problem is implemented with an in-house code written in Python and C++, using FFTW [27] and PETSc [3] libraries for the Fast Fourier Transform and the matrix-free conjugated gradient method respectively. SciPy [40] and Numdifftools [10] packages have been used to solve the optimization problem (23) and to compute the hessian (24) respectively.

4.1 Fatigue analysis of a real material

We reconstruct an aluminum–silicon alloy, A319-LFC, studied in [19], using six 2D CT images of the alloy (Figure 4, top). Using the Laplace approximation for Bayesian inference (see Section 2.3 for the procedure) of the design parameters of the surrogate material model introduced in Section 2, we obtain

[TABLE]

The estimator $\tau_{\operatorname{MLE}}$ corresponds to porosity ${\phi_{0}}_{\operatorname{MLE}}\approx 0.014$, $\nu$ varies in $(1.0,\,1.26)$. The correlation length is $\ell\approx 0.256\,mm$. With the design parameters in hand, we can generate the surrogate microstructure samples in 2D or in 3D. The original CT images and the surrogate 2D and 3D samples are presented in Figure 4. The two-point correlation function $S_{2}(r)$ obtained with MLE parameters is traced in Figure 5 (left) and compared with two-point correlations of original images (computed as cross-correlations using FFT). In [19], the experimentally computed distribution has been fitted with a lognormal and exponential distributions. In Figure 5 (right), we compare these reference distributions to the surrogate pore-size distribution, obtained from $10,000$ 3D samples of $2^{8\cdot 3}\approx{1.6\cdot 10^{7}}$ voxels. We observe that both distributions are in very good agreement.

For the surrogate material, we compute the quantity of interest $Q(\theta)$ for different loading types $\theta$, using $20,000$ 3D samples of $2^{7\cdot 3}\approx{2\cdot 10^{6}}$ voxels. The angle $\theta$ is discretized with $45$ equidistant points in $[0,\pi/2]$. The probability distributions of the quantity of interest $Q(\theta)$ are shown in Figure 6 (left). The ”worst” loading case $\theta_{\max}$, maximizing the average of $Q(\theta)$, is given by $\theta_{\max}\approx 0.3\,\pi$. The associated distribution is also depicted in Figure 6 (right), along with its lognormal fit.

4.2 Sensitivity to the design parameters

In what follows, $Q$ is associated with the load $\theta_{\max}=\max_{0\leq\theta\leq\frac{\pi}{2}}\operatorname{\mathbb{E}}\left[Q(\theta)\right]$, which corresponds to $0.3\,\pi<\theta_{\max}<0.32\,\pi$. We want to study the sensitivity of the quantity of interest $Q$ (47) to the two design parameters: porosity $\phi_{0}$ and pore regularity $\nu$. To be representative, a volume element has to be of size much larger than the correlation length $\ell$. So we consider a unit square as RVE and fix $\ell=0.05$. For each $\phi_{0}$ and $\nu$, we solve the plane-strain problem (31) on $10,000$ 2D samples with $2^{8\cdot 2}\approx 6.5\cdot 10^{4}$ voxels. Young’s modulus is $1\,GPa$ for the material matrix and [math] in the pores, material Poisson ratio is $0.3$. Solution examples for $\nu=0.5$ and $\nu=10$ with porosity $\phi_{0}=0.1$ are shown in Figure 7. Figure 8 depicts the average of $Q$ as function of porosity $\phi_{0}$ and of pore regularity $\nu$. The standard deviation is denoted with the error bars. We observe that the damage parameter grows with the porosity and with the pore regularity. We note the asymptotic behavior of $Q$ when $\nu$ goes to infinity. This limit case corresponds to infinitely differentiable pore interfaces, when the Matérn covariance of the intensity field becomes a squared exponential covariance. We also note that the influence of the parameter $\nu$ increases with the porosity. The standard deviation of $Q$ also grows with the porosity.

Finally, in addition, we want to estimate the homogenized elastic moduli, $K$ (bulk) and $G$ (shear). Figure 9 shows their averages as functions of porosity, compared to the corresponding upper Hashin-Shtrikman bounds [35]. The standard deviation is presented there with error bars. As before, we observe the standard deviation growing with the porosity.

Figure 10 shows the scattering of the homogenized moduli over $100$ samples for each value of $\nu$ and fixed $\phi_{0}=0.2$. We remark that the parameter $\phi_{0}$ is the average volume fraction, and thus the computed average $K$ and $G$ respect the bound. However, a particular realization has random volume fraction and can, therefore, violate the bound (see Figure 10). For each sample, the associated homogenized moduli are approximated as

[TABLE]

where $\boldsymbol{\Sigma}_{\theta}$ is the imposed macro stress, and $\boldsymbol{E}_{\theta}$ is the associated computed macro strain (see Subsection 3.2).

5 Conclusion

In this work, we have demonstrated a framework for uncertainty quantification in fatigue analysis using a surrogate microstructure model. The proposed microstructure model presents a trade off between computational performance, realistic random microstructure reconstruction and the small number of design parameters.

The approach presented here is quite general, although for specifity we have focused on the simple case of statistically isotropic media. The anisotropy can be included through the shape operator $\boldsymbol{\Theta}$ in (11). Modifying the ratio of different correlation lengths, we change the pore aspect ratio. In more sophisticated cases, when $\boldsymbol{\Theta}(\boldsymbol{x})$ depends on the position, the sampling procedure can be still based on the solution of the SPDE (14), however, the FFT can not be employed any more.

In this work, we have considered covariances of Matérn class [71, 48, 34], which presents a large class of covariance functions and provides enough flexibility in material reconstruction. In addition to the traditional morphology descriptors as porosity, size and aspect ratio, this covariance class provides the regularity of the inclusion interface, which is also related to the pore sphericity [12, 42]. In order to introduce more design parameters and thus extend the class of possible reconstructed materials, one can consider other more general covariance kernels [48, 65, 44, 67].

We also demonstrated the use of the model in the fatigue analysis, and we have studied the influence of the regularity of the pore interface on the statistical properties of the simplified fatigue criteria and of the homogenized elastic moduli.

Appendix A Proof of Lemma 2.1

Proof A.1**.**

We show the formulas (5)-(7) for $S_{1}$ and $S_{2}$ in two steps. First, we consider $S_{1}$. Since the intensity $m(\boldsymbol{x};\omega)$ is a Gaussian random field, the mean of the phase $\chi(\boldsymbol{x};\omega)$, defined by (3), satisfies

[TABLE]

Let us consider the volume fraction of all inclusions

[TABLE]

Then, its average over samples $\phi_{0}=\operatorname{\mathbb{E}}\left[\varphi\right]$ is

[TABLE]

The mean volume fraction $\phi_{0}$ can be also directly related to $\tau$ by the Gauss error function $\operatorname{erf}(\cdot)$, i.e.,

[TABLE]

Now let us focus on $S_{2}$. We denote the bivariate Gaussian covariance matrix by

[TABLE]

where $\mathcal{C}(\boldsymbol{x},\boldsymbol{y})$ is the covariance function of the Gaussian field $m(\boldsymbol{x};\omega)$. Then,

[TABLE]

Hence, the two-point correlation function of the phase $\chi$ is given by

[TABLE]

where $a=\frac{1}{2\pi\sqrt{1-g^{2}}}$. Let us denote with $S_{2}^{-}$ and $S_{2}^{+}$ the integrals

[TABLE]

Thus, $S_{2}=S_{2}^{-}+S_{2}^{+}$. Let us first consider the integral $S_{2}^{-}$. After a change of variable $\hat{\xi}_{1}=\frac{\xi_{1}-g\xi_{2}}{\sqrt{1-g^{2}}}$, it becomes

[TABLE]

which can be written as

[TABLE]

where

[TABLE]

is Owen’s T function [58, 60]. Similarly, changing the variable $\hat{\xi}_{1}=\frac{\xi_{1}+g\xi_{2}}{\sqrt{1-g^{2}}}$ in $S_{2}^{+}$ leads to

[TABLE]

Thus, we have

[TABLE]

Moreover, from

[TABLE]

we have that $S_{2}(\boldsymbol{x},\boldsymbol{x})=\phi_{0}$ and $\lim_{\lVert\boldsymbol{x}-\boldsymbol{y}\rVert\rightarrow\infty}S_{2}(\boldsymbol{x},\boldsymbol{y})=\phi_{0}^{2}$.

Alternatively, given $T(\tau/\sigma,1)=\frac{1}{4}\phi_{0}(1-\frac{1}{2}\phi_{0})$, changing the variable $z=\sqrt{\frac{1-t}{1+t}}$ in (69) and $z=\sqrt{\frac{1+t}{1-t}}$ in (73) leads respectively to

[TABLE]

and

[TABLE]

Hence,

[TABLE]

*Similar but more general formula can be found in [6, eq.(33)]. *