# A Completely Discrete Heterogeneous Multiscale Finite Element Method for Multiscale Richards’ Equation of van Genuchten-Mualem Model

Authors: Haitao Cao; Tao Yu; Xingye Yue
Journal: Journal of Applied Mathematics
 (2014)

---

## Abstract

We propose a fully discrete method for the multiscale Richards’ equation of van Genuchten-Mualem model which describes the flow transport in unsaturated heterogenous porous media. Under the framework of heterogeneous multiscale method (HMM), a fully discrete scheme combined with a regularized procedure is proposed. Including the numerical integration, the discretization is given byC
0 piecewise finite element in space and an implicit scheme in time. Error estimates between the numerical solution and the solution of homogenized problem are derived under the assumption that the permeability is periodic. Numerical experiments with periodic and random permeability are carried out for the van Genuchten-Mualem model of Richards’ equation to show the efficiency and accuracy of the proposed method.


---

## Body

## 1. Physical Model

Richards’ equation is most often used to model the movement of groundwater flow in saturated-unsaturated porous media. It was formulated by Richards in 1931 [1]. It is a nonlinear partial differential parabolic equation. Depending on the saturation and pressure, three main forms of Richards’ equation are usually presented: pressure-based form, saturation-based form, or mixed form. Combining the continuity equation with Darcy’s law, the mixed form can be expressed as
(1)
∂
t
Θ
(
p
)
-
∇
(
K
(
x
ϵ
,
Θ
(
p
)
)
∇
(
p
+
z
)
)
=
0
,
where the reduced saturation is defined by Θ
=
(
θ
-
θ
r
)
/
(
θ
s
-
θ
r
)
∈
[
0,1
]. θ
r (the residual fluid content) and θ
s (the fluid content at saturation) are constants relying on the porous medium. θ denotes the saturation and p denotes the pressure. Due to the complex heterogeneity of natural medium, the permeability K
(
x
/
ϵ
,
Θ
(
p
)
) of the medium oscillates rapidly with large contrast. ϵ
(
>
0
) is the characteristic length presenting the small scale variability of the media. The coordinate in the direction of gravity is denoted by z.In this paper, we adopt the retention curves proposed by van Genuchten [2] and Mualem [3]. The relations of the reduced saturation Θ, pressure p, and permeability K are read as
(2)
Θ
(
p
)
=
{
1
(
1
+
|
c
p
|
n
)
m
for
p
≤
0
,
1
for
p
>
0
,
K
(
x
ϵ
,
Θ
)
=
K
s
(
x
ϵ
)
K
r
(
Θ
)
=
K
s
(
x
ϵ
)
·
Θ
1
/
2
[
1
-
(
1
-
Θ
1
/
m
)
m
]
2
,
where c, m
∈
(
0,1
) and n
=
1
/
(
1
-
m
) are the parameters of porous medium. K
s and K
r can be presented as the absolute and relative permeability, respectively. When p
>
0, Θ
=
1, the medium is in the saturated state; otherwise, it is in the unsaturated case.Under the unsaturated case, Richards’ equation can also be rewritten in terms of reduced saturation by using the above relations. Consider(3)
∂
t
Θ
-
∇
(
D
(
x
ϵ
,
Θ
)
∇
Θ
+
K
z
→
(
x
ϵ
,
Θ
)
)
=
0
,
with the moisture diffusivity D defined by
(4)
D
(
x
ϵ
,
Θ
)
=
-
K
(
x
ϵ
,
Θ
)
∂
p
∂
Θ
=
(
1
-
m
)
K
s
(
x
/
ϵ
)
α
m
(
θ
s
-
θ
r
)
Θ
(
1
/
2
)
-
(
1
/
m
)
×
[
(
1
-
Θ
1
/
m
)
m
+
(
1
-
Θ
1
/
m
)
-
m
-
2
]
≐
D
s
(
x
ϵ
)
D
r
(
Θ
)
.
Here, under the relations (2), we should point out that (1) can model both saturated and unsaturated cases. In saturated region, (1) is an elliptic equation; in unsaturated region, the equation is parabolic. On the other hand, in unsaturated case, (1) can be degenerated because K
r may approximate to zero when Θ
→
0 (in almost completely dry region). Richards’ equation (3) can only model the unsaturated case. But (3) may also be degenerated, because the diffusivity D may vanish or explode (see Figure 1(a)). For Θ
=
0 (in almost completely dry region), the diffusivity D
r vanishes, while for Θ going to 1 (in almost full saturated area) D
r goes to infinity. The degeneracy of (1) and (3) leads to the fact that the solutions must be understood in the sense of distribution as proposed in [4–7].(a)D
r
(
x
) function. (b) ψ
(
x
) function.
(a)
(b)Numerous papers have been published on numerical scheme for the non-multiscale degenerate parabolic problem (including Richards’ equation). For example, in [5–9], the authors firstly regularized the degenerated problem then used the FEM or mixed FEM for spatial discretization. In [10–13], considering the time discretization aspect, the authors developed some schemes for time discretization, such as linear time discretization, adaptive time stepping, or relaxation scheme. The relaxation-iteration scheme analyzed in [13] for the fast diffusion case is a fixed point iteration which was proposed in [8] (also see [12] for the mixed finite element context).Here, we should point out the work of [6, 7]. In [6], the authors considered a class of degenerate parabolic equations including Stefan problem, porous-medium problem, and nonstationary filtration problem. Combining with a regularization approach (when necessary), the authors developed their fully discrete scheme including FEM in spatial and Backward-Euler semi-implicit scheme in time. At the same time, the effects of numerical integration and domain change were taken into account. In [7], under the relations (2), Backward-Euler implicit scheme in time with a regularization step for non-multiscale Richards’ equation (3) was established and analyzed.On the other hand, for this kind of multiscale problem, it is impossible to account explicitly for the spatial variability at fine scale because of the computational resource limitations in realistic situation. So, traditional numerical methods such as standard finite element methods (FEM) and finite difference methods (FDM) are generally not capable of solving this multiscale problem directly. In recent years, a number of multiscale numerical methods, such as multiscale finite element method (MsFEM) [14], heterogeneous multiscale method (HMM) [15], and upscaling method [16], have been proposed to solve the general multiscale problems based on similar ideas. Among them, the so-called heterogeneous multiscale method (HMM) has proved to be an efficient tool to assemble information from microscale problems in order to perform macroscale simulations. In [17], the author discussed the modeling and analysis of finite element methods (such as hybrid methods, coupling spectral or discontinuous Galerkin methods with FEM) for multiscale problems constructed in the framework of the HMM for multiscale PDEs (such as elliptic and parabolic equations).To our knowledge, there are few papers about the numerical analysis of multiscale Richards’ equation of van Genuchten-Mualem model. Based on the frame of HMM-FEM, we will use the idea of [6, 7] to develop a fully discrete method for multiscale Richards’ equation (5). In order to overcome the above difficulties, we regularize the equation firstly; secondly, a fully discrete multiscale numerical method based on heterogeneous multiscale method (HMM) [15] and FEM is developed on a macroscale mesh; at the same time, we use the numerical quadrature for computing the integral over each element. In the paper, we are only interested in (3).Noticing thatD
(
x
/
ϵ
,
Θ
) is separable, we set ψ
(
Θ
)
=
∫
0
Θ
D
r
(
s
)
d
s (see Figure 1(b)). By Kirchhoff transform [4], we rewrite (3) as
(5)
∂
t
Θ
-
∇
(
D
s
(
x
ϵ
)
∇
ψ
(
Θ
)
+
K
z
→
(
x
ϵ
,
Θ
)
)
=
0
with ψ
′
(
Θ
)
=
D
r
(
Θ
). So, when ψ
′
(
Θ
)
→
0
(
Θ
→
0
), (5) may degenerate to hyperbolic equation; on the other hand, ψ
′
(
Θ
)
→
+
∞
(
Θ
→
1
), and (5) may transfer to elliptic equation.Notations. Ω
⊂
R
d
(
d
≥
1
) is a bounded polyhedral domain with boundary ∂
Ω. Set Q
=
Ω
×
[
0
,
T
], where T
>
0 is fixed. For defining a solution in weak sense, we let (
·
,
·
) stand for the inner product on L
2
(
Ω
) or the duality pairing between H
-
1
(
Ω
) and H
0
1
(
Ω
), ∥
·
∥
0 stand for the norm in L
2
(
Ω
), and ∥
·
∥
-
1 and ∥
·
∥
1 stand for the norm in H
-
1
(
Ω
) and H
0
1
(
Ω
).


## 2. Formulation of the Problem

In this paper, we will consider the following general nonlinear parabolic equation which includes Richards’ equation (5):
(6)
∂
t
u
ϵ
-
∇
(
a
(
x
ϵ
)
∇
β
(
u
ϵ
)
+
g
(
x
ϵ
,
u
ϵ
)
)
=
0
,
in
Q
,
u
ϵ
(
t
=
0
)
=
u
0
(
x
)
,
in
Ω
,
u
ϵ
=
0
,
on
⁡
∂
Ω
×
(
0
,
T
]
.Assumption 1.
(H1)β is a maximal monotone graph in [
0,1
]
×
R and β
(
0
)
=
0.
(H2) For alls
∈
(
0,1
], β
′
(
s
)
>
0 and β
′
(
s
) has the following asymptotic behavior:
(7)
∃
α
>
0
,
0
<
γ
<
1
,
s
.
t
.
β
′
(
s
)
~
O
(
s
α
)
(
s
⟶
0
)
,
s
.
t
.
β
′
(
s
)
~
O
(
(
1
-
s
)
-
γ
)
(
s
⟶
1
)
.
(H3)0
≤
u
0
(
x
)
≤
1 almost everywhere.
(H4) The matrixa
(
x
/
ϵ
)
=
(
a
i
j
(
x
/
ϵ
)
)
i
,
j
=
1
d is bounded and uniformly elliptic. There exist constants C
a
_ and C
a
¯ such that C
a
_
∥
ξ
∥
2
≤
(
a
(
y
)
ξ
,
ξ
)
≤
C
a
¯
∥
ξ
∥
2, ∀
ξ
∈
R
d.
(H5)g
(
y
,
s
)
:
R
n
×
R
↦
R
n is bounded, uniformly Lipschitz continuous and satisfying
(8)
|
g
(
·
,
u
1
)
-
g
(
·
,
u
2
)
|
2
≤
C
(
u
1
-
u
2
)
(
β
(
u
1
)
-
β
(
u
2
)
)
.
A weak solution of problem (6) is defined as follows.Definition 2.
Findu
ϵ
∈
H
1
(
0
,
T
;
H
-
1
(
Ω
)
) such that u
ϵ
(
0
)
=
u
0 and v
ϵ
=
β
(
u
ϵ
)
∈
L
2
(
0
,
T
;
H
0
1
(
Ω
)
) and, ∀
φ
∈
L
2
(
0
,
T
;
H
0
1
(
Ω
)
),
(9)
(
∂
t
u
ϵ
,
φ
)
Q
+
(
a
ϵ
(
x
ϵ
,
u
ϵ
)
∇
β
(
u
ϵ
)
,
∇
φ
)
Q
+
(
g
ϵ
(
x
ϵ
,
u
ϵ
)
,
∇
φ
)
Q
=
0
.Remark 3.
The existence, uniqueness, and regularity results for the above problem can be found in [4]. The solution of problem (6) may denote the reduced saturation and therefore should be bounded by 0 and 1. In the setting stated above, the maximum principle holds for problem (6), so the solution remains bounded by 0 and 1.Remark 4.
According to the constitution relations of the van Genuchten-Mualem, the parameter in (H2) is determined asα
=
(
1
/
2
)
+
(
1
/
m
) and γ
=
m (m
∈
(
0,1
)). In Richards’ equation, the convective function g
(
y
,
u
) denotes the permeability K
s
(
y
)
K
r
(
u
)
e
→
z. And it is easy to verify that (
K
r
)
′
(
u
)
~
O
(
u
(
2
/
m
)
-
(
1
/
2
)
) and β
′
(
u
)
~
O
(
u
(
1
/
m
)
+
(
1
/
2
)
) when u approaches to zero. Also, it is easy to verify that condition (8) is valid when the second variable approaches to zero. The other cases are obvious.Let0
<
δ
<
1 be a small perturbed parameter. We approximate β by, for all s
∈
[
0,1
],
(10)
β
δ
(
s
)
=
∫
0
s
‍
max
⁡
{
δ
,
min
⁡
{
1
δ
,
β
′
(
s
)
}
}
d
s
.
So, we get
(11)
β
δ
′
(
s
)
=
{
δ
,
if
0
≤
s
≤
C
δ
1
/
α
,
β
′
(
s
)
,
otherwise
,
1
δ
,
if
1
-
C
δ
1
/
γ
≤
s
≤
1
.
Hence, β
δ
-
1 is also a maximal monotone graph in R
×
[
0,1
] and it follows that
(12)
δ
≤
(
β
δ
-
1
)
′
(
x
)
≤
1
δ
,
|
β
δ
-
1
(
x
)
-
β
-
1
(
x
)
|
≤
C
δ
μ
,
μ
=
min
⁡
{
1
γ
,
1
α
}
.LetT
H
=
{
K
}
H be a regular triangle decomposition of Ω, where H stands for the mesh size. Define the functional space
(13)
X
H
=
{
χ
is
linear
for
all
K
∈
T
H
}
,
X
H
0
=
{
χ
is
linear
for
all
K
∈
T
H
;
χ
=
0
on
⁡
∂
Ω
}
.In order to simplify the computation, we will use the numerical integration scheme as in [6]. Let Π
H be the local linear interpolant operator; then,
(14)
(
ϕ
,
φ
)
H
,
K
≐
∫
K
Π
H
(
ϕ
φ
)
d
x
,
for any piecewise uniformly continuous functions ϕ and φ. So, for all W, V
∈
X
H and every element K, the numerical integration scheme satisfies that [6]
(15)
E
(
W
,
V
)
≐
|
(
W
,
V
)
K
-
(
W
,
V
)
H
,
K
|
≤
C
H
∥
W
∥
L
2
(
K
)
∥
V
∥
H
1
(
K
)
,
(16)
E
(
W
,
V
)
≐
|
(
W
,
V
)
K
-
(
W
,
V
)
H
,
K
|
≤
C
H
2
∥
W
∥
H
1
(
K
)
∥
V
∥
H
1
(
K
)
,
(17)
C
1
∥
W
∥
L
2
(
K
)
2
≤
(
W
,
W
)
H
,
K
≤
C
2
∥
W
∥
L
2
(
K
)
2
,
(18)
|
(
W
,
V
)
H
,
K
|
≤
C
2
∥
W
∥
L
2
(
K
)
∥
V
∥
L
2
(
K
)
and, for any b
(
x
)
∈
W
1
,
∞
(
K
) [18],
(19)
E
(
b
W
,
V
)
≤
C
H
∥
b
∥
W
1
,
∞
(
K
)
∥
W
∥
H
1
(
K
)
∥
V
∥
H
1
(
K
)
.
I
H stands for the piecewise linear interpolant operator and satisfies
(20)
∥
I
H
φ
-
φ
∥
L
2
(
Ω
)
≤
C
H
∥
∇
φ
∥
L
2
(
Ω
)
,
φ
∈
H
1
(
Ω
)
,
and E
H stands for the discrete L
2 projection operator; more precisely, for any z
∈
L
2
(
Ω
), E
H
z
∈
X
H and satisfies
(21)
∥
E
H
z
-
z
∥
H
-
1
(
Ω
)
≤
C
H
∥
z
∥
L
2
(
Ω
)
.
Let x
r and ω
K
,
r (r
=
1
,
L
¯) denote the quadrature nodes and weights in K. Let N be an integer and let τ
=
T
/
N. Our fully discrete formulation for problem (6) is read as follows: find v
i
,
H
∈
X
H
0,  u
i
,
H
=
I
H
β
δ
-
1
(
v
i
,
H
)
∈
X
H (i
=
1
,
N
¯), such that
(22)
∑
K
‍
(
u
i
,
H
-
u
i
-
1
,
H
τ
,
χ
)
H
,
K
+
∑
K
‍
A
H
(
v
i
,
H
,
χ
)
+
∑
K
‍
G
H
(
u
i
,
H
,
χ
)
=
0
,
∀
χ
∈
X
H
0
,
v
0
,
H
=
I
H
β
(
u
0
)
,
u
0
,
H
=
E
H
u
0
,
where
(23)
A
H
(
v
i
,
H
,
χ
)
≐
∑
r
=
1
L
‍
ω
K
,
r
|
I
l
(
x
r
)
|
∫
I
l
(
x
r
)
‍
a
(
x
ϵ
)
∇
R
1
(
v
i
,
H
)
∇
R
1
(
χ
)
d
x
≐
(
a
^
∇
v
i
,
H
,
∇
χ
)
H
,
K
,
(24)
G
H
(
u
i
,
H
,
χ
)
≐
∑
r
=
1
L
‍
ω
K
,
r
|
I
l
(
x
r
)
|
∫
I
l
(
x
r
)
‍
g
(
x
ϵ
,
u
i
,
H
)
∇
R
1
(
χ
)
d
x
≐
(
g
^
(
u
i
,
H
)
,
∇
χ
)
H
,
K
,
and the operator R
1 is defined by the following problem:
(25)
-
∇
(
a
ϵ
(
x
)
∇
R
1
(
χ
)
)
=
0
in
I
l
(
x
r
)
,
hhhhhh
R
1
(
χ
)
=
χ
r
on
⁡
∂
I
l
(
x
r
)
,
where χ
r
=
χ
(
x
r
)
+
χ
′
(
x
r
)
(
x
-
x
r
) and the cell I
l
(
x
r
) is a square of size l centered at x
r. So far, our fully discrete HMM-FEM is settled down and the porous media do not have to be periodic.


## 3. Some Results on the Homogenized Problem

In this section, we will review the results of homogenization and the analysis of the discrete homogenized problem. Under the assumptions (H1)–(H5) and if the functiona
(
y
) and g
(
y
,
s
) are periodic in y with unit square Y. In [19], the authors have established the homogenization theory for (6). It has been shown that, in suitable topology space, the solutions u
ϵ of problem (6) converge to the solution of the following problem as ϵ
→
0.
(26)
∂
t
u
-
∇
(
a
*
∇
β
(
u
)
+
g
*
(
u
)
)
=
0
,
in
Q
,
u
(
t
=
0
)
=
u
0
(
x
)
,
in
Ω
,
u
=
0
,
on
⁡
∂
Ω
×
(
0
,
T
]
,
where a
*
=
(
a
i
j
*
)
1
≤
i
,
j
≤
d and g
*
=
(
g
i
*
)
i
=
1
d are defined as
(27)
(
a
*
)
i
j
=
∫
Y
‍
(
a
i
j
(
y
)
+
a
i
k
(
y
)
∂
Λ
k
(
y
)
∂
y
j
)
d
y
,
(28)
g
*
(
s
)
=
∫
Y
‍
(
g
(
y
,
s
)
+
a
(
y
)
∇
y
η
(
y
,
s
)
)
d
y
,
and Λ
k and η are the periodic solutions of following cell problems, respectively.
(29)
-
∇
y
·
(
a
(
y
)
∇
y
Λ
k
(
y
)
)
=
∇
y
·
(
a
(
y
)
·
e
k
)
,
in
Y
,
∫
Y
Λ
k
d
y
=
0
,
-
∇
y
·
(
a
(
y
)
∇
y
η
(
y
,
s
)
)
=
∇
y
·
(
g
(
y
,
s
)
)
,
in
Y
,
∫
Y
‍
η
d
y
=
0
,
and a
* and g
* inherit the property of a and g [16].Remark 5.
According to Alt and Luckhaus [4], we have at least u
∈
L
∞
(
0
,
T
;
L
1
(
Ω
)
) and u
t
∈
L
2
(
0
,
T
;
H
-
1
(
Ω
)
). Moreover, the maximum principle leads to u
∈
L
∞
(
0
,
T
;
L
∞
(
Ω
)
). We therefore conclude that u
∈
C
(
0
,
T
;
L
2
(
Ω
)
). This gives us u
(
t
,
·
) pointwise for every t
∈
[
0
,
T
].The problem (26) can be read as follows.Definition 6.
For every time interval[
t
i
-
1
,
t
i
]
⊂
[
0
,
T
], find u
(
t
i
)
∈
L
2
(
Ω
) such that ∫
t
i
-
1
t
i
‍
β
(
u
(
t
)
)
d
t
∈
H
0
1
(
Ω
) and, for all φ
∈
H
0
1
(
Ω
),
(30)
(
u
(
t
i
)
-
u
(
t
i
-
1
)
,
φ
)
Ω
+
(
∫
t
i
-
1
t
i
a
*
∇
β
(
u
)
d
t
,
∇
φ
)
Ω
+
(
∫
t
i
-
1
t
i
g
*
(
u
)
d
t
,
∇
φ
)
Ω
=
0
.We also give out the fully discrete scheme for the problem (26) by a regularization procedure.Definition 7.
FindV
i
∈
X
H
0 such that U
i
=
I
H
β
δ
-
1
(
V
i
)
∈
X
H and, for all W
∈
X
H
0,
(31)
(
U
i
-
U
i
-
1
,
W
)
H
+
τ
(
a
*
∇
V
i
,
∇
W
)
H
+
τ
(
I
H
g
*
(
U
i
)
,
∇
W
)
H
=
0
,
V
0
=
I
H
β
(
u
0
)
,
U
0
=
E
H
u
0
,
where (
U
,
W
)
H
=
∑
K
∈
T
H
‍
(
U
,
W
)
H
,
K.In the rest of this section, we will give out a prior estimate for problem (31) and the error estimate of V
i
-
β
(
u
(
t
i
)
) and U
i
-
u
(
t
i
).Theorem 8.
Under the assumptions (H1) to (H5), there exist a positive constantC independent of H
,
δ
,
ϵ such that
(32)
max
⁡
1
≤
i
≤
N
∥
U
i
∥
L
∞
(
Ω
)
+
δ
max
⁡
1
≤
i
≤
N
∥
V
i
∥
L
2
(
Ω
)
2
+
∑
i
=
1
N
‍
τ
∥
V
i
∥
H
1
(
Ω
)
2
≤
C
,
∑
i
=
1
N
‍
δ
(
∥
V
i
-
V
i
-
1
∥
0
2
+
∥
U
i
-
U
i
-
1
∥
0
2
)
≤
C
.Proof.
LetW
=
V
i in (31) and sum it over i from 1 to N.
(33)
∑
i
=
1
N
‍
(
U
i
-
U
i
-
1
,
V
i
)
H
+
∑
i
=
1
N
‍
τ
(
a
*
∇
V
i
,
∇
V
i
)
H
+
∑
i
=
1
N
‍
τ
(
I
H
g
*
(
U
i
)
,
∇
V
i
)
H
=
0
.
We will estimate each term separately. Following the idea in [6], we have
(34)
∑
i
=
1
N
‍
(
U
i
-
U
i
-
1
,
V
i
)
H
=
∑
i
=
1
N
‍
(
I
H
β
δ
-
1
(
V
i
)
-
I
H
β
δ
-
1
(
V
i
-
1
)
,
V
i
)
H
≥
δ
∥
V
N
∥
0
2
-
C
.
Using (H4), (H5) and Poincaré inequality, we have(35)
∑
i
=
1
N
‍
τ
(
a
*
∇
V
i
,
∇
V
i
)
H
≥
C
∑
i
=
1
N
‍
τ
∥
∇
V
i
∥
0
2
≥
C
∑
i
=
1
N
‍
τ
∥
V
i
∥
1
2
,
|
∑
i
=
1
N
‍
τ
(
I
H
g
*
(
U
i
)
,
∇
V
i
)
H
|
≤
∑
i
=
1
N
‍
τ
∥
I
H
g
*
(
U
i
)
∥
0
∥
∇
V
i
∥
0
≤
C
+
η
∑
i
=
1
N
‍
τ
∥
∇
V
i
∥
0
2
.
Since U
i
=
I
H
β
δ
-
1
(
V
i
), we have max
⁡
1
≤
i
≤
N
∥
U
i
∥
L
∞
(
Ω
)
≤
C. Combining all the terms and choosing η properly, we complete the first part.
Again, we letW
=
(
V
i
-
V
i
-
1
) in (31) and sum it over i from 1 to N.
(36)
∑
i
=
1
N
‍
(
U
i
-
U
i
-
1
,
V
i
-
V
i
-
1
)
H
+
∑
i
=
1
N
‍
τ
(
a
*
∇
V
i
,
∇
(
V
i
-
V
i
-
1
)
)
H
+
∑
i
=
1
N
‍
τ
(
I
H
g
*
(
U
i
)
,
∇
(
V
i
-
V
i
-
1
)
)
H
=
0
.
Considering that V
i
=
I
H
β
δ
-
1
(
U
i
), we have
(37)
∑
i
=
1
N
‍
(
U
i
-
U
i
-
1
,
V
i
-
V
i
-
1
)
H
≥
C
∑
i
=
1
N
‍
δ
(
∥
U
i
-
U
i
-
1
∥
0
2
+
∥
V
i
-
V
i
-
1
∥
0
2
)
.
By (H4), (H5), and (32), it is easy to get
(38)
|
∑
i
=
1
N
‍
τ
(
a
*
∇
V
i
,
∇
(
V
i
-
V
i
-
1
)
)
H
h
+
∑
i
=
1
N
‍
τ
(
I
H
g
*
(
U
i
)
,
∇
(
V
i
-
V
i
-
1
)
)
H
|
≤
C
.
So, we finish the second part of the theorem.Define a regular G-operatorH
-
1
→
H
0
1 and G
H-operator H
-
1
→
X
H
0 as in [6], for all φ
∈
H
-
1,
(39)
(
a
*
∇
G
φ
,
∇
ψ
)
=
(
φ
,
ψ
)
,
∀
ψ
∈
H
0
1
,
(
a
*
∇
G
φ
,
∇
χ
)
=
(
a
*
∇
G
H
φ
,
∇
χ
)
,
∀
χ
∈
X
H
0
,
and there is
(40)
∥
(
G
-
G
H
)
φ
∥
H
s
≤
C
H
2
-
r
-
s
∥
φ
∥
H
-
r
,
0
≤
r
,
s
≤
1
.
Firstly, we have
(41)
∥
φ
∥
-
1
=
sup
⁡
ψ
∈
H
0
1
|
(
φ
,
ψ
)
|
∥
∇
ψ
∥
,
(42)
C
a
_
∥
∇
G
φ
∥
2
≤
(
a
*
∇
G
φ
,
∇
G
φ
)
=
(
φ
,
G
φ
)
≤
∥
φ
∥
-
1
∥
∇
G
φ
∥
.
On the other hand,
(43)
∥
φ
∥
-
1
=
sup
⁡
ψ
∈
H
0
1
|
(
φ
,
ψ
)
|
∥
∇
ψ
∥
=
sup
⁡
ψ
∈
H
0
1
|
(
a
*
∇
G
φ
,
∇
ψ
)
|
∥
∇
ψ
∥
≤
C
a
¯
∥
∇
G
φ
∥
,
so we have ∥
φ
∥
-
1 equivalent to ∥
∇
G
φ
∥. We also introduce the following two identities:
(44)
∑
i
=
1
n
‍
(
a
i
-
a
i
-
1
,
b
i
)
=
a
n
b
n
-
a
0
b
0
-
∑
i
=
1
n
‍
a
i
-
1
(
b
i
-
b
i
-
1
)
,
(45)
∑
i
=
1
n
‍
(
a
i
-
a
i
-
1
,
a
i
)
=
(
a
n
)
2
-
(
a
0
)
2
+
∑
i
=
1
n
‍
(
a
i
-
a
i
-
1
)
2
.Theorem 9.
u
i and U
i are the solutions of problems (26) and (31) at time t
i, respectively. Denote that u
¯
i
=
(
1
τ
)
∫
I
i
u
d
t, e
u
i
=
U
i
-
u
¯
i, and e
v
i
=
V
i
-
(
1
/
τ
)
∫
I
i
β
(
u
)
d
t. Under the assumptions (H1) to (H5), there exists a positive constant C independent of H, δ, ϵ, τ such that
(46)
max
⁡
1
≤
i
≤
N
∥
e
u
i
∥
-
1
2
+
∑
i
=
1
N
‍
∥
e
u
i
-
e
u
i
-
1
∥
-
1
2
+
∑
i
=
1
N
‍
∫
I
i
‍
(
V
i
-
β
(
u
)
,
β
-
1
(
V
i
)
-
u
(
t
)
)
Ω
d
t
≤
C
(
H
+
H
2
+
δ
2
μ
+
δ
μ
+
τ
+
H
δ
+
H
4
τ
δ
2
+
H
2
δ
2
+
H
2
τ
δ
)
.Remark 10.
The proof method of Theorem9 can be found in [6, 12]. However, there is a little difference between our result and the result in [6, 12] because we consider the two degenerate points case.


## 4. Main Result

In this section, the main task is to give out the error betweenV
i and v
i
,
H. Here, V
i and v
i
,
H are the solutions of problems (31) and (22), respectively. Before the proof of our main result, the following useful lemma will be introduced [15, 16].Lemma 11.
a
*, g
* are defined in (27) and (28) and a
^, g
^ are defined in (23) and (24), respectively. Then, for all s
∈
R,
(47)
|
a
*
(
s
)
-
a
^
(
s
)
|
≤
C
ϵ
l
,
|
g
*
(
s
)
-
g
^
(
s
)
|
≤
C
ϵ
l
.Then, we have the main result of the paper.Theorem 12.
u
i
,
H, v
i
,
H and U
i, V
i are the solutions of problems (22) and (31) at time t
i, respectively. Denote e
v
,
i
H
=
v
i
,
H
-
V
i
∈
X
H
0 and e
u
,
i
H
=
u
i
,
H
-
U
i
∈
X
H. Then, under assumptions (H1)–(H5), there exists a positive constant C independent of H, δ, ϵ, τ such that
(48)
max
⁡
1
≤
i
≤
N
∥
e
u
,
i
H
∥
-
1
2
+
∑
i
=
1
N
‍
∥
e
u
,
i
H
-
e
u
,
i
-
1
H
∥
-
1
2
+
τ
∑
i
=
1
N
‍
(
v
i
,
H
-
V
i
,
β
-
1
(
v
i
,
H
)
-
β
-
1
(
V
i
)
)
Ω
≤
C
(
H
+
H
δ
+
H
2
δ
2
+
δ
μ
+
δ
2
μ
+
ϵ
2
l
2
)
.Proof.
Subtract (31) from (22) and notice the notation ∑
K
‍
(
·
,
·
)
H
,
K
=
(
·
,
·
)
H; we have
(49)
(
u
i
,
H
-
u
i
-
1
,
H
-
U
i
+
U
i
-
1
,
χ
)
H
+
τ
(
A
H
(
v
ϵ
,
H
,
χ
)
-
(
a
*
∇
V
i
,
∇
χ
)
H
)
+
τ
(
G
H
(
u
ϵ
,
H
,
χ
)
-
(
g
*
(
U
i
)
,
∇
χ
)
H
)
=
0
;
that is,
(50)
(
u
i
,
H
-
u
i
-
1
,
H
-
U
i
+
U
i
-
1
,
χ
)
H
+
τ
(
(
a
^
∇
v
i
,
H
,
∇
χ
)
H
-
(
a
*
∇
V
i
,
∇
χ
)
H
)
+
τ
(
(
g
^
(
u
i
,
H
)
,
∇
χ
)
H
-
(
g
*
(
U
i
)
,
∇
χ
)
H
)
=
0
.
Let χ
=
G
H
e
u
,
i
H
∈
X
H
0 and sum i
=
1
⋯
N; then we have
(51)
∑
i
=
1
N
‍
(
u
i
,
H
-
u
i
-
1
,
H
-
U
i
+
U
i
-
1
,
G
H
e
u
,
i
H
)
H
+
τ
∑
i
=
1
N
‍
(
(
a
^
∇
v
i
,
H
-
a
*
∇
V
i
,
∇
G
H
e
u
,
i
H
)
H
)
+
τ
∑
i
=
1
N
‍
(
(
g
^
(
u
i
,
H
)
-
g
*
(
U
i
)
,
∇
G
H
e
u
,
i
H
)
H
)
=
0
.
Denote the above equality by T
1
+
T
2
+
T
3
=
0.
For the termT
1, noticing that a
* is also positive and bounded, we use (15), (45), and the definition of G
H to get
(52)
T
1
=
∑
i
=
1
N
‍
(
e
u
,
i
H
-
e
u
,
i
-
1
H
,
G
H
e
u
,
i
H
)
H
=
∑
i
=
1
N
‍
(
e
u
,
i
H
-
e
u
,
i
-
1
H
,
G
H
e
u
,
i
H
)
+
∑
i
=
1
N
‍
∑
K
‍
E
(
e
u
,
i
H
-
e
u
,
i
-
1
H
,
G
H
e
u
,
i
H
)
≥
∑
i
=
1
N
‍
(
a
*
∇
G
H
(
e
u
,
i
H
-
e
u
,
i
-
1
H
)
,
∇
G
H
e
u
,
i
H
)
-
C
H
∑
i
=
1
N
‍
∥
e
u
,
i
H
-
e
u
,
i
-
1
H
∥
0
∥
G
H
e
u
,
i
H
∥
1
≥
(
a
*
∇
G
H
e
u
,
N
H
,
∇
G
H
e
u
,
N
H
)
+
∑
i
=
1
N
‍
(
a
*
∇
G
H
(
e
u
,
i
H
-
e
u
,
i
-
1
H
)
,
∇
G
H
(
e
u
,
i
H
-
e
u
,
i
-
1
H
)
)
-
(
a
*
∇
G
H
e
u
,
0
H
,
∇
G
H
e
u
,
0
H
)
-
η
∑
i
=
1
N
‍
∥
e
u
,
i
H
-
e
u
,
i
-
1
H
∥
-
1
2
-
C
(
η
)
H
2
∑
i
=
1
N
‍
∥
e
u
,
i
H
∥
-
1
2
≥
C
a
_
∥
e
u
,
N
H
∥
-
1
2
-
C
a
¯
∥
e
u
,
0
H
∥
-
1
2
+
(
C
a
_
-
η
)
∑
i
=
1
N
‍
∥
e
u
,
i
H
-
e
u
,
i
-
1
H
∥
-
1
2
-
C
(
η
)
H
2
∑
i
=
1
N
‍
∥
e
u
,
i
H
∥
-
1
2
.
Rewrite the term T
2 as
(53)
T
2
=
τ
∑
i
=
1
N
‍
(
a
*
∇
e
v
,
i
H
,
∇
G
H
e
u
,
i
H
)
+
τ
∑
i
=
1
N
‍
∑
K
‍
E
(
a
*
∇
e
v
,
i
H
,
∇
G
H
e
u
,
i
H
)
+
τ
∑
i
=
1
N
‍
(
(
a
^
-
a
*
)
∇
v
i
,
H
,
∇
G
H
e
u
,
i
H
)
H
=
τ
∑
i
=
1
N
‍
(
e
v
,
i
H
,
e
u
,
i
H
)
+
τ
∑
i
=
1
N
‍
∑
K
‍
E
(
a
*
∇
e
v
,
i
H
,
∇
G
H
e
u
,
i
H
)
+
τ
∑
i
=
1
N
‍
(
(
a
^
-
a
*
)
∇
v
i
,
H
,
∇
G
H
e
u
,
i
H
)
H
≐
T
21
+
T
22
+
T
23
.
Applying (20), (12), and (32), we get
(4)
T
21
=
τ
∑
i
=
1
N
‍
(
v
i
,
H
-
V
i
,
β
-
1
(
v
i
,
H
)
-
β
-
1
(
V
i
)
)
+
τ
∑
i
=
1
N
‍
(
e
v
,
i
H
,
I
H
β
δ
-
1
(
v
i
,
H
)
-
β
-
1
(
v
i
,
H
)
+
τ
∑
i
=
1
N
‍
h
h
+
I
H
β
δ
-
1
(
V
i
)
-
β
-
1
(
V
i
)
)
≥
τ
∑
i
=
1
N
‍
(
v
i
,
H
-
V
i
,
β
-
1
(
v
i
,
H
)
-
β
-
1
(
V
i
)
)
-
C
τ
∑
i
=
1
N
‍
∥
e
v
,
i
H
∥
0
(
H
∥
∇
β
δ
-
1
(
v
i
,
H
)
∥
-
C
τ
∑
i
=
1
N
‍
∥
e
v
,
i
H
∥
0
h
h
+
H
∥
∇
β
δ
-
1
(
V
i
)
∥
+
δ
μ
)
≥
τ
∑
i
=
1
N
‍
(
v
i
,
H
-
V
i
,
β
-
1
(
v
i
,
H
)
-
β
-
1
(
V
i
)
)
-
C
τ
∑
i
=
1
N
‍
∥
e
v
,
i
H
∥
0
(
H
δ
∥
∇
v
i
,
H
∥
+
H
δ
∥
∇
V
i
∥
+
δ
μ
)
≥
τ
∑
i
=
1
N
‍
(
v
i
,
H
-
V
i
,
β
-
1
(
v
i
,
H
)
-
β
-
1
(
V
i
)
)
-
C
(
H
δ
+
δ
μ
)
,
|
T
22
|
≤
C
H
τ
∑
i
=
1
N
‍
∥
a
*
∥
1
,
∞
∥
∇
e
v
,
i
H
∥
0
∥
∇
G
H
e
u
,
i
H
∥
0
≤
C
∑
i
=
1
N
‍
τ
(
H
2
∥
e
v
,
i
H
∥
1
2
+
∥
e
u
,
i
H
∥
1
2
)
≤
C
(
H
2
+
∑
i
=
1
N
‍
τ
∥
e
u
,
i
H
∥
1
2
)
,
|
T
23
|
≤
C
ϵ
l
τ
∑
i
=
1
N
‍
∥
∇
v
i
,
H
∥
0
∥
e
u
,
i
H
∥
-
1
≤
C
(
ϵ
l
)
2
+
C
τ
∑
i
=
1
N
‍
∥
e
u
,
i
H
∥
-
1
2
.
For the terms T
22, T
23, we have used (19), (18), and Lemma 11.
Similar toT
2, the term T
3 can be estimated as the following:
(4)
T
3
=
τ
∑
i
=
1
N
‍
(
g
^
(
u
i
,
H
)
-
g
*
(
U
i
)
,
∇
G
H
e
u
,
i
H
)
H
=
τ
∑
i
=
1
N
‍
(
g
^
(
u
i
,
H
)
-
g
*
(
u
i
,
H
)
,
∇
G
H
e
u
,
i
H
)
H
+
τ
∑
i
=
1
N
‍
(
g
*
(
u
i
,
H
)
-
g
*
(
U
i
)
,
∇
G
H
e
u
,
i
H
)
H
≐
T
31
+
T
32
,
|
T
31
|
≤
C
ϵ
l
τ
∑
i
=
1
N
‍
∥
e
u
,
i
H
∥
-
1
≤
C
(
ϵ
l
)
2
+
τ
∑
i
=
1
N
‍
∥
e
u
,
i
H
∥
-
1
2
,
T
32
=
τ
∑
i
=
1
N
‍
(
g
*
(
u
i
,
H
)
-
g
*
(
U
i
)
,
∇
G
H
e
u
,
i
H
)
+
τ
∑
i
=
1
N
‍
E
(
g
*
(
u
i
,
H
)
-
g
*
(
U
i
)
,
∇
G
H
e
u
,
i
H
)
,
|
T
32
|
≤
τ
∑
i
=
1
N
‍
(
∥
g
*
(
u
i
,
H
)
-
g
*
(
U
i
)
∥
0
≤
τ
∑
i
=
1
N
‍
hh
+
C
H
(
∥
∇
g
*
(
u
i
,
H
)
∥
0
+
∥
∇
g
*
(
U
i
)
∥
0
)
)
×
∥
e
u
,
i
H
∥
-
1
.
Noticing u
i
,
H
=
I
H
β
δ
-
1
(
v
i
,
H
) (also U
i) and (H5), it follows that
(56)
∥
g
*
(
u
i
,
H
)
-
g
*
(
U
i
)
∥
0
2
≤
∥
g
*
(
β
-
1
(
v
i
,
H
)
)
-
g
*
(
β
-
1
(
V
i
)
)
∥
0
2
+
∥
g
*
(
u
i
,
H
)
-
g
*
(
β
-
1
(
v
i
,
H
)
)
∥
0
2
+
∥
g
*
(
β
-
1
(
V
i
)
)
-
g
*
(
U
i
)
∥
0
2
≤
(
v
i
,
H
-
V
i
,
β
-
1
(
v
i
,
H
)
-
β
-
1
(
V
i
)
)
Ω
+
C
H
2
δ
2
(
∥
∇
v
i
,
H
∥
2
+
∥
∇
V
i
∥
2
)
+
C
δ
2
μ
≤
(
v
i
,
H
-
V
i
,
β
-
1
(
v
i
,
H
)
-
β
-
1
(
V
i
)
)
Ω
+
C
(
H
2
δ
2
+
δ
2
μ
)
,
∥
∇
g
*
(
u
i
,
H
)
∥
0
=
∥
∇
g
*
(
I
H
β
δ
-
1
(
v
i
,
H
)
)
∥
0
≤
C
1
δ
∥
∇
v
i
,
H
∥
0
.
So, we have
(57)
|
T
32
|
≤
η
τ
∑
i
=
1
N
‍
(
v
i
,
H
-
V
i
,
β
-
1
(
v
i
,
H
)
-
β
-
1
(
V
i
)
)
Ω
+
C
(
H
2
δ
2
+
δ
2
μ
)
+
τ
∑
i
=
1
N
‍
∥
e
u
,
i
H
∥
-
1
2
.
Combining all the terms and choosing the parameter η properly, we have
(58)
∥
e
u
,
N
H
∥
-
1
2
+
∑
i
=
1
N
‍
∥
e
u
,
i
H
-
e
u
,
i
-
1
H
∥
-
1
2
+
τ
∑
i
=
1
N
‍
(
v
i
,
H
-
V
i
,
β
-
1
(
v
i
,
H
)
-
β
-
1
(
V
i
)
)
Ω
≤
C
(
H
2
+
H
δ
+
H
2
δ
2
+
δ
μ
+
δ
2
μ
+
ϵ
2
l
2
)
+
τ
∑
i
=
1
N
‍
∥
e
u
,
i
H
∥
-
1
2
.
At last, we use the Gronwall inequality to the above formulation to finish the proof completely.


## 5. Numerical Example

In this subsection, we consider the Richards equation under the van Genuchten-Mualem model [2]. Consider
(59)
∂
t
θ
-
∇
·
(
D
(
x
,
θ
)
∇
θ
+
K
(
x
,
θ
)
·
e
→
x
3
)
=
0
∂
t
θ
(
D
(
x
,
θ
)
∇
θ
+
h
)
in
[
0,1
]
2
×
(
0
,
T
)
,
θ
=
θ
s
on
⁡
Γ
t
=
{
x
1
∈
(
0,1
)
,
x
3
=
1
}
,
θ
=
θ
r
on
⁡
Γ
b
=
{
x
1
∈
(
0,1
)
,
x
3
=
0
}
,
(
D
(
x
,
θ
)
∇
θ
+
k
(
x
,
θ
)
e
→
x
3
)
·
n
→
=
0
on
⁡
Γ
l
r
=
{
x
1
=
0,1
;
x
3
∈
(
0,1
)
}
,
θ
(
x
,
0
)
=
θ
r
in
[
0,1
]
2
,
and, here, the constitutive relations are
(60)
θ
(
u
)
=
θ
r
+
(
θ
s
-
θ
r
)
(
1
+
α
s
|
u
|
n
)
-
m
,
K
(
x
,
θ
)
=
K
s
(
x
1
,
x
3
)
η
1
/
2
(
1
-
(
1
-
η
1
/
m
)
m
)
2
,
D
(
x
,
θ
)
=
(
1
-
m
)
K
s
(
x
1
,
x
3
)
α
s
m
(
θ
s
-
θ
r
)
η
(
1
/
2
)
-
m
×
(
(
1
-
η
1
/
m
)
-
m
+
(
1
-
η
1
/
m
)
m
-
2
)
,
where Θ
=
(
θ
-
θ
r
)
/
(
θ
s
-
θ
r
) and K
s, α
s, m, and n are the media parameters. Throughout this subsection, we set m
=
0.5, n
=
2, θ
r
=
0.05, and θ
s
=
0.489. The heterogeneity comes from absolute permeability K
s
(
x
) and α
(
x
).Notice that the Richards equation under the van Genuchten-Mualem model (59) is degenerated at θ
=
θ
r (D
=
0) and θ
=
θ
s (D
=
+
∞) and sharp interface would be developed between saturated and unsaturated regions. Our HMM-FEM scheme can be also used to treat this kind of problems.In our computation, in order to avoid the difficulty of degeneration, the coefficientD
(
x
,
θ
) is replaced by a regularized one as
(61)
D
¯
(
x
,
θ
)
=
{
D
(
x
,
θ
(
Θ
)
)
Θ
∈
[
0.001
,
0.999
]
,
D
(
x
,
θ
(
0.001
)
)
Θ
∈
[
0
,
0.001
)
,
D
(
x
,
θ
(
0.999
)
)
Θ
∈
(
0.999
,
1
]
,
where θ
(
Θ
)
=
(
θ
s
-
θ
r
)
Θ
+
θ
r.


### 5.1. van Genuchten-Mualem Model with Periodic Coefficients

For periodic case, we set(62)
K
s
(
x
1
,
x
3
)
=
1
/
[
2
+
1.8
sin
⁡
(
2
π
(
2
x
3
-
x
1
)
/
ϵ
)
]
117.4
,
α
s
(
x
1
,
x
3
)
=
1
/
[
2
+
1.5
sin
(
2
π
(
2
x
3
-
x
1
)
/
ϵ
)
]
7.6
.
In our simulation, ϵ is chosen to be 1
/
256.The numerical solution of HMM-FEM is obtained on a macroscale grid32
×
32. For each local sample cell I
l, we choose l
=
4
ϵ and solve the microproblems (25) on a grid with size of ϵ
/
16. We compare the HMM-FEM solution with the fine-scale solution which is obtained on a grid 4096
×
4096 by FVM. The time step is also chosen to be τ
=
1
/
1000 in this subsection. We compare the HMM-FEM solution with the fine-scale solution along the cross-section of x
1
=
0.5 (Figure 2(a)).The HMM-FEM solution (star) compared with fine-scale solution (solid line) along the cross-sectionx
1
=
0.5 for periodic case (a) and random case (isotropic) (b), respectively.
(a)
(b)


### 5.2. van Genuchten-Mualem Model with Random Coefficients

For the random model, we only consider the isotropic heterogeneities case. We generate the random log-normal permeability fieldsK
s
(
x
) and α
(
x
) by the same methods used in [16]. The corresponding normal distribution of log
⁡
(
K
s
) is N
(
-
3,0.8
) and the corresponding normal distribution of log
⁡
(
α
s
) is N
(
-
2,0.8
). The correlation lengths of both k
s and α
s are l
1
=
l
3
=
0.01 in x
1 and x
3 directions. The solution of HMM-FEM is obtained on a macroscale grid 32
×
32. For each local sample cell I
l, we choose l
=
1
/
64 and solve the microproblems (25) on a grid with size of δ
/
16. We compare the HMM-FEM solution with the fine-scale solution which is obtained by solving (59) on a grid 1024
×
1024 by FVM. We compare the HMM-FEM solution with the fine-scale solution along the cross-section x
1
=
0.5 (Figure 2(b)).


## 5.1. van Genuchten-Mualem Model with Periodic Coefficients

For periodic case, we set(62)
K
s
(
x
1
,
x
3
)
=
1
/
[
2
+
1.8
sin
⁡
(
2
π
(
2
x
3
-
x
1
)
/
ϵ
)
]
117.4
,
α
s
(
x
1
,
x
3
)
=
1
/
[
2
+
1.5
sin
(
2
π
(
2
x
3
-
x
1
)
/
ϵ
)
]
7.6
.
In our simulation, ϵ is chosen to be 1
/
256.The numerical solution of HMM-FEM is obtained on a macroscale grid32
×
32. For each local sample cell I
l, we choose l
=
4
ϵ and solve the microproblems (25) on a grid with size of ϵ
/
16. We compare the HMM-FEM solution with the fine-scale solution which is obtained on a grid 4096
×
4096 by FVM. The time step is also chosen to be τ
=
1
/
1000 in this subsection. We compare the HMM-FEM solution with the fine-scale solution along the cross-section of x
1
=
0.5 (Figure 2(a)).The HMM-FEM solution (star) compared with fine-scale solution (solid line) along the cross-sectionx
1
=
0.5 for periodic case (a) and random case (isotropic) (b), respectively.
(a)
(b)


## 5.2. van Genuchten-Mualem Model with Random Coefficients

For the random model, we only consider the isotropic heterogeneities case. We generate the random log-normal permeability fieldsK
s
(
x
) and α
(
x
) by the same methods used in [16]. The corresponding normal distribution of log
⁡
(
K
s
) is N
(
-
3,0.8
) and the corresponding normal distribution of log
⁡
(
α
s
) is N
(
-
2,0.8
). The correlation lengths of both k
s and α
s are l
1
=
l
3
=
0.01 in x
1 and x
3 directions. The solution of HMM-FEM is obtained on a macroscale grid 32
×
32. For each local sample cell I
l, we choose l
=
1
/
64 and solve the microproblems (25) on a grid with size of δ
/
16. We compare the HMM-FEM solution with the fine-scale solution which is obtained by solving (59) on a grid 1024
×
1024 by FVM. We compare the HMM-FEM solution with the fine-scale solution along the cross-section x
1
=
0.5 (Figure 2(b)).



---

*Source: 657816-2014-03-20.xml*