Finite Volume Methods G. Manzini 1 1 IMATI - CNR, Pavia Padova - - PowerPoint PPT Presentation

Finite Volume Methods

• G. Manzini1

1IMATI - CNR, Pavia

Padova 16/4/2007

Outline

Contents:

• part 1: an informal introduction to Finite Volume Methods

(FVMs);

• part 2: application to steady transport problems in

convection-dominated regimes and layer calculation;

Other collaborations:

• application to strongly anisotropic diffusion models: the locking effect

(joint work with M. Putti);

• application to partially saturated multi-layered soils: mass conservation

(joint work with. S. Ferraris);

The basic model problem: the convection-diffusion equation

• Find a function u(x, t) defined on Ω ×

∂u ∂t + div (βu − K∇u) = G, in Ω ×

with Ω ⊂

• u, conserved quantity;
• β, velocity field, K, diffusion tensor;
• f, production/consumption source term;

The basic model problem: the convection-diffusion equation

• u0(x) initial condition,

u = u0 in Ω, t=0,

• gD, gN, boundary condition data defined on ∂Ω = ΓD ∪ ΓN

u = gD,

• n

ΓD ×

n · K∇u = gN,

• n

ΓN ×

with Ω ⊂

Sub-problems

∂u ∂t + div (βu − K∇u) = G, in Ω ×

u = u0 in Ω, t=0, u = gD,

• n

ΓD ×

n · K∇u = gN,

• n

ΓN ×

• ∂u/∂t = 0, steady-state case;
• |β|>

>K, advection-dominated case, layers;

• K = diag(1, δ), with δ<

<1, or δ> >1 (and possibly |β| = 0), strongly anisotropic diffusion tensors.

• ΓN (ΓD) = ∅, pure Dirichlet (Neumann) problem.

Four steps towards finite volume discretizations:

• introduce a mesh partitioning of the domain Th =
• T
• , T,

control volumes;

• re-formulate the equation by integrating over each control

volume T as a balance of physical fluxes;

• mimic the conservation property of the underlying

divergence operator by a discrete balance of properly defined numerical fluxes;

• solve for an approximation of cell averages of the

conserved quantity and eventually reconstruct piece-wise linear approximate solutions.

Deeper inside. . .

• we reformulate the model problem as an integral equation on

each control volume;

• we represent the integral equation for each control volume as an
• rdinary differential equation (ODE);
• collecting ODEs from all control volumes we get a system of

(possibly non-linear) algebraic equations. . .

• . . . that are (hopefully) amenable to a solution using

computational methods.

Thus, we need to approximate volume integrals and surface integrals to form algebraic expressions.

Domain, Zone, Grid, and Cell

Prior to discussing the Finite Volume approximation, let us examine the control volumes on which volume and surface integrals will be

• approximated. . .

The control volumes exists at several levels:

• flow domain, extent of CFD analysis
• zone, divide domain for convenience, if needed
• grid, divides each zone into cells
• cell, smallest control volume, but finite

Domain, Zone, Grid, and Cell

domain zone 1 zone 2 grid 2, 32 × 16 cells grid 1, 16 × 8 cells

General Classification of Grids

Grids can be classified as

• structured,
• unstructured,
• hybrid.
• Structured grids use a structured grid topology:
• cells are arranged in an array structure
• location of neighbor items (cells, faces, vertices) is implicit

in the array indices (i,j,k)

Use of structured grids allows efficient storage and book-keeping of grid items information.

Example of a structured grid

i,j+1

• i,j

i+1,j i+1,j+1 i,j+1

i−1,j i,j i+1,j i,j−1

General Classification of Grids

• unstructured grids do not have an underlying “regular”

grid topology:

• no inherent ordering of grid items,
• the arrangement of grid items must be specified explicitly

(by using lists).

Use of unstructured grids allows greater flexibility in generating and adapting grids at the expense of greater storage of cell information.

• hybrid grids contain both structured and unstructured

grids, usually in separate zones

Example of an unstructured grid

14

• 2

27 38 51

37 24 41 6

• Control volume T14 is formed by vertices labeled by {37, 24, 41, 6}

and has neighbours labeled by {2, 51, 38, 27}.

Remarks

Grids of quadrilaterals may be considered as structured or unstructured, but when control volumes of more general shape are considered, we must use unstructured grids:

Anatomy of a finite volume cell

As integral equation are approximated on the finite volume cell to form an algebraic relation, we point out that:

• any cell can take on a generalized shape, but they are

usually

• non-overlapping;
• convex shaped (not always required);
• any cell contains a finite (positive) volume;
• the size of cells, usually denoted by the symbol h, indicates

the level of computational resolution of error analysis;

For 2-D finite volume cells:

• a volume cell is a polygon defined by a control poly-line;
• the control poly-line is subdivided into a finite number of

edges (2-D);

• edges are usually straight lines;
• an edge is bounded by two non-coincident vertices.

For 3-D finite volume cells:

• a volume cell is a polyhedron defined by a control surface;
• the control surface is faceted into a finite number of faces

(3-D);

• faces can take on a variety of shapes;
• a face is bounded by edges;
• edges are usually straight lines.

Integral conservation formulation

• We take K = ν

control volume T,

• T

∂u ∂t dV +

• T

div (βu − ν∇u) dV =

• T

G dV,

• assume that control-volume shapes be time-independent,

d dt

• T

u dV +

• T

div (βu − ν∇u) dV =

• T

G dV,

• and apply the Gauss divergence theorem,

d dt

• T

u dV +

• ∂T

nT· (βu − ν∇u) dS =

• T

G dV.

Integral conservation formulation

• Some remarks on the equation structure:

d dt

• T

u dV

• time variation of

u in volume T

+

• ∂T

nT · (βu − ν∇u) dS

• flux balance of u through

surface ∂T

=

• T

G dV

source of u in volume T

.

• We divide by |T| and split face contributions to fluxes

d dt 1 |T|

• T

u dV + 1 |T|

• e∈∂T

|e|

• 1

|e|

• e

ne · βu dS − 1 |e|

• e

ne · ν∇u dS

• = 1

|T|

• T

G(t) dV

• to obtain

d dt uT(t) + 1 |T|

• e∈∂T

|e|

• Fe(u) − Ge(u)
• = GT,

Summarizing:

• finite volume methods search for an approximation of

d dt uT(t) + 1 |T|

• e∈∂T

|e|

• Fe(u) − Ge(u)
• = GT(t),

for all T ∈ Th

where

uT = 1 |T|

• T

u dV, cell average of u over T Fe = 1 |e|

• e

ne,T·βu dS, advective flux Ge = 1 |e|

• e

ne,T·ν∇u dS, diffusive flux GT = 1 |T|

• T

G dV, cell average of source term over T

Some general remarks

d dt uT(t) + 1 |T|

• e∈∂T

|e|

• Fe(u) − Ge(u)
• = GT,

for all T ∈ Th

• uT is the conserved quantity representing the flow;
• the integral conservation equation applies for each control

volume;

• control surface bounds the control volume;
• flow is driven through the control surface by fluxes Fe and

Ge;

• flow can be time-varying (unsteady), i.e. Fe and Ge depend
• n uT(t);
• further, we may consider time-dependent control volumes

and control surfaces as well.

Finite Volume Approximation

In the integral equation

d dt uT(t) + 1 |T|

• e∈∂T

|e|

• Fe(u) − Ge(u)
• = GT,

for all T ∈ Th

we approximate the continuous flux balance summation term by the discrete flux balance

• e∈∂T

|e|

• Fe(u) − Ge(u)
• ≈
• e∈∂T

|e|

• Fe(uh) − Ge(uh)
• ,

where

• uh ≈ {uT} collects the approximation of cell averages of u;
• Fe(uh) ≈

1 |e|

• e

ne,T·βu dS, numerical advective flux;

• Ge(uh) ≈

1 |e|

• e

ne,T·ν∇u dS, numerical diffusive flux.

Underlying assumptions on solution approximation

Note that within the finite volume framework we are approximating uh(t)|Ti = ui(t) ≈ uTi(t) = 1 |Ti|

• Ti

u dV, which is a constant (in space) quantity over Ti. Thus,

• uh = {ui} is a piecewise constant function over Ω;
• uh = {ui} is (also) a vector in

forming Th;

• the position of the solution point in the cell is not yet defined

(compare with FEM or FDM methods). Two possible interpretations:

• we are approximating cell averages;
• we are approximating the value of u at some point inside T

(barycenters work up to second-order of accuracy).

Location of the solution in the finite volume cell

The location of the flow solution and geometry of the finite volume cell with respect to the grid can be of two types:

• cell-centered cell: the flow solution is located at the

centroid of the cell volume defined by the grid lines (primary grid);

• cell-vertex cell (also node-centered cell):
• the flow solution is located at the vertices of the grid;
• the finite-volume cell is formed about the vertex (dual grid).

Each approach has its advantages and disadvantages, but if things are done right, both approaches do well.

Cell-centered grid

Vertex-centered grid

Underlying assumptions on flux approximation

Note that within the finite-volume framework we are approximating Fe(uh) ≈ 1 |e|

• e

ne,T·βu dS, numerical advective flux; Ge(uh) ≈ 1 |e|

• e

ne,T·ν∇u dS, numerical diffusive flux. Thus,

• the flux is uniform over the surface of each face of the cell;

Remark also that:

• computing the flux on the face may be one of the most difficult

and computationally intensive operations of a finite volume code.

Thus, how do we compute numerical fluxes?

Numerical convective flux: first-order scheme

• Ti

Tj ui uj

• We want to approximate numerically this integral

Feij(uh) ≈ 1 |eij|

• eij

n·βu dS by taking into account that:

• the face eij is shared by the control volumes Ti and Tj;
• available information is represented by ui and uj;
• information is travelling transported by the velocity field β

⇒ upwinding!

Numerical convective flux: first-order scheme

• Ti

Tj ui uj eij

Feij(uh) = β+

ij ui + β− ij uj ≈

1 |eij|

• eij

n·βu dS upwind velocities: β±

ij =

• βij ± |βij|
• /2

with βij = 1 |eij|

• eij

β · nij dS. This approximation is exact for constant functions!

Numerical diffusive flux

• uij

uji hij uα uβ ui uj Ti Tj

• We want to approximate numerically this integral

Gij(uh) ≈ 1 |eij|

• eij

n · ν∇u dS. so that the flux approximation is exact for linear functions, i.e. constant gradients.

Numerical diffusive flux

• uij

uji hij uα uβ ui uj Ti Tj

Gij(uh) ≈ 1 |eij|

• eij

n · ν∇u dS.

• if G⋄

ij(uh) is a constant approximation of ∇u on the diamond Dij

defined by xi, xj, xα and xβ, we set Gij(uh) = nij · νijG⋄

ij(uh).

• we must express G⋄

ij(uh) as function of uh = {uT}.

Derivation of gradient formula on diamonds, (i)

• uij

uji hij uα uβ ui uj Ti Tj

The flux approximation must be exact for linear functions:

• for any integrable function u, we define G⋄

ij(u) as the gradient of

u averaged on the “diamond” Dij: G⋄

ij(u) ≡ ∇uDij =

1 |Dij|

• Dij

∇u dV

• transform the volume integral into a boundary integral by

applying the Gauss-Green theorem: 1 |Dij|

• Dij

∇u dV = 1 |Dij|

• ∂Dij

ηu dS.

Derivation of gradient formula on diamonds, (ii)

• we split the boundary integral in edge contributions,
• ∂Dij

=

• eiα

+

• eαj

+

• ejβ

+

• eβi

= (⋆);

• we apply the trapezoidal rule, which is exact for linear functions,

(⋆) = eiαηiα ui + uα 2 + eαjηαj uα + uj 2 +ejβηjβ uj + uβ 2 + eβiηβi uβ + ui 2 ;

• using ηiα+ηαj +ηjβ+ηβi =0, ηiα+ηβi+ηij =0, and similar

geometrical relations yields (⋆) = uj − ui |xi − xj| ηαβ ηαβ·τij + uβ − uα |xβ − xα| ηij ηij ·ταβ = G⋄

ij,

(formula by Coudière, Vila, Villedieu, M2AN, (1999)).

Derivation of gradient formula on diamonds, (iii)

• uij

uji hij uα uβ ui uj Ti Tj

Bertolazzi-M., M3AS, (2004) proposed an alternative derivation that yields the following (equivalent) definition.

• We start by defining the one-sided edge gradients:

Gij(uh) = uij − ui hij nij + uβ − uα |eij| tij, Gji(uh) = uji − uj hji nji + uα − uβ |eij| tji

Derivation of gradient formula on diamonds, (iv)

• uij

uji hij uα uβ ui uj Ti Tj

• Edge gradients are averaged to obtain

G⋄

ij(uh) = WijGij(uh) + WjiGji(uh),

Wij = |Ti| |Ti| + |Tj|.

• Note that

Wij = |Ti| |Ti| + |Tj| = hij Hij where Hij = hij + hji,

Derivation of gradient formula on diamonds, (v)

• uij

uji hij uα uβ ui uj Ti Tj

Finally, some calculations yields the final formula for diamond gradients: G⋄

ij(uh) =

uj − ui Hij + uij − uji Hij

• nij + uβ − uα

|eij| tij.

The gradient formula on diamonds

Coudière, Vila, Villedieu, M2AN, (1999) : G⋄

ij(uh) = uj − ui

|xi − xj| ηαβ ηαβ·τij + uβ − uα |xβ − xα| ηij ηij ·ταβ ; Bertolazzi & M., M3AS, (2004) : G⋄

ij(uh) = Wij

uij − ui hij nij + Wji uji − uj hji nji + uβ − uα |eij| tij. It can be noted that:

• for Wij = hij/Hij the two gradient formulas coincide. . .
• . . . but a different choice of edge weights gives a different

gradient formula, and thus a different numerical flux;

• furthermore, non-linear weights can be easily introduced;
• the formulation by Bertolazzi-M., 2004 can be readily extended

to n − D problems with n > 2.

Gradient formula in 3-D

• uji

uij hij hji uα uβ uγ ui uj

Ti Tj

G⋄

ij(uh) = Wij

uij − ui hij nij + Wji uji − uj hji nji

• normal terms

+uβ − uα |eij| t(αβ)

ij

+ uγ − uα |eij| t(αγ)

ij

• tangential terms

.

Mumble Mumble. . .

So, we are very happy with these gradient formulas. . . . . . but there is a problem. . . . . . these formulas for edge gradient require the knowledge

• f vertex values!!!

Reconstruction of vertex values by Least-Squares

v T

• u(x) = av+bv·(x −xv),

x ∈∪T∈σvT The reconstruction weights wv,T uv = u(xv) = av =

• T wv,TuT,

v internal vertex, g(xv), v boundary vertex. are given by taking (av, bt

v) as minimizers of

J (av, bt

v) =

• T
• u(xT) − uT

2 .

Reconstruction of vertex values by Least-Squares

Linear Consistency

• T

wv,T = 1,

• T

wv,T(xT − xv) = 0. Uniform Boundedness max

T

• wv,T
• < C.

Using these two properties, a weaker form of consistency for the original scheme has been shown (Coudière et al. (1999)). Alternative choice of weights for vertices are possible. . .

Numerical convective flux: second-order scheme

• Gi

Gj Ti Tj b uij b uji ui uj

We want an approximation of the convective flux integral that is exact for linear functions: Feij(uh) ≈ 1 |eij|

• eij

n·βu dS by using the cell reconstructions

• uij = ui +(xij −xi)·Gi(uh)
• uji = uj +(xij −xj)·Gj(uh)

and upwind velocities β±

ij =

• βij ± |βij|
• /2

with βij = 1 |eij|

• eij

β · nij dS.

Numerical convective flux: second-order scheme

• Gi

Gj Ti Tj b uij b uji ui uj b u|Ti (x)=ui +(x −xi)·Gi(uh) b uij = b u|Ti (xij) b u|Tj (x)=uj +(x −xj)·Gj(uh) b uji = b u|Tj (xij)

Feij(uh) = β+

ij

uij + β−

ij

uji ≈ 1

• eij
• eij

n·βu dS Feij(uh) = β+

ij ui + β− ij uj

• first-order term

+ β+

ij (xij − xi) · Gi(uh) + β− ij (xij − xj) · Gj(uh)

• second-order correction

,

Reconstruction of cell gradients and slope limiters

To complete the scheme, here is how we define the cell gradients Gi and Gj: (i) define a constant gradient GT(uh) within the cell T:

• Triangles: take the gradient of the linear function going

through the vertex values {(xv, uv)};

• Convex polygons: take the gradient of the linear function

approximating the vertex values {(xv, uv)} in the Least-Square sense; (ii) reconstruct the solution within the cell T as ˜ uT(x) = uT + (x − xT) · GT(uh) (iii) introduce a slope limiter factor φT(uh) to reduce slopes GT(uh) = φT(uh) GT(uh) (à la Barth-Jespersen) to obtain

• uT(x) = uT + (x − xT) · GT(uh),

x ∈ T.

“Diamond” scheme with non-linear upwind

The shock-capturing technique using limiters introduces a non-linear character into the scheme.

• Some theoretical results:
• Well-posedness: there exists at least one solution of the

non-linear discrete problem

• Quasi-uniqueness: the numerical solutions are within a ball
• f diameter decreasing as h2

“Diamond” scheme with non-linear upwind

• a very important implementation detail:
• A0 + A1(uh)
• uh = f

→ A0uk+1

h

= f − A1(uk

h)uk h

solving by a direct method (ma41, HSL) asks for

• a single factorization of the “first-order” matrix A0
• several evaluations of the non-linear correction A1(uk

h)uk h

• and corresponding resolution of two triangular systems

using the factors of A0. ⇒ non-linear iterations have a very low cost!

Stencils: reconstruction of vertex values

uk uV

• uk
• ≡ set of cell-averaged data needed to reconstruct the

vertex value uV

Stencils: numerical definition of the edge flux

uk eij

• uk
• ≡ set of cell-averaged data needed to evaluate/approximate

the numerical flux integrals Fij, Gij across the edge eij.

Stencils: cell balance of finite volume fluxes

The assembly process yields the system of (non-)linear equations

• F
• adv.

− G

• diff.
• u(l+1) =

G

• source

− f(u(l))

nlin adv.

+ g

• diff

ui uk Sparsity of F − G: row i − → flux balance of Ti

• ui
• ∪
• uk
• k −

→ non-zeroes of row i

Limiter #1

ui uk ˜ uV

• We reduce the slope of the reconstructed gradients as follows

Gi(uh) → ΦGi(uh), where Φ is the maximum value in [0, 1] such that min{uk} ≤ min{˜ uV} ≤ max{˜ uV} ≤ max{uk} and ˜ uV = ui + (xV − xi) · ΦGi.

Limiter #2

ui uj ˜ uij

• Additionally, we may require that

|˜ uij − ui| ≤ |uj − ui| & sign ˜ uij − ui

• = sign
• uj − ui
• ,

where ˜ uij = ui + (xij − xi) · ΦGi, and xij is the mid-point of eij.

Upwind-biased Least-Squares vertex reconstruction

v T

Let us consider the linear solution repre- sentation:

• u(x) = av+bv·(x −xv),

x ∈∪T∈σvT

We solve the ΛT(β)-weighted Least Squares problem for (av, bt

v) by

minimizing J (av, bt

v) =

• T

ΛT(β)

• u(xT) − uT

2 . The weights ΛT(β) are chosen to incorporate upwind information into the vertex values.

v T

Recall: J (av, bt

v) =

• T

ΛT(β)

• u(xT) − uT

2 .

With upwind-biased Least-Squares coeffiicients ΛT(β), the reconstruction weights wv,T still read as uv = u(xv) = av =

• T wv,TuT,

v internal vertex, g(xv), v boundary vertex. but now upwind information is incorporated into the vertex values.

Final Remarks

Finite Volume Methods are based on

• the conservative formulation of the model equations;
• a discrete form of flux balance for any control volume;
• the numerical discretization of flux integrals;
• some reconstruction algorithm for getting vertex values

from cell averages.