Finite Volume Methods: the steady convection-diffusion equation G. - - PowerPoint PPT Presentation

Finite Volume Methods: the steady convection-diffusion equation

• G. Manzini1

1IMATI - CNR, Pavia

Padova 16/4/2007

The reference problem

• Steady convected-dominated diffusion-transport problems are

ruled by the following equation: Find a function u(x) such that div (βu − ν∇u) = G, in Ω ∈

u = gD,

• n

∂ΩD, n · ν∇u = gN,

• n

∂ΩN.

• Challenging/interesting cases:

|β| > > ν

• r (almost equivalently)

ν → 0, the very big issue is calculation of layers (steep gradient regions)!

The one-dimensional case

Let us consider βu′ − νu′′ = 1, for x ∈ (0, 1), u(0) = u(1) = 0 for β, ν ∈ R+. From the differential equation, we expect that:

• for very small values of ν, i.e ν <

<β, the second derivative term νu′′ is negligible in a large part of the domain (0, 1) = ⇒ u grows almost linearly far from x = 1, i.e. u(x) ≈ 1/β;

• an abrupt decrease of u, the boundary layer, occurs near x = 1

because u must reach zero at this boundary point. The behavior of u for x close to 1 is governed by the second derivative.

The one-dimensional case: exact solution u(x)

0.2 0.4 0.6 0.8 1 0.2 0.4 0.6 0.8 1

u(x) = 1

β

• x − 1−eβx/ν

1−eβ/ν

• β = 1

1/ν = 1 1/ν = 10 1/ν = 100

The exponential layer calculation

1 Y 1 X This is what we want to approximate.

The exponential layer calculation

• 1

Y 1 X Sketch of a nodal approximation that we would like to see.

The exponential layer calculation

• 1

Y 1 X

• vershoot

Sketch of a solution with an overshoot: this might be still

• acceptable. . .

The exponential layer calculation

• 1

Y 1 X Sketch of a solution with oscillations: this is likely not to be acceptable as an approximation.

The reduced problem

div (βu − ν∇u) = G,

• reduced problem: take ν = 0;

the equation becomes hyperbolic: div (βu) = G or β · ∇u = G (assuming divβ = 0);

• streamlines: parameterized curves w(s) such that

β(w(s)) is the vector tangent to w at w(s); formally: β(w(s)) = d dsw(s); if β is a constant vector streamlines are straight lines.

The reduced problem

Note that: β(w(s)) = d ds w(s) d ds u(w(s)) = d dw u(w(s)) · d ds w(s) = β(w(s)) · ∇u

• n w = G.

This implies that

• the solution of the reduced problem u is the solution of an ODE;
• if G = 0, then du(w(s))/ds = 0, u is constant along the

streamline!

• if G = 0, then u along the streamline changes accordingly to the

streamline derivative and starting from its initial (inflow) value.

Consequences (from reduced problem)

• an inflow value may propagate along the streamline and

changes accordingly with the “inviscid” derivative, but takes a different value at the outflow boundary, forming an exponential layer;

• an inflow discontinuity may propagate along the

streamline, forming in the diffusive case an internal layer;

• a “parabolic layer” may form along walls that coincide with

characteristic lines, e.g. lines with β · n = 0. These layers, steep gradient regions, are difficult to be solved numerically!

Notation (Peclet Number)

• L, characteristic length scale in boundary, e.g. length of inflow

boundary

• ˆ

x = x/L in normalized domain

• |β|, normalization for the velocity (wind) β, e.g. β = |β| ˆ

β, where ˆ β = 1 In normalized variables, the advection-diffusion equation reads as divˆ

x

|β| L ν ˆ βˆ u − ∇ˆ

xˆ

u

• =

L2 ν

• ˆ

G We characterizes the relative contributions of convection and diffusion terms by the Pe ≡ |β| L ν , Peclet number,

Finite Volume Scheme: basic terminology

Main idea: mimic the conservation property of the underlying divergence operator on Th =

• T
• , suitable triangulation of Ω.

We approximate:

• ∂T

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

• T

G dV for each T ∈ Th by the discrete flux balance

• e∈∂T

|e|

• Fe(uh) − Ge(uh)
• = GT |T|

for each T ∈ Th where uh ≡ {uT} are the cell averages of u, and, for each e ∈ ∂T,

|e| Fe(uh) ≈

• e

ne,T·βu dS, |e| Ge(uh) ≈

• e

ne,T·ν∇u dS

First-order Finite Volume Scheme

(Formulation by Coudière, Vila, Villedieu, M2AN, (1999-2000)) The numerical integrator of fluxes is built by using cell averages and vertex values: cell averages ⇒ vertex values cell averages vertex values

• ⇒

face gradients ⇒ numerical diffusive flux cell-averages ⇒ numerical convective flux by upwinding

Extension to second-order accuracy: structure of the numerical scheme

(Formulation by Bertolazzi & M., M3AS, (2004))

The numerical integrator of fluxes is built by using cell averages and vertex values: cell averages ⇒ vertex values vertex values ⇒ cell gradients cell averages vertex values

• ⇒

face gradients ⇒ numerical diffusive flux cell averages cell gradients

• ⇒

numerical convective flux by non-linear upwinding (slope limiters)

The reference problem.

A solution with exponential boundary layers in the convection-dominated regime (ν = 10−4, β = (2, 3)t)

1 1.4 1 1

front view

1 1.4 1

lateral view (X-Axis)

u(x, y) =

• x − e

2(x−1) ν

• y2 − e

3(y−1) ν

• ,

(x, y) ∈ Ω =]0, 1[×]0, 1[.

What shall we try to do?

• 1

1

u ≈ xy2 Errors are measured by: barycenter error

u − uhΩ,h = h X

T

|T| ˛ ˛u(xT) − uT) ˛ ˛2i 1

2

L2-norm error

u − uhL2(Ω) = " X

T

Z

T

˛ ˛u(x) − uT(x) ˛ ˛2 dV # 1

2

Basically, we want an approximation of u that is

• second-order accurate away from the boundary layers (blue

region);

• do not show overshoots, undershoots, numerical oscillations, in

the domain region where the exponential layers develop (rose-red region).

First-order scheme

1 1.4 1 1

Exact Solution

1 1.4 1 1

First-order approximation

mumble mumble. . . it does not look too bad. . .

First-order scheme

1 1.4 1

Exact Solution

1 1.4 1

First-order approximation

• Mmm. . . is it really a good approximation?

Let us have a look at the errors. . .

Error tables on the reduced domain

Barycenters (L2-norm) NT hmax Eabs (uh) Erel (uh) Rate 400 5.00 10−2 6.098 10−3 5.778 10−2 −− 1600 2.50 10−2 3.020 10−3 2.857 10−2 1.015 6400 1.25 10−2 1.505 10−3 1.423 10−2 1.005 25600 6.25 10−3 7.507 10−4 7.099 10−3 1.003 Piecewise constant approximation (L2-norm) NT hmax Eabs (uh) Erel (uh) Rate 400 5.00 10−2 6.047 10−3 5.718 10−2 −− 1600 2.50 10−2 3.007 10−3 2.844 10−2 1.007 6400 1.25 10−2 1.501 10−3 1.420 10−2 1.002 25600 6.25 10−3 7.499 10−4 7.091 10−3 1.001 Convergence rate is first order.

Second-order scheme without limiters

1 1.4 1 1

Exact Solution

1 1.4 1 1

Second-order approximation

An overshoot is visible near outflow boundaries!

Second-order scheme without limiters

1 1.4 1

Exact Solution

1 1.4 1

Second-order approximation

So we need one more ingredient! We should try a limiter, but first. . .

Error tables on the reduced domain

Barycenters (L2-norm) NT hmax Eabs (uh) Erel (uh) Rate 400 5.00 10−2 1.607 10−4 1.523 10−3 −− 1600 2.50 10−2 4.093 10−5 3.872 10−4 1.975 6400 1.25 10−2 1.025 10−5 9.697 10−5 1.997 25600 6.25 10−3 2.555 10−6 2.416 10−5 2.005 Piecewise constant approximation (L2-norm) NT hmax Eabs (uh) Erel (uh) Rate 400 5.00 10−2 2.117 10−4 2.002 10−3 −− 1600 2.50 10−2 4.843 10−5 4.579 10−4 2.128 6400 1.25 10−2 1.206 10−5 1.140 10−4 2.005 25600 6.25 10−3 3.004 10−6 2.841 10−5 2.005 Errors are not bad!

A slope limiter driven by vertex values. . .

1 1.4 1 1

Exact Solution

1 1.4 1 1

Second-order approximation with slope limiter

The overshoot is greatly reduced!

A slope limiter driven by vertex values. . .

1 1.4 1

Exact Solution

1 1.4 1

Second-order approximation with slope limiter

• Mmm. . . it looks nice, but let us try a stronger limiter. . .

The slope limiter driven by vertex and mid-point face

• values. . .

1 1.4 1 1

Exact Solution

1 1.4 1 1

Second-order approximation with slope limiter

The overshoot is no more visible. . .

The slope limiter driven by vertex and mid-point face

• values. . .

1 1.4 1

Exact Solution

1 1.4 1

Second-order approximation with slope limiter

. . . but these nice plots are not for free. . .

Error tables on the reduced domain, limiter #1

Barycenters (L2-norm) NT hmax Eabs (uh) Erel (uh) Rate 400 5.000 10−2 9.531 10−4 9.031 10−3 −− 1600 2.500 10−2 4.119 10−4 3.897 10−3 1.212 6400 1.250 10−2 1.627 10−4 1.538 10−3 1.341 25600 6.250 10−3 6.062 10−5 5.732 10−4 1.424 Piecewise constant approximation (L2-norm) NT hmax Eabs (uh) Erel (uh) Rate 400 5.00 10−2 1.163 10−3 1.100 10−2 −− 1600 2.50 10−2 4.888 10−4 4.622 10−3 1.250 6400 1.25 10−2 1.905 10−4 1.801 10−3 1.359 25600 6.25 10−3 7.065 10−5 6.680 10−4 1.431 Accuracy is reduced!

Error tables on the reduced domain, limiter #2

Barycenters (L2-norm) NT hmax Eabs (uh) Erel (uh) Rate 400 5.00 10−2 9.560 10−4 9.059 10−3 −− 1600 2.50 10−2 4.111 10−4 3.889 10−3 1.219 6400 1.25 10−2 1.623 10−4 1.535 10−3 1.341 25600 6.25 10−3 6.042 10−5 5.713 10−4 1.425 Piecewise constant approximation (L2-norm) NT hmax Eabs (uh) Erel (uh) Rate 400 5.00 10−2 1.163 10−3 1.100 10−2 −− 1600 2.50 10−2 4.877 10−4 4.611 10−3 1.254 6400 1.25 10−2 1.902 10−4 1.798 10−3 1.358 25600 6.25 10−3 7.047 10−5 6.663 10−4 1.432 Accuracy is reduced!

Upwinding vertex reconstruction for cell-gradients. . .

1 1.4 1 1

Exact Solution

1 1.4 1 1

Second-order approximation No slope limiters!

There is no overshoot. . .

Upwinding vertex reconstruction for cell-gradients

1 1.4 1

Exact Solution

1 1.4 1

Second-order approximation No slope limiters!

And the errors?

Error tables on the reduced domain

Barycenters (L2-norm) NT hmax Eabs (uh) Erel (uh) Rate 400 5.00 10−2 1.543 10−4 1.462 10−3 −− 1600 2.50 10−2 3.993 10−5 3.777 10−4 1.952 6400 1.25 10−2 1.012 10−5 9.575 10−5 1.980 25600 6.25 10−3 2.539 10−6 2.401 10−5 1.995 Piecewise constant approximation (L2-norm) NT hmax Eabs (uh) Erel (uh) Rate 400 5.00 10−2 1.921 10−4 1.816 10−3 −− 1600 2.50 10−2 4.802 10−5 4.541 10−4 2.000 6400 1.25 10−2 1.201 10−5 1.135 10−4 1.999 25600 6.25 10−3 2.997 10−6 2.834 10−5 2.001 Wow!!!

Diamond scheme on a very general mesh

The computational grids.

• convex polygons: mainly octagons and parallelograms, but also hexagons near

the boundaries and pentagons in the corners;

• grids are built by dualization of refined criss-cross triangulations;
• nodes are remapped by a suitable coordinate transformation (Lipnikov-Shaskov).

#polys=221, #edges=700 #vertices=480 #polys=841, #edges=2600 #vertices=1760

Diamond scheme on a very general mesh

Numerical Solution (plot on 1760 vertices):

−0.1 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 1.1 −0.1 1 1.1 1 1

front view

−0.1 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 1.1 −0.1 1 1.1 1 1

rear view

Diamond scheme on a very general mesh

Behavior in the diffusive regime, ν = 1, β = (2, 3)t, Ω =]0, 1[×]0, 1[. 1st-order upwind 2nd-order upwind NT h Barycenters Rate Barycenters Rate 221 1.093 10−1 5.482 10−2 −− 2.006 10−2 −− 841 5.598 10−2 3.212 10−2 0.798 6.668 10−3 1.645 3281 2.810 10−2 1.771 10−2 0.863 1.967 10−3 1.771 12961 1.406 10−2 9.355 10−3 0.922 5.341 10−4 1.883 51521 7.033 10−3 4.815 10−3 0.958 1.391 10−4 1.941 205441 3.517 10−3 2.443 10−3 0.978 3.548 10−5 1.971

Diamond scheme on a very general mesh

Behavior in the convection-dominated regime, ν = 10−4, β = (2, 3)t, ˆ Ω =]0, 0.8[×]0, 0.8[. 1st-order upwind 2nd-order upwind NT h Barycenters Rate Barycenters Rate 221 1.093 10−1 1.915 10−1 −− 1.490 10−2 −− 841 5.598 10−2 9.351 10−2 1.071 3.599 10−3 2.122 3281 2.810 10−2 4.632 10−2 1.019 9.071 10−4 1.999 12961 1.406 10−2 2.324 10−2 0.996 2.350 10−4 1.951 51521 7.033 10−3 1.157 10−2 1.006 5.912 10−5 1.991 205441 3.517 10−3 5.760 10−3 1.006 1.473 10−5 2.005

The barycenter error is measured by u − uhˆ

Ω,h =

" X

T⊂ˆ Ω

|T| ˛ ˛u(xT) − uT) ˛ ˛2 # 1

2

Diamond scheme on a very general mesh

0.001 0.01 0.1 10

• 5

10

• 4

10

• 3

10

• 2

10

• 1

10

mesh size h

u − uhL2(Ω) / uL2(Ω)

linear SD-FEM 1st-ord UPW-FVM

ν = 1 (circles), ν = 10−1 (squares), ν = 10−4 (diamonds) ν = 10−4, ˆ Ω =]0, 0.8[×]0.8[, (triangles)

The second reference problem

−1.68 −1 1 2 2.52 1 1

r(x, y) = βyx − βxy + 1; ǫ = ν |βx| + |βy|; u(x, y) = 1 2 − arctan 1 ǫ r(x, y)

• We take β = (1, 3) and ν = 10−4

.

First-order scheme

−1.68 −1 1 2 2.52 1 1

Exact Solution

−1.68 −1 1 2 2.52 1 1

First-order approximation

Stable scheme, but very diffusive. . .

Second-order scheme without limiters

−1.68 −1 1 2 2.52 1 1

Exact Solution

−1.68 −1 1 2 2.52 1 1

Second-order approximation without limiters

Still stable, sharp transition, but. . .

Second-order with limiter driven by vertex values

−1.68 −1 1 2 2.52 1 1

Exact Solution

−1.68 −1 1 2 2.52 1 1

Second-order approximation using limiter #1

No overshoots/undershoots but quite diffusive scheme.

Second-order with limiter driven by vertex values and mid-point face values

−1.68 −1 1 2 2.52 1 1

Exact Solution

−1.68 −1 1 2 2.52 1 1

Second-order approximation using limiter #2

No overshoots/undershoots but quite diffusive.

Upwinding vertex reconstruction for cell-gradients

−1.68 −1 1 2 2.52 1 1

Exact Solution

−1.68 −1 1 2 2.52 1 1

Second-order approximation with upwind-biased gradient reconstruction and no limiters!

Gosh, I really hoped something else!

. . . and imposing a local maximum principle

−1.68 −1 1 2 2.52 1 1

Exact Solution

−1.68 −1 1 2 2.52 1 1

Second-order approximation with upwind-biased gradient reconstruction and a local maximum principle

That’s fine, but it is not better than the others.

Non-linear diffusion: structure of the numerical scheme∗

The numerical integrator of fluxes is built by using cell averages and vertex values: cell averages ⇒ vertex values vertex values ⇒ cell gradients cell averages vertex values

• ⇒

non-linear face gradients ⇒ numerical diffusive flux cell averages cell gradients

• ⇒

numerical convective flux by non-linear upwinding (slope limiters)

∗(Bertolazzi & M., SINUM, 2005)

Non-linear Diamond Scheme

• isolate the term (uj − ui) in the one-sided gradients:

nij · Gij(uh) =

• k

wij

k uk = ˜

Dij uj − ui Hij + gij(uh) nji · Gji(uh) =

• k

wji

k uk = ˜

Dji ui − uj Hij + gji(uh)

• Dij = min{˜

Dij, ˜ Dji} ≥ C > 0, minimal contribution to edge gradients

• introduce non-linear weights for averaging edge gradients:

nij · G⋄

ij(uh) = Dij

uj − ui Hij + Wij(uh)gij(uh) + Wji(uh)gji(uh)

• non-linear average

Non-linear Diamond Scheme

When averaging the two one-sided edge gradients: nij · G⋄

ij(uh) = Dij

uj − ui Hij + Wij(uh)gij(uh) + Wji(uh)gji(uh)

• non-linear average

the non-linear weights are chosen so that the remainder terms are averaging to zero if their contributions have opposite signs. A possible choice is Wij(uh) = |gji| /(|gij| + |gji|).

Non-linear Diamond Scheme:

Positive weights for reconstructing vertex values:

v T

• main idea: use convex linear combination of barycenters to

express the vertex position;

• linear consistency:
• T wv,T = 1,
• T wv,T(xT − xv) = 0;
• positivity: there exists a constant Cgrid such that

0 < Cgrid ≤ wv,T < 1.

Non-linear Diamond Scheme

• we are solving the non-linear problem

A(w)uh = b(u; w) w = ψ(uh) (edge weights)

• with

A(w) = F + G(uh) + B b(uh; w) = s + r(uh) + g(w) + ˆ b

• main properties:
• there esists at least one solution (uh, w) of the non-linear

discrete problem

• the non-linear discrete operator A(w) is an M-Matrix for

every possible set of admissible edge weights w

Non-linear Diamond Scheme

• fixed-point formulation

uh w

• =
• A(w)−1b(uh; w)

ψ(w)

• ;
• all fixed-point solutions satisfy a Maximum Principle:

G ≤ 0 in Ω ⇒ maxT∈ThuT ≤ max

• 0, maxv∈∂Ω g(xv)
• remark: solving iteratively, only the limit solution satisfies the

maximum principle, but we experienced that intermediate steps are reducing possibly overshoots and undershoots.

• optimal convergence rate from experiments:
• u −

uh

• L2 = O(h2).

Non-linear Diamond Scheme

Approximation of a smooth solution: u(x, y) = x cos (πy), (x, y) ∈ Ω =]0, 1[×]0, 1[.

h Error (L2 norm) ν = 1 (circles), ν = 10−1 (squares), ν = 10−5 (diamonds) β = (−1, −1)t

Non-linear Diamond Scheme

Approximation of a solution with boundary layers: u(x, y) =

• x − e

2(x−1) ν

• y2 − e

3(y−1) ν

• ,

(x, y) ∈ Ω =]0, 1[×]0, 1[.

h Error (L2 norm) ν = 1 (circles), ν = 10−1 (squares), ν = 10−5 (diamonds) ν = 10−5 (stars) on ˆ Ω = [0, 0.8] × [0, 0.8] β = (2, 3)t

Non-linear Diamond Scheme

Approximation of a solution with a boundary and an internal layer: (Burman and Ern, CMAME (2002)) Regular mesh (40 × 40 × 2 triangles) Scheme Overshoot Undershoot Normalized Cost Diamond 2.42E-4 1.03E-6 1 NL Diamond 5.65E-12 3.01E-12 1.12 Unstructured mesh (4352 triangles) Scheme Overshoot Undershoot Normalized Cost Diamond 8.87E-5 2.34E-6 1 NL Diamond 9.48E-11 6.46E-11 1.13

Final remarks

• first variant: non-linear definition of numerical convective

flux using cell-gradient reconstructions and slope limiters;

• second variant: non-linear definition of numerical diffusive

flux to get a Maximum Principle for the discrete solution;

• no a-priori estimates, but experimental evidences that

these schemes may work optimally;

• good shock-capturing behavior.