- Home
- Documents
- Multigrid Methods for A Mixed Finite Element Method chenlong/Papers/ · Multigrid Methods for A Mixed…

Published on

30-Aug-2018View

213Download

1

Transcript

Noname manuscript No.(will be inserted by the editor)

Multigrid Methods for A Mixed Finite Element

Method of The Darcy-Forchheimer Model

Jian Huang Long Chen Hongxing Rui

Received: date / Accepted: date

Abstract An ecient multigrid method for a mixed finite element method of theDarcy-Forchheimer model is constructed in this paper. A Peaceman-Rachford typeiteration is used as a smoother to decouple the nonlinearity from the divergenceconstraint. The nonlinear equation can be solved element wise with a closed for-mulae. The linear saddle point system for the constraint is further reduced intoa symmetric positive definite system of Poisson type. Furthermore an empiricalchoice of the parameter used in the splitting is found. By comparing the numberof iterations and CPU time of dierent solvers in several numerical experiments,our method is shown to convergent with a rate independent of the mesh size andthe nonlinearity parameter and the computational cost is nearly linear.

Keywords Darcy-Forchheimer model Multigrid method Peaceman-Rachforditeration

1 Introduction

Darcys law

u = Krp,

The work of the authors Jian Huang and Hongxing Rui was supported by the National NaturalScience Foundation of China Grant No. 91330106, 11171190. Long Chen was supported by NSFGrant DMS-1418934 and in part by NIH Grant P50GM76516. The work of the author JianHuang was supported by 2014 China Scholarship Council (CSC).

J. Huang H. RuiSchool of Mathematic, Shandong University, Jinan, Shandong 250100, ChinaE-mail: yghuangjian@sina.com

L.ChenDepartment of Mathematics, University of California at Irvine, Irvine, CA 92697, USAE-mail: chenlong@math.uci.edu

H. RuiE-mail: hxrui@sdu.edu.cn

2 Jian Huang et al.

describes the linear relationship between the velocity of the creep flow and thegradient of the pressure, which is valid when the Darcy velocity is extremely small[2]. Forchheimer in [8] carried out flow experiments and pointed out that whenthe velocity of the flow is relatively high, Darcys law should be replaced by theso-called Darcy-Forchheimer (DF) equation by adding a quadratic nonlinear termto the velocity, shown as follows:

K

1u+

|u|u+rp = 0.

A theoretical derivation of the Darcy-Forchheimer equation can be found in [17].In recent years, many numerical methods of the Darcy-Forchheimer equation

have been developed. A mixed finite element method with a semi-discrete schemefor the time dependent problem was considered by Park in [12]. Pan and Ruiin [11] gave a mixed element method for the Darcy-Forchheimer problem basedon the Raviart-Thomas (RT) element or the Brezzi-Douglas-Marini (BDM) ele-ment. Rui and Pan in [15] proposed the block-centered finite dierence method forthe Darcy-Forchheimer model, which was thought of as the lowest-order Raviart-Thomas mixed element with proper quadrature formula. Rui, Zhao and Pan in [16]presented a block-centered finite dierence method for Darcy-Forchheimer modelwith variable Forchheimer number. Wang and Rui in [19] constructed a stabilizedCrouzeix-Raviart (CR) element for the Darcy-Forchheimer equation. Rui and Liuin [14] introduced a two-grid block-centered finite dierence method for the Darcy-Forchheimer model.

Girault and Wheeler in [9] proved the existence and uniqueness of the solutionof the Darcy-Forchheimer model. They approximated the velocity and the pressureby piecewise constant and nonconforming CR elements, respectively. They alsosuggested an alternating directions iterative method to solve the nonlinear system.The convergences of both the iterative algorithm and the mixed finite elementscheme are also proved. Lopez, Molina and Salas in [10] carried out numericaltests of the methods studied in [9], and compared with the Newtons method forsolving this problem.

Since it is a nonlinear system, an iterative method should be used, whichcould be very computationally expensive if the convergence of nonlinear iterationis slow. Multigrid method is one of the most ecient methods on solving thenonlinear elliptic systems. It should be clarified that we no longer have a simplelinear residual equation, which is the most significant dierence between linear andnonlinear systems. The multigrid scheme we used here is the most commonly usednonlinear version of multigrid. It is called the full approximation scheme (FAS) [3]because the problem in the coarse grid is solved for the full approximation ratherthan the correction; see Section 5 for details.

We use piecewise constant and piecewise linear polynomial to discretize thevelocity and the pressure, respectively. We shall apply FAS to construct an ef-ficient V-cycle multigrid method for the nonlinear Darcy-Forchheimer equationand demonstrate the eciency of our multigrid method. We use the Peaceman-Rachford (PR) type iteration developed in [9] as a smoother. The most relevantwork is [10] and our improvement are: 1. We choose a smaller value of the stop-ping criterion by achieving a better approximation of the pressure accuracy. 2. Wereport a better choice of the parameter for decoupling the nonliearity from theconstraint rather than the suggested value = 1 for dierent values of , and

Multigrid Method of The Darcy-Forchheimer Model 3

show the advantage of our choice by comparing the number of iterations and CPUtime. 3. We carry out some experiments to show the eciency of our multigridsolver by comparing with the PR iterative solvers. 4. We reduce the saddle pointsystem into a Poisson system and demonstrate the eciency of our approach.

The remainder of this article is organized as follows: The model problem isdemonstrated in Section 2. The mixed weak formulation and the discrete weakformulation are presented in Section 3. The PR iteration and an ecient solverfor the linear saddle point systems are posted in Section 4. We construct a V-cyclemultigrid scheme by applying FAS for the nonlinear problem in Section 5. Somenumerical experiments using our multigrid method are carried out in Section 6 toverify that the eciency of our method in comparison with solving this nonliearproblem using the other iterative methods. Finally, conclusions, and further ideasare presented in Section 7.

2 The Problem and Some Notations

We consider the steady Darcy-Forchheimer flow of a single phase fluid in a porousmedium in a two-dimensional bounded domain, with Lipschitz continuous bound-ary @:

K

1u+

|u|u+rp = f in , (2.1)

with the divergence constraint

divu = g in , (2.2)

and Neumann boundary conditions,

u n = gN on @, (2.3)

where u and p are the velocity vector and the pressure, respectively; , and aregiven positive constants that represent the viscosity of the fluid, its density and itsdynamic viscosity, respectively; || denotes the Euclidean vector norm |u|2 = u u,n is the unit exterior normal vector to the boundary of the given domain ; K isthe permeability tensor, assumed to be uniformly positive definite and bounded.According to the divergence theorem, g and gN are given given functions satisfyingthe compatibility condition

Z

g (x)dx =

Z

@gN ()d. (2.4)

We use the standard notations of the Sobolev spaces and the associated norms,see e.g.[1]. We also use the space

L20 () =

v 2 L2 () :Z

v (x) dx = 0

.

4 Jian Huang et al.

3 Weak Formulation

In this section, we define the function spaces as follows:

X = L3()2,

M = W 1,32 () \ L20 () .

where the zero mean value condition is added because p is only defined by the (2.1)-

(2.3) up to an additive constant. Given f 2 L3(), g 2 L32 (@), the variational

formulation is: find a pair (u, p) in X M such that

Z

K

1u

' dx+

Z

|u| (u ') dx

+

Z

rp ' dx =

Z

f ' dx, 8' 2 X,

(3.1)

Z

rq u dx =

Z

gq dx+

Z

@gNq dx, 8q 2M. (3.2)

The variational formulation (3.1), (3.2) and the original problem (2.1)-(2.3)are equivalent by using the Greens formula:

Z

v rq dx =

Z

q div v dx+ hq, v ni@ , 8q 2M, 8v 2 H, (3.3)

whereH =

n

v 2 L3()2 : div v 2 L65 ()

o

.

If the given functions g and gN satisfy the compatibility condition (2.4), thenthe problem has a unique solution (u, p) in X M [9].

Let be a polygon in two dimensions which can be completely triangulated bytriangles. Let T1 be a triangulation of , and the triangulations Tk (k = 2, 3, . . .) beobtained form T1 via regular subdivision, i.e. edge midpoints in Tk1 are connectedby new edges to form Tk. Therefore, Tk is a family of conforming triangulations of,

=[

T2Tk

T for k = 1, 2, 3, . . . ,

The family Tk is regular (also called non-degenerate) in the sense of Ciarlet [6].We discretize u and p in dierent finite element spaces. The velocity u is

approximated in the following space:

Xk =n

v 2 L2()2 : 8T 2 Tk, v|T 2 P20o

, (3.4)

and the pressure p is approximated in the following space:

Mk = Qk \ L20 () , (3.5)

where Pm denotes the space of polynomials of degree m, and Qk is the linear finiteelement space.

Qk =n

q 2 C0

: 8T 2 Tk, q|T 2 P1o

.

Multigrid Method of The Darcy-Forchheimer Model 5

With these spaces, we can have the k-th level discrete formulation of the prob-lem (3.1),(3.2):

Z

K

1uk

'k dx+

Z

|uk| (uk 'k) dx

+X

T2Tk

Z

Trpk 'k dx =

Z

f 'k dx, 8'k 2 Xk,

(3.6)

X

T2Tk

Z

Trqk uk dx =

Z

gqk dx+

Z

@gNqk dx, 8qk 2Mk. (3.7)

In universal practice,

hk1 = 2hk, for k = 2, 3, . . . .

Note that Tk are nested meshes, and thus

Xk1 Xk,Mk1 Mk.

In [18], the authors demonstrated that the discrete problem has a unique solution.Moreover, if Tk is shape regular with mesh size h and the solution u belongsto W 1,4() and p belongs to W 2,

32 (), then the following error estimations are

obtained in [18](Theorem 4.10):

ku uhkL2() Ch|u|W 1,4(), (3.8)

kr (p ph)kL 32 (T ) Ch

|p|W 2,

32 ()

+ kukW 1,4()

. (3.9)

4 Nonlinear Iteration

In this section, we present the Peaceman-Rachford (PR) type method developedin [9] to decouple the nonlinearity and the constraint.

First, choose an initial guess

u

0k, p

0k

by solving a linear Darcy equation:

Z

K

1u

0k

'k dx+X

T2Tk

Z

Trp0k 'k dx =

Z

f 'k dx, 8'k 2 Xk, (4.1)

X

T2Tk

Z

Trqk u0k dx =

Z

gqk dx+

Z

@gNqk dx, 8qk 2Mk. (4.2)

The linear Darcy system (4.1) and (4.2) can be rewritten in the matrix formas

A BBT 0

u

p

=

fd

w

, (4.3)

where A is the symmetric and positive definite matrix associated to the term

Z

K

1uk

'k dx,

6 Jian Huang et al.

B is the matrix corresponding to

X

T2Tk

Z

Trpk 'k dx,

and fd and w represent the right hand side of (4.1) and (4.2), respectively.Then, knowing

u

0k, p

0k

, construct a sequence

u

n+1k , p

n+1k

for n 0 in twosteps: Let be a positive parameter chosen to enhance the convergence.

1. A nonlinear step without constraint: knowing (unk , pnk ) compute the inter-

mediate velocity un+ 1

2k by solving the following equation:

1

Z

u

n+ 12

k unk

'kdx+

Z

u

n+ 12

k

u

n+ 12

k 'k

dx =

Z

f 'k dx

Z

K

1u

nk

'kdxX

T2Tk

Z

Trpnk 'k dx, 8'k 2 Xk.

(4.4)

2. A linear step with constraint: compute

u

n+1k , p

n+1k

with the known un+ 1

2k

1

Z

u

n+1k u

n+ 12

k

'kdx+

Z

K

1u

n+1k

'kdx+X

T2Tk

Z

Trpn+1k 'k dx

=

Z

f 'k dx

Z

u

n+ 12

k

u

n+ 12

k 'k

dx, 8'k 2 Xk,

(4.5)

X

T2Tk

Z

Trqk un+1k dx =

Z

gqk dx+

Z

@gNqk dx, 8qk 2Mk. (4.6)

A key observation in [9] is that because the test functions 'k, the solution

u

n+ 12

k , and rpnk are constant in each element T , the nonlinear step (4.4) can be

solved as follows:

u

n+ 12

T =1F

n+ 12

T (4.7)

where

F

n+ 12

T =1u

nT

K

1T u

nT rT pnk + fT ,

K

1T =

1|T |

Z

TK

1 (x) dx,

=12

+12

s

12

+ 4

F

n+ 12

T

.

In the second step, the linear system (4.5) and (4.6) can be rewritten in thefollowing matrix form:

A BBT 0

u

p

=

fn+ 12w

, (4.8)

Multigrid Method of The Darcy-Forchheimer Model 7

where A is the matrix corresponding to

1

Z

u

n+1k

'kdx+

Z

K

1u

n+1k

'kdx,

and fn+ 12 is the vector corresponding to

Z

f 'k dx+

1

Z

u

n+ 12

k

'kdx

Z

u

n+ 12

k

u

n+ 12

k 'k

dx.

In [9], the authors proved that (4.1), (4.2) and (4.5), (4.6) have a unique so-lution. The iterative method is convergent for an arbitrary choice of the initialguess

u

0k, p

0k

and an arbitrary positive . Numerically, choice of will aectthe convergence of the nonlinear iteration. We shall report a good choice of thisparameter.

We can reduce the linear saddle point system into a SPD system. Because ofA and A are symmetric positive definite operators, without loss of generality, wetake (4.8) as an example to expound an idea as follows:

Eliminate u from the first equation of (4.8), i.e.

u = A1

fn+ 12Bp

, (4.9)

and then, substituting to the second equation of (4.8), we get

Mp = b, (4.10)

where M = BTA1 B, b = BTA1 fn+ 12 w. After solving (4.10), we can get u

by solving (4.9).Since A is block-diagonal, A

1 can be formed easily. Equation (4.10) is the

linear finite element dscretization of an elliptic equation in primary formulation.Solving the SPD system (4.10) is much easier than the saddle point system (4.8)and many fast solvers are available. The equivalence between (4.9),(4.10) and (4.8)is obvious.

5 Multigrid Algorithm

In this section, we consider a generic system of nonlinear equations,

L (z) = s

where z, s 2 Rn. Suppose that v is an approximation to the exact solution z.Define the error e and the residual r:

e = z v,r = s L (v) .

We can solve this nonlinear system by using some iterative solver,

Lk (z) = s,

here k means the discrete problem on the k-th level mesh Tk.

8 Jian Huang et al.

Because of the iterative nature, multigird ideas should be eective on the non-linear problem. The multigrid scheme here we used for this nonlinear problem isthe most commonly used nonlinear version of multigrid. It is called the full approx-imation scheme (FAS) [3] because the problem in the coarse grid is solved for thefull approximation zk1 = I

k1k vk + ek1 rather than the error ek1. A V-cycle

multigrid scheme is described as follows:

Full Approximation Scheme (FAS).

Pre-smoothing: For 1 j m, relax m times with an initial guess v0 byvj = Rkvj1. The current approximation vk = vm.

Restrict the current approximation and its fine grid residual to the coarse grid:rk1 = I

k1k (sk Lk (vk)) and vk1 = I

k1k vk.

Solve the coarse grid problem: Lk1 (zk1) = Lk1 (vk1) + rk1.

Compute the coarse grid approximation to the error: ek1 = zk1 vk1.

Interpolate the error approximation up to the fine grid and correct the currentfine grid approximation: vm+1 vk + Ikk1ek1.

Post-smoothing: For m+ 2 j 2m+ 1, relax m times by vj = Rk0vj1.then we get the approximate solution v2m+1. Here m denotes the number of pre-smoothing and post-smoothing steps, Rk denotes the chosen relaxation method,and Ik1k is an intergrid transfer operator from the fine grid to the coarse grid.

We apply the PR iteration as the smoother Rk. We need to switch the orderingof the linear and nonlinear steps in the post-smoothing step in order to keep thesymmetry of the V-cycle. It is worth pointing out that although the chosen finiteelement spaces are nested, the constrained subspaces are non-nested when weinterpolated the correction of the velocity, which was obtained in the coarser space,to the finer space. Namely, if we directly interpolated the correction obtained onthe coarser grid to the finer grid, the approximation we go...