DeepAI

# An efficient multi-modes Monte Carlo homogenization method for random materials

In this paper, we propose and analyze a new stochastic homogenization method for diffusion equations with random and fast oscillatory coefficients. In the proposed method, the homogenized solutions are sought through a two-stage procedure. In the first stage, the original oscillatory diffusion equation is approximated, for each fixed random sample w, by a spatially homogenized diffusion equation with piecewise constant coefficients, resulting a random diffusion equation. In the second stage, the resulted random diffusion equation is approximated and computed by using an efficient multi-modes Monte Carlo method which only requires to solve a diffusion equation with a constant diffusion coefficient and a random right-hand side. The main advantage of the proposed method is that it separates the computational difficulty caused by the spatial fast oscillation of the solution and that caused by the randomness of the solution, so they can be overcome separately using different strategies. The convergence of the solution of the spatially homogenized equation (from the first stage) to the solution of the original random diffusion equation is established and the optimal rate of convergence is also obtained for the proposed multi-modes Monte Carlo method. Numerical experiments on some benchmark test problems for random composite materials are also presented to gauge the efficiency and accuracy of the proposed two-stage stochastic homogenization method.

• 2 publications
• 5 publications
• 19 publications
• 2 publications
08/12/2021

### Multilevel Monte Carlo estimators for elliptic PDEs with Lévy-type diffusion coefficient

General elliptic equations with spatially discontinuous diffusion coeffi...
03/30/2020

### Multiscale stochastic reduced-order model for uncertainty propagation using Fokker-Planck equation with microstructure evolution applications

Uncertainty involved in computational materials modeling needs to be qua...
07/03/2021

### An asymptotically compatible probabilistic collocation method for randomly heterogeneous nonlocal problems

In this paper we present an asymptotically compatible meshfree method fo...
10/24/2022

### An efficient Monte Carlo scheme for Zakai equations

In this paper we develop a numerical method for efficiently approximatin...
10/26/2019

### The Vlasov-Fokker-Planck Equation with High Dimensional Parametric Forcing Term

We consider the Vlasov-Fokker-Planck equation with random electric field...
06/25/2022

### A two-stage method for reconstruction of parameters in diffusion equations

Parameter reconstruction for diffusion equations has a wide range of app...
09/21/2019

### Multithreaded Filtering Preconditioner for Diffusion Equation on Structured Grid

A parallel and nested version of a frequency filtering preconditioner is...

## 1 Introduction

This paper is concerned with numerical solutions of the following diffusion equation with random coefficients and data encountered in materials science: equationparentequation

 =f(x,ω) in D×Ω, (1a) uε(x,ω) =0 on ∂D×Ω. (1b)

Here is a bounded domain and

denotes a sample point which belongs to a probability (sample) space

. The coefficient matrix and the right hand side term are random fields with continuous and bounded covariance functions. The parameter represents the size of microstructure for the composite materials, which is usually very small, that is, .

The random diffusion equation (1) has many applications in mechanics, hydrology and thermics (see [15, 19, 24, 26]). A direct accurate numerical solution of (1) is difficult to obtain because it requires a very fine mesh and large-scale sampling of , and thus a prohibitive amount of computation time. In the case of absence of the randomness (i.e., the dependence on in (1) is dropped ), the homogenization method has been successfully developed for solving the diffusion equation with periodic deterministic coefficients (cf. [5, 8, 25]), in which the homogenized coefficients are obtained by solving a cell problem defined in the unit cell. For the diffusion equation (1) with random coefficients, a stochastic homogenization theory has also been well developed, see [6, 7, 10, 9, 14, 18, 17, 21, 23, 27, 30]. Similar to the deterministic case, the homogenized diffusion equation is constructed by solving a certain cell problem. However, a fundamental difference is that the stochastic cell problem is a random second order elliptic problem, which is posed in the whole space (see (7)). Solving such an infinite domain problem numerically is not only challenging but also very expensive, and the homogenization methods quoted above did not give any practical recipe for numerically approximating the cell problem and the homogenized equation.

To circumvent the above difficulties, some localized approximations of the effective (or homogenized) coefficients by using “periodization” and “cut-off” procedures were introduced in [2, 7, 22, 31]. A big benefit of the localized approximations is that the resulting cell problem is now posed on a bounded domain and it was proved that the approximated coefficients converge to the effective (or homogenized) coefficients as the size of the bounded domain goes to infinity. We also note that the localized approximation methods, such as the representative volume element (RVE) method, have also been used to compute the effective parameters of highly heterogeneous materials and to calculate the effective coefficients related to random composite materials, by utilizing a possibly large number of realizations [9, 13, 18, 29]. After having constructed the approximated effective (or homogenized) coefficients, the main task then reduces to solve the approximated (random) diffusion equation with the constructed coefficients. As a direct application of Monte Carlo or stochastic Galerkin method, since solving this random diffusion equation is computationally expensive, other more efficient methods have been developed for the job. In [3, 4], a perturbative model for weakly random materials was proposed and the first-order and second-order asymptotic expansions were established by means of an ergodic approximation based on the weak randomness assumption. In [11], a multi-modes Monte Carlo (MMC) finite element method was proposed to solve random under the assumption that the medias (or coefficients) are weakly random in the sense that they can be expressed as small random perturbations of a deterministic background. However, to the best of our knowledge, there is still no efficient method for solving the homogenized equation for problem (1) in general setting.

The goal of this paper is to develop and analyze an efficient and practical two-stage homogenization method for problem (1). In the first stage, we construct a piecewise homogeneous (i.e., piecewise constant) material as an approximation to the original composite material for each fixed sample by solving several cell functions on bounded domains. Under the stationary (process) assumption, we are able to prove that the coefficient matrix of the piecewise homogeneous material can be rewritten as a small random perturbation of a deterministic matrix, which then sets the stage for us to adapt the MMC framework. In the second stage, we utilize the MMC finite element method to solve random diffusion problem with the piecewise coefficients obtained from the first stage and provide a complete convergence analysis for the MMC method. Since the first stage of the proposed method is similar to the RVE method, the work of this paper can be regarded as a mathematical interpretation and justification for the RVE method for problem (1) and introduces a numerical framework for an efficient implementation of the RVE method.

The rest of the paper is organized as follows. In Section 2, we present a few preliminaries including notations and assumptions. In Section 3, we introduce our two-stage stochastic homogenization method and characterize the piecewise constant approximate coefficients. In Section 4, we present the convergence analysis of the proposed two-stage stochastic homogenization method under the stationary assumption on the diffusion coefficients. In Section 5, we propose a finite element discretization and a detailed implementation algorithm for the proposed method. In Section 6, we present several benchmark numerical experiments to demonstrate the efficiency of the proposed method and to validate the theoretical results. Finally, the paper is completed with some concluding remarks given in Section 7.

## 2 Preliminaries

### 2.1 Notations and assumptions

Standard notations will be adopted in paper. denotes a probability space and

stands for the expectation value of random variable

. Let be the unit cell and for with . Let be a bounded domain which can be written as , where . For a positive integer , set and . Let denote the set of invertible real-valued matrices with entries in and satisfying -a.s

 α|ξ|2≤(Aξ,ξ)≤β|ξ|2for any ξ∈Rd and a.e. in D. (2)

Here denotes the standard inner product in and .

Similar to [2, 3], we also assume that is stationary in the sense that

 A(xε+k,ω)=A(xε,τkω)for any k∈Zd, a.e. in D and P-a.s., (3)

where is a mapping which is ergodic and preserves the measure , that is

 τkE=E∀E∈F% implies thatP(E)=0 or 1. (4)

For the ease of presentation, we set in the rest of the paper, that is, is a deterministic function.

### 2.2 Elements of the classical stochastic homogenization theory

It is well known that [2, 6, 7, 17, 21, 23] as , the solution of equation (1) converges to the solution of the following homogenized problem: equationparentequation

 −div(A⋆∇u⋆(x)) =f(x) in D, (5a) u⋆(x) =0 on ∂D, (5b)

where the entry of the homogenized matrix (or effective coefficient) is defined by

 a⋆ij=E(∫Q(ei+∇Nei(y,ω))TA(y,ω)ejdy). (6)

denotes the canonical basis of and the cell function is defined as the solution of the following cell problem: equationparentequation

 −div[A(y,ω)(ei+∇Nei(y,ω))] =0in Rd, (7a) E(∫Q∇Nei(y,⋅)dy) =0, (7b) ∇Nei(y,ω) is~{}stationary~% {}in ~{}the~{}sense~{}of~{}(???). (7c)

As shown above, the classical stochastic homogenization method obtains the homogenized solution in one step, see Fig. 1 for a schematic explanation. We note that the cell problem (7) is random elliptic problem which is posed on the whole space and can not be reduced to a cell problem on a bounded domain due to the global constraint . Solving problem (7) is the main computational challenge for implementing the classical stochastic homogenization method. A natural and widely used approach (cf. [2, 14]) is to approximate by a truncated cubic domain with size by using “periodization” and “cut-off” techniques and then to solve the truncated problem equationparentequation

 −div[A(y,ω)(ei+∇Nei,N(y,ω))]=0 in QN, (8a) Nei,N(y,ω) is QN-% periodic. (8b)

Consequently, the deterministic homogenized coefficient matrix

can be practically approximated by a random matrix

whose entry is defined as

 a⋆ij,N(ω)=1|QN|(∫QN(ei+∇Nei,N(y,ω))TA(y,ω)(ej+∇Nej,N(y,ω))dy). (9)

Then, the solution of problem (5) is approximated as with being the solution of the following equation equationparentequation

 −div(A⋆N(ω)∇u⋆N(x,ω)) =f(x) in D, (10a) u⋆N(x,ω) =0 on ∂D. (10b)

We refer the reader to [2, 4, 14] for a detailed account about the above classical numerical approach.

## 3 Two-stage stochastic homogenization method

In this section, we shall present a detailed formulation of our two-stage stochastic homogenization method for problem (1). The new method only needs to solve similar diffusion equations with constant diffusion coefficients and random right-hand sides.

### 3.1 Formulation of the two-stage stochastic homogenization method

As explained earlier, the main difficulty for solving problem (1) is due to the oscillatory nature of its solution which is caused by the oscillatory coefficient matrix of the problem. Recall that the classical numerical homogenization methods approximate effective (or homogenized) coefficient matrix by matrix which is formed by solving the cell problem (8), which is often expensive to solve numerically. Motivated by this difficulty, the main idea of our method is to propose a different procedure to construct an approximation to (in the first stage) whose corresponding homogenized problem can be solved efficiently (in the second stage).

Specifically, our proposed method aims to construct a homogenized solution by the following two stages as illustrated in Fig. 1. In the first stage, for each given sample , the composite material with micro-structure is equivalently transformed to a piecewise homogeneous material with coefficient matrix (see Fig. 2), referred to as the equivalent matrix in each block , where

 ^akij(ω)=1|QM|∫QM(ei+∇Nkei(y,ω))TA(y+Mk,ω)(ej+∇Nkej(y,ω))dy, (11)

and the cell function is defined as the solution of the following cell problem on block : equationparentequation

 −div[A(y,ω)(ei+∇Nkei(y,ω))]=0 ~{}in QkM, (12a) Nkei(y,ω) is QkM-periodic. (12b)

Here is a parameter used to balance the efficiency and accuracy of the proposed method. In numerical simulations, we usually choose . Notice that the equivalent matrix in each block is a constant matrix. The equivalent matrix can be regarded as a coarse-grained approximation of the original matrix . In the coarsening process, the equivalent material with coefficient matrix is homogeneous in each block , but still maintain the heterogeneity between different blocks (see Fig. 2-(b)). It should be pointed out that the equivalent matrix is usually different from obtained by “periodization” procedure. In fact, by comparing cell problems (12) and (8) with , it is easy to see that is the same as and may be different from for due to the possible heterogeneity between different blocks (recall that denotes the equivalent matrix in block ).

In the second stage, we intend to solve the random diffusion problem with the equivalent (piecewise constant) coefficient matrix , namely, equationparentequation

 −div(^A(x,ω)∇^u(x,ω)) =f(x) in D, (13a) ^u(x,ω) =0 on ∂D. (13b)

In other words, the original oscillatory random coefficient matrix is approximated by the equivalent matrix which is a constant matrix in each block . For a given , the computational cost for solving the homogenized problem (13) is less than that for solving the original problem (1). However, the equivalent matrix fluctuates on different blocks due to the non-periodicity. The fluctuation leads to expensive computational costs for solving the homogenized problem (13) with small parameter and , because the computational mesh size must be proportional to . To overcome the difficulty, we adapt the MMC finite element method of [11] to solve (13) in an efficient way. This is possible thanks to our discovery which shows that the equivalent matrix has a nice structure, that is, it can be rewritten as a small random perturbation of the deterministic matrix , see Section 3.2. This then sets an ideal stage for us to solve problem (13) by using the MMC method. The leading term in the MMC approximation will be defined as , which is the after-sought approximate solution alluded earlier, see Section 3.3.

It is important to point out that the MMC method presented in [11] can not be directly applied to a random diffusion problem (1) because its diffusion coefficient does not satisfy the weak media assumption of the MMC method (see Appendix A).

### 3.2 Characterization of the equivalent coefficient matrix ^A(x,ω)

In this subsection, we show that the equivalent matrix can be rewritten as a (small) random perturbation of a deterministic matrix. The precise statement is given in the following theorem.

###### Theorem 1

Suppose satisfies stationary hypothesis (3). Then the equivalent matrix can be rewritten as:

 ^A(x,ω)=E(^A0(ω))+δA1(x,ω), (14)

where denotes the equivalent matrix in any block , with , and is a (small) parameter which depends on and .

###### Proof

By the stationary assumption (3), the -entry of can be written as

 ^akij(ω)=1|QM|∫QM(ei+∇Nkei(y,ω))TA(y+Mk,ω)(ej+∇Nkej(y,ω))\,dy, (15) =1|QM|∫QM(ei+∇Nkei(y,τMkω))TA(y,τMkω)(ej+∇Nkej(y,τMkω))\,dy,

where the cell function satisfies equationparentequation

 −div[A(y,τMkω)(ei+∇Nkei(y,τMkω))]=0 ~{}in QM, (16a) Nkei(y,τMkω) is QM-periodic. (16b)

Thus, , which shows that the equivalent matrix coincides with on each block and at each sample. Since the ergodic mapping preserves the measure , then we have

 E(^Ak(ω))=∫Ω^Ak(ω)dP(ω)=∫Ω^A0(τMkω)dP(τMkω)=E(^A0(ω)), (17)

and

 Var(^Ak(ω)) =∫Ω(^Ak(ω)−E(^Ak(ω)))2dP(ω), (18) =∫Ω(^A0(τMkω)−E(^A0(τMkω)))2d% P(τMkω)=Var(^A0(ω)).

Next, we derive a (small) random perturbation form for the equivalent matrix . For a given , and the autocorrelation function of , which is a constant, is defined by

 Covk=Var(^Ak(ω))=Var(^A0(ω)). (19)

Introduce the self-adjoint covariant operator as

 Tkv(⋅):=∫εQkMCovkv(x)dx=Covk∫εQkMv(x)dx∀v∈L2(εQkM). (20)

Let denote a complete eigen-set of the operator with and

 ∫εQkMφl(x)φm(x)d% x=δlml,m=1,2,⋯.

By the Karhunen-Loéve expansion, we obtain

 ^akij(ω)=E(^akij(ω))+∞∑l=1√λlZkl(ω)φl(x)=E(^akij(ω))+√λ1Zk1(ω)φ1(x), (21)

where and is a standard normal variable given by

 Zk1(ω)=^akij(ω)−E(^akij(ω))√λ1∫εQ1φ1(x)dx=√|εQk|^akij(ω)−E(^akij(ω))√λ1. (22)

The principal eigenvalue

satisfies

 λ1 =∫εQkMTkφ1⋅φ1(x)dx∫εQkMφ1(x)φ1(x)dx=Covk(∫εQkMφ1(x)dx)2=∣∣εQkM∣∣Covk, (23)

which implies that

 ^akij(ω)=E(^akij(ω))+√CovkZk1(ω)=E(^a0ij(ω))+√Var(^a0ij(ω))Zk1(ω). (24)

Thus, -component of satisfies

 ^aij(x,ω)=E(^a0ij(ω))+√Var(^a0ij(ω))∑kZk1(ω)Ξk(x). (25)

where

stands for the characteristic function of the domain

.

Finally, setting , then the equivalent matrix can be rewritten as a (small) random perturbation as stated in (14).

###### Remark 1

The smallness of and the precise relationship between and as well as will be given in Theorem 3 later in the next section.

### 3.3 An efficient multi-modes Monte Carlo method for solving the homogenized problem (13)

By Theorem 1 we know that the equivalent matrix can be rewritten as a (small) random perturbation of as given in (14). This then sets the stage for us to solve homogenized problem (13) by using the MMC method. To proceed, we first notice that with satisfying

 P{ω∈Ω;∥∥a1ij(ω)∥∥L∞(D)≤¯a}=1.

In [11], the perturbation term is assumed to be in . However, the perturbation term in (14) is not in because is a piecewise constant function which is discontinuous in . Nevertheless, we show below that the MMC method can be easily extended to the case.

Due to the linear nature of the equivalent problem Eq. 13 and the small random perturbation structure of the equivalent matrix , we can postulate the following multi-modes expansion for :

 ^u(x,ω)=∞∑n=0δnu0n(x,ω). (26)

Substituting Eq. 14 and Eq. 26 into Eq. 13 and matching the coefficients of order terms for we get equationparentequation

 −∇⋅(E[^A0(ω)]∇u00(x,ω)) =f(x), (27a) −∇⋅(E[^A0(ω)]∇u0n(x,ω)) =∇⋅(A1(x,ω)∇u0n−1(x,ω))∀n≥1, (27b) u0n(x,ω) =0on∂D∀n≥0. (27c)

Clearly, the first mode function satisfies a diffusion equation with a deterministic coefficient matrix and a deterministic source term . Thus, is independent of and we relabel it as . Moreover, the mode functions satisfy a family of diffusion equations that have the same deterministic diffusion operator but different right-hand side source terms. Furthermore, are defined recursively with the current mode function being only dependent directly on the proceeding mode function . The well-posedness of multi-mode functions and the corresponding error estimates will be discussed in the next section (see Theorem 4 and Theorem 5).

In this paper, we shall use the first term of the multi-modes expansion, that is, as an approximation for . How to efficiently compute the mode functions is important and challenging, since in the right-hand side of (27b) fluctuates in probability space and oscillates in different blocks . A natural but expensive approach is to solve (27b) for each sample by using a fine mesh with mesh size proportional to , which is easily implemented but low efficient. According to Theorem 5 in the next section, taking as an approximation of only enjoys the first order convergence rate. The convergence rate can be further improved by using more mode functions . Thus, developing more efficient numerical methods and algorithms for computing the mode functions are important and will be addressed in future work.

## 4 Convergence analysis

In this section, we analyze the convergence and error estimates for the proposed two-stage stochastic homogenization method. We first show the convergence of the equivalent matrix to the homogenization matrix as in the next theorem.

###### Theorem 2

Suppose satisfies the stationary hypothesis (3). Let be the equivalent matrix defined by (11) and be defined by (6), then there holds

 limM→∞∥^aij(⋅,ω)−a⋆ij∥L∞(D)=0a.s.. (28)

###### Proof

By Theorem 1 in [7], we have

 limM→∞∥^a0ij(⋅,ω)−a⋆ij∥L∞(D)=0a.s.. (29)

For any , it follows that

 (30)

The proof is complete.

To derive the rate of convergence of , we introduce the uniform mixing condition as in [7]. For a given random field in , let denote the -algebra . The uniform mixing coefficient of is defined as

 γ(s)=sup~D1,~D2⊂Rd, dist(~D1,~D2)≥ssup~D1∈F~D1,~D2∈F~D2∣∣P(~D1∩~D2)−P(~D1)P(~D2)∣∣. (31)

In the rest of this section, we assume satisfies the following growth condition:

 γ(s)≤c(1+s)−θfor~{}some θ>0 and ∀s>0. (32)

Then we have

###### Theorem 3

Suppose satisfies the stationary hypothesis (3). Assume that the uniform mixing coefficient of satisfies (32). Let and , then there holds

 E[(^aij(x,⋅)−a⋆ij)2]≤CM−ζfor a.e. x∈D, (33)

and there exists such that equationparentequation

 Var(^a0ij(ω)) ≤CM−ζ, (34a) δ ≤CM−ζ/2. (34b)

###### Proof

By Theorem 5 of [7], we have

 E[(^a0ij(⋅)−a⋆ij)2]≤CM−ζfor a.e. x∈D. (35)

For any , similar to (18), we have

 E[(^akij(⋅)−a⋆ij)2]=E[(^a0ij(⋅)−a⋆ij)2]≤CM−ζfor a.e. x∈D, (36)

which completes the proof of (33). To show (34), using the Hölder’s inequality and (33), we get

 (E[^a0ij]−a⋆ij)2=(E[^a0ij−a⋆ij])2≤E[(^a0ij−a⋆ij)2]≤CM−ζ, (37)

which implies

 Var(^a0ij) =E[(^a0ij−E[^a0ij])2], (38) ≤CM−ζ.

Thus, the proof is complete.

###### Remark 2

One immediate corollary of Theorem 3 is that as after taking with . This verifies that is indeed a small parameter.

Following Theorem 3.1 of [11], we can show the unique existence and the stability estimate for each mode in (26).

###### Theorem 4

Assume that and . There exists a unique solution to problem (27) for each which satisfies

 E(∥u0n∥2H1(D))≤Cn+10∥f∥2L2(D). (39)

for some independent of and .

###### Proof

Using Hashin-Shtrikman bounds [16], we have . For , the existence of a unique weak solution follows immediately from the Lax-Milgram Theorem (see, e.g. [12]) and there holds

 E(∥u00∥2H1(D))≤C10∥f∥2L2(D). (40)

The proof for the cases can be done by the induction argument. Assume that (39) holds for , then (which is the source term in (27)). Using Theorem 3.3 in [20], there exists a unique which solves (27) for and fulfills the following estimate:

 E(∥u0l∥2H1(D)) ≤C(D)E(∥A1(x,ω)∇u0n−1∥2[L2(D)]d), (41) ≤dC(D)¯a2E(∥∇u0n−1∥2[L2(D)]d), ≤dC(D)¯a2E(∥u0n−1∥2H1(D)), ≤Cl+10∥f∥2L2(D).

The proof is complete.

The next theorem establishes an estimate for the error function , which is an analogue to Theorem 3.2 of [11].

###### Theorem 5

Suppose that satisfies stationary hypothesis (3). Assume that and , we have

 E(∥∥^u(x,ω)−u00(x)∥∥2H1(D))⩽C1δ2∥f∥2L2(D). (42)

for some positive constant independent of .

###### Proof

Let . Substituting (14) into (13) and combining it with the first equation of Eq. 27, we obtain equationparentequation

 −div(E(^A(x,ω))∇r(x,ω)) =δdiv(A1(x,ω)∇^u(x,ω)) in D, (43a) r(x,ω) =0 on ∂D. (43b)

By Theorem 3.3 of [20], we get

 E(∥r∥2H1(D)) ⩽C(D)δ2E(∥A1(x,ω)∇^u