Numerical Simulation of Transport Processes in Porous Media

Numerical Simulation of Transport Processes in Porous Media

Finite-Volume Methods

Olaf Ippisch

Interdisziplin¨ ares Zentrum f¨ ur Wissenschaftliches Rechnen Universit¨ at Heidelberg

Im Neuenheimer Feld 368 D-69120 Heidelberg Telefon: 06221/54-8252

E-Mail: olaf.ippisch@iwr.uni-heidelberg.de

November 3, 2009

Cell-Centred Finite-Volume Method

We want to discretise the steady-state ground-water equation

∇ · ~ J w (~ x) + r w (~ x) = 0 with

J w = −K s (~ x )∇p w

with the Cell-Centred Finite-Volume method.

Divide Domain into Grid Cells

First we divide grid into rectangular grid cells g ij

x_0,0 x_4,0

x_4,0.5

N=5

x_1,0 x_2,0 x_3,0 x_0,1

x_0,2 x_0,3 x_0,4

x_3.5,1 x_4.5,1

x_4,1.5

Transformation of Volume Integral into Boundary Integral

We demand that the integral of the partial differential equation over each grid cell is fulfilled:

Z

g

∇ · ~ J _w dx dy = Z

g

r(~ x) dx dy

and use the Satz of Gauss to transform the volume integral over the divergence of the flux into a boundary integral over the flux normal to the boundary:

⇔

|{z}

Satz of Gauss

Z

∂g

~ J w · ~ n ds = Z

g

r (~ x) dx dy

Inner Grid Cell

(

)

(

) (

)

(

)

(i,j−1) (i,j)

(i−1,j) (i−0.5,j) (i+0.5,j) (i+1,j)

(i,j+0.5)

(i,j−0.5)

i+0.5,j

F (i,j+1)

Finite Volume Discretisation: Split into Sum over Faces

For our rectangular cell, we can split the integral over the boundary of the cell into integrals over each face

Z

∂g

~ J _w · ~ n ds = X

k=i±0.5

Z

F

~ J _w · ~ n ds + X

l=j±0.5

Z

F

~ J _w · ~ n ds

and approximate the integral over each face with the Midpoint rule

≈

|{z}

Midpointrule

X

k =i±0.5

~ J _w (~ x _k,j ) · ~ n · h

|{z}

Face Area

+ X

l=j±0.5

~ J _w (~ x _i,l ) · ~ n · h

|{z}

Face Area

If the permeability is a diagonal matrix the flux over a face depends only on the

gradient in the normal direction

Finite Volume Discretisation: Insert Flux law

The multiplication with the (normalised) normal vector only influences the sign of the flux integral

X

k=i±0.5

~ J w (~ x k,j ) · ~ n · h + X

l=j±0.5

~ J w (~ x i,l ) · ~ n · h =

−K xx (~ x _i−0.5,j ) · ∂p

∂x (~ x _i−0.5,j ) · (−1)

| {z }

from n

·h

−K xx (~ x i+0.5,j ) · ∂p

∂x (~ x i+0.5,j ) · (1)

|{z}

from n

·h

−K yy (~ x _i,j−0.5 ) · ∂p

∂x (~ x _i,j−0.5 ) · (−1)

| {z }

from n

·h

−K yy (~ x i,j+0.5 ) · ∂p

∂x (~ x i,j+0.5 ) · (1)

|{z}

from n

·h

Finite Volume Discretisation: Approximate Derivatives

The gradient at the face midpoint is approximated by a central difference quotient

≈

|{z}

approx.Derivative

+K xx (~ x _i−0.5,j ) · p(~ x i,j ) − p(~ x _i−1,j )

h · h

−K _xx (~ x _i+0.5,j ) · p(~ x i+1,j ) − p(~ x i,j )

h · h

+K yy (~ x _i,j−0.5 ) · p(~ x i,j ) − p(~ x i,j−1 )

h · h

−K _yy (~ x _i,j+0.5 ) · p(~ x _i,j+1 ) − p(~ x _i,j )

h · h

Finite Volume Discretisation: Approximate Derivatives

The integration of the source/sink term is also done with the midpoint rule:

Z

g

r (~ x)dx dy ≈ h ² r (~ x i,j )

Matrix Contribution of each Grid Cell

We get one line of a linear equation system for each grid cell:

−K

(~ x

) · p

− K

(~ x

) · p

−K

(~ x

) · p

− K

(~ x

) · p

+ [K

(~ x

) + K

(~ x

) + K

(~ x

) + K

(~ x

)] · p

= h

r(~ x

)

Dirichlet Boundary Conditions

Let us assume that at x = 0 there is a Dirichlet boundary:

h/2 x_0,j

h

The derivative between the face midpoint and the element midpoint can be approximated by a difference quotient (only first order):

∂p

∂x (~ x _−0.5,j ) ≈ p(~ x _0,j ) − p _d (0, y _j )

h/2

Matrix Contribution at Dirichlet Boundary x = 0

The constant term −K xx (~ x _i−0.5,j ) · p d (0, y j ) is brought to the right-hand side of the equation:

−K

(~ x

) · p

−K

(~ x

) · p

− K

(~ x

) · p

+ [2K

(~ x

) + K

(~ x

)

+K

(~ x

) + K

(~ x

)] · p

= h

r(~ x

) + 2K

(~ x

) · p

(0, y

)

Neumann Boundary Conditions

To integrate Neumann boundary conditions we go back to the point before the integration of the face fluxes with the midpoint rule. For each face we had to determine

Z

F

~ J w · ~ n ds

At a Neumann boundary ~ J w · ~ n is given directly by the boundary condition φ n (~ x ), we can therefore use

Z

F

~ J w · ~ n ds ≈

|{z}

Midpointrule

h · φ ~ N (~ x ) · ~ n

at each Neumann boundary.

Matrix Contribution at Neumann Boundary x = 0

We transfer the constant term h · φ N (~ x −0.5,j ) to the right hand side:

−K

(~ x

) · p

−K

(~ x

) · p

− K

(~ x

) · p

+ [K

(~ x

) + K

(~ x

) + K

(~ x

)] · p

= h

r(~ x

) − h · φ

(~ x

)

Different Grid Spacing in x and y direction

If the grid spacing in x and y direction is different, the h factors do not clear:

− h

h

(K

(~ x

) · p

− K

(~ x

) · p

)

− h

h

(K

(~ x

) · p

− K

(~ x

) · p

) +

» h

h

(K

(~ x

) + K

(~ x

)) + h

h

(K

(~ x

) + K

(~ x

)) –

· p

= h

h

r(~ x

)

Example: 3 × 3 Grid

Let us perform a Finite-Volume discretisation of the stead-state groundwater equation on a 3 × 3 grid with a homogeneous permeability field and Dirichlet boundary condition on the north and south side and no-flux boundary conditions at the left and right:

3 5

6 7 8

noflux 4 noflux

dirichlet north

K (~ x ) =

K 0

0 K

Example: 3 × 3 Grid

0 1 2

3 5

6 7 8

noflux

4

dirichlet north

dirichlet south

The resulting linear equation system is:

0 B B B B B B B B B B

@

4K −K 0 −K 0 0 0 0 0

−K 5K −K 0 −K 0 0 0 0

0 −K 4K 0 0 −K 0 0 0

−K 0 0 3K −K 0 −K 0 0

0 −K 0 −K 4K −K 0 −K 0

0 0 −K 0 −K 3K 0 0 −K

0 0 0 −K 0 0 4K −K 0

0 0 0 0 −K 0 −K 5K −K

0 0 0 0 0 −K 0 −K 4K

1 C C C C C C C C C C A

0 B B B B B B B B B B

@ p

p

p

p

p

p

p

p

p

1 C C C C C C C C C C A

= 0 B B B B B B B B B B B

@

2Kp

2Kp

2Kp

0 0 0 2Kp

2Kp

2Kp

north

1

C

C

C

C

C

C

C

C

C

C

C

A

Effective Permeability

We assume that the permeability is a diagonal Tensor, which is depending on the position, but constant on each grid cell g ij .

We need to evaluate K at the cell boundaries x i±0.5,j±0.5 .

What is the correct value of K if it is not homogeneous but element-wise constant?

1D-Example

The steady-state groundwater flow equation is:

dJ _w

dx = 0 in Ω = (0, `

|{z}

length

)

J w = −K (x) dp dx

p

K(x)

p

with the Dirichlet boundary conditions p(0) = p ₀

p(`) = p `

because of ^dJ _dx

= 0 in Ω ⇔ J w (x) = J 0 ∈ R this means

J 0 = −K (x) dp dx ⇔ dp

dx = − J 0

K (x)

1D-Example

By integration of both sides over the domain dp

dx = − J ₀ K (x)

⇔

`

Z

0

dp

dx dx = [p(x)] <sup>`</sup> <sub>0</sub> = p <sub>`</sub> − p ₀ = −J 0

`

Z

0

1 K (x) dx

we get the flux depending on the boundary conditions and the permeability distribution:

⇔ J 0 = − `

`

R 1 dx

· p ` − p 0

`

| {z }

1D-Example, cell-wise constant Permeability

If we divide the domain into two halves with constant permeability if K (x) =

K l x ≤ ₂ ^{`</sup> K r x > <sup>`} ₂

K _r

0 l

K l

we can perform the integration and get the effective permeability

K _eff = `

`

R

0 1 K(x ) dx

= `

` 2

1 K

+ ^` _{2 K} ¹

r

= 2

1 K

+ _K ¹

r

We therefore choose for cell-wise constant permeabilities the harmonic mean K (~ x _i±0.5,j ) = 2

1

K(~ x

) + _{K(~ x} ¹

i±1,j

)

Finite-Volume Method for tensor-product Grids

− 2h

K

(~ x

) h

+ h

· p

− 2h

K

(~ x

) h

+ h

· p

− 2h

K

(~ x

) h

+ h

· p

− 2h

K

(~ x

) h

+ h

· p

+

» 2h

K

(~ x

)

+ 2h

K

(~ x

)

Complexer Grids with Cell-Centred Finite Volumes

With the Cell-Centred Finite Volume Method it is also possible to use some kind of unstructured grids:

Nested Grids

I M J K A B

γ

x

Voronoi Grids

A B

C

D

E F G

x i

x _k

C _i

Summary Cell-Centred Finite-Volume Method

• Only the integral of the partial differential equation over each grid cell must fullfill the equation.

• Implementation of Dirichlet Boundary and Neumann Boundary conditions straight forward

• Structured and unstructured grids possible

• Dirichlet boundary conditions can easily be integrated by rearranging the equation systems and bringing them to the right side of the equation.

• Neumann boundary conditions can easily be integrated in the flux integrals

• Convergence order can differ dependent on the concrete method.

Properties of the Cell-Centred Finite-Volume Method

• Advantages:

• well suited for structured grids

• locally mass conservative

• good approximation of average permeability

• limited variety of unstructured grids possible

• limited local adaptivity possible

• cheap for simple problems

• Problems:

• Only linear convergence rate on non-equidistant grids

• grid generation can be complicated (must fulfil rather strong conditions)

The Vertex-Centred Finite-Volume Method

vi

vj bi

bj

• The unknowns are located at the edges of the elements (vertices)

• Base functions are used on each element, which are parameterised with the values at the vertices

• A secondary mesh is constructed connecting the face centres and the

barycenter of the element

The Vertex-Centred Finite-Volume Method (2)

vj vi

b_i^k xikf γikf

nikf

x_ij^k n_ij^k γijk

element e_k

vi vj

b^k_i b_j^k

element e_k xijk

n_ij^k γijk

γjkf xjkf

njkf

• Material properties are assumed to be constant for each element

• The volume integrals are calculated as a sum over the subcontrol-volumes using the midpoint rule and the material properties valid for the specific control-volume. P

i

b _i ^k · r _i ^k

• The face integrals are calculated as a sum over all subcontrol-volume faces with the midpoint rule P

ij

γ _ij ^k ~ J _ij ^k ~ n ^k _ij

• The gradient at the face centres is given by the base functions.

Properties of the Vertex-Centred Finite-Volume Method

• Advantages:

• can be used for domains with complicated shape

• well suited for unstructured grids

• local adaptivity possible

• locally mass conservative

• Problems:

• grid generation can be complicated (must often fullfill certain conditions)

• more computationally expensive for simple problems

• bad approximation of average permeability