## 1 Introduction

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

(1a) | |||||

(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 partial differential equations (PDEs) 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) |

Here denotes the standard inner product in and .

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

(3) |

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

(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

(5a) | |||||

(5b) |

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

(6) |

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

(7a) | ||||

(7b) | ||||

(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

(8a) | ||||

(8b) |

Consequently, the deterministic homogenized coefficient matrix

can be practically approximated by a random matrix

whose entry is defined as(9) |

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

(10a) | |||||

(10b) |

## 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

(11) |

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

(12a) | ||||

(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

(13a) | |||||

(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.

### 3.2 Characterization of the equivalent coefficient matrix

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:

(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

(15) | ||||

where the cell function satisfies equationparentequation

(16a) | ||||

(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

(17) |

and

(18) | ||||

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

(19) |

Introduce the self-adjoint covariant operator as

(20) |

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

By the Karhunen-Loéve expansion, we obtain

(21) |

where and is a standard normal variable given by

(22) |

The principal eigenvalue

satisfies(23) |

which implies that

(24) |

Thus, -component of satisfies

(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

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 :

(26) |

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

(27a) | ||||

(27b) | ||||

(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

###### Proof

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

(31) |

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

(32) |

Then we have

###### Theorem 3

###### Proof

###### 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

(39) |

for some independent of and .

###### Proof

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

(42) |

for some positive constant independent of .