SLIDE 1
Development of high-order well-balanced schemes for geophysical flows
Development of high-order well-balanced schemes for geophysical flows
Development of high-order well-balanced schemes for geophysical - - PowerPoint PPT Presentation
Development of high-order well-balanced schemes for geophysical - - PowerPoint PPT Presentation
Development of high-order well-balanced schemes for geophysical flows Development of high-order well-balanced schemes for geophysical flows Victor Michel-Dansac Thursday, September 29th, 2016 PhD advisors: Christophe Berthon and Franoise
SLIDE 2
SLIDE 3
Development of high-order well-balanced schemes for geophysical flows Introduction and motivations
1 Introduction and motivations 2 A well-balanced scheme 3 1D numerical experiments 4 Two-dimensional and high-order extensions 5 2D numerical experiments 6 Conclusion and perspectives
SLIDE 4
Development of high-order well-balanced schemes for geophysical flows Introduction and motivations Geophysical flows
Several kinds of destructive geophysical flows
Dam failure (Malpasset, France, 1959) Tsunami (T¯
- hoku, Japan, 2011)
Flood (La Faute sur Mer, France, 2010) Mudslide (Madeira, Portugal, 2010)
1 / 46
SLIDE 5
Development of high-order well-balanced schemes for geophysical flows Introduction and motivations The shallow-water equations
The shallow-water equations and their source terms
∂th + ∂x(hu) = 0 ∂t(hu) + ∂x
- hu2 + 1
2gh2
- = −gh∂xZ − kq|q|
hη (with q = hu) We can rewrite the equations as ∂tW + ∂xF(W) = S(W), with W = h q
- .
x h(x, t)
water surface channel bottom
u(x, t) Z(x)
η = 7/3 and g is the gravitational constant k ≥ 0 is the so-called Manning coefficient: a higher k leads to a stronger Manning friction
2 / 46
SLIDE 6
Development of high-order well-balanced schemes for geophysical flows Introduction and motivations Steady state solutions
Steady state solutions
Definition: Steady state solutions W is a steady state solution iff ∂tW = 0, i.e. ∂xF(W) = S(W). Taking ∂tW = 0 in the shallow-water equations leads to ∂xq = 0 ∂x q2 h + 1 2gh2
- = −gh∂xZ − kq|q|
hη . The steady state solutions are therefore given by q = cst = q0 ∂x q2 h + 1 2gh2
- = −gh∂xZ − kq0|q0|
hη .
3 / 46
SLIDE 7
Development of high-order well-balanced schemes for geophysical flows Introduction and motivations Steady state solutions
Steady states for the friction source term
Assume a flat bottom (Z = cst): the steady states are given by ∂x q2 h + 1 2gh2
- = −kq0|q0|
hη . Assuming smooth steady state solutions and integrating this relation between some x0 ∈ R and x ∈ R yields (with h = h(x) and h0 = h(x0)): − q2 η − 1
- hη−1 − hη−1
- +
g η + 2
- hη+2 − hη+2
- + kq0|q0|(x − x0) = 0.
next step: study of the above nonlinear equation, denoted by χ(h; x, h0, x0, q0) = 0
4 / 46
SLIDE 8
Development of high-order well-balanced schemes for geophysical flows Introduction and motivations Steady state solutions
Steady states for the friction source term
1 We show that ∂χ
∂h(h; x, h0, x0, q0) < 0 if and only if h < hc, where hc = q2 g 1
- 3
. As a consequence, χ(h) is: decreasing for h < hc; increasing for h > hc.
2 In the context of a steady state solution, the Froude number is
defined, using the sound speed c = √gh, by: Fr(h) = u c = q2
- gh3 .
Therefore, Fr(hc) = 1 and the steady state is: supercritical if h < hc; subcritical if h > hc.
5 / 46
SLIDE 9
Development of high-order well-balanced schemes for geophysical flows Introduction and motivations Steady state solutions
Steady states for the friction source term, assuming q0 < 0
h < h c supercritical solution h > hc subcritical solution hc
Sketches of χ(h; x) for h ∈ [0, 0.41], for different values of x, and for hc = 0.25. We are interested in the solutions of χ(h) = 0.
6 / 46
SLIDE 10
Development of high-order well-balanced schemes for geophysical flows Introduction and motivations Steady state solutions
Steady states for the friction source term, assuming q0 < 0
h(x)
blue curve: subcritical solution; red curve: supercritical
7 / 46
SLIDE 11
Development of high-order well-balanced schemes for geophysical flows Introduction and motivations Objectives
Objectives
1 Derive a scheme that:
is well-balanced for the shallow-water equations with friction and/or topography, i.e.:
preservation of all steady states with k = 0 and Z = cst, preservation of all steady states with k = 0 and Z = cst, preservation of steady states with k = 0 and Z = cst;
preserves the non-negativity of the water height; is able to deal with wet/dry transitions, where the friction source term is stiff.
2 Provide two-dimensional and high-order extensions of this
scheme, while keeping the above properties.
8 / 46
SLIDE 12
Development of high-order well-balanced schemes for geophysical flows A well-balanced scheme
1 Introduction and motivations 2 A well-balanced scheme 3 1D numerical experiments 4 Two-dimensional and high-order extensions 5 2D numerical experiments 6 Conclusion and perspectives
SLIDE 13
Development of high-order well-balanced schemes for geophysical flows A well-balanced scheme Non-exhaustive state of the art
Non-exhaustive state of the art
Well-balanced schemes for the shallow-water equations introduction of the well-balance property: Bermudez- Vazquez (1994), Greenberg-LeRoux (1996) preservation of the lake at rest: Audusse et al. (2004), Berthon-Foucher (2012), Audusse et al. (2015) 1D fully well-balanced schemes: Gosse (2000), Castro et al. (2007), Fjordholm et al. (2011), Xing et al. (2011), Berthon-Chalons (2016) 1D high-order schemes that preserve steady states: Castro et al. (2006), Castro Díaz et al. (2013) 2D schemes preserving the lake at rest on unstructured meshes: Duran et al. (2013), Clain-Figueiredo (2014) for the friction source term: Liang-Marche (2009), Chertock et al. (2015)
9 / 46
SLIDE 14
Development of high-order well-balanced schemes for geophysical flows A well-balanced scheme Structure of the scheme
The HLL scheme
To approximate solutions of ∂tW + ∂xF(W) = 0, the HLL scheme (Harten, Lax, van Leer (1983)) may be chosen; it uses the approximate Riemann solver W, displayed on the right.
W HLL WL WR λL x t λR −∆x/2 ∆x/2
The consistency condition (as per Harten and Lax) holds if: 1 ∆x ∆x/2
−∆x/2
- W(∆t, x; WL, WR)dx =
1 ∆x ∆x/2
−∆x/2
WR(∆t, x; WL, WR)dx, which gives WHLL = λRWR − λLWL λR − λL − F(WR) − F(WL) λR − λL = hHLL qHLL
- .
Note that, if hL > 0 and hR > 0, then hHLL > 0 for |λL| and |λR| large enough.
10 / 46
SLIDE 15
Development of high-order well-balanced schemes for geophysical flows A well-balanced scheme Structure of the scheme
Modification of the HLL scheme
With Y (t, x) = x, we rewrite the shallow-water equations with a generic source term S as follows: ∂th + ∂xq = 0, ∂tq + ∂x q2 h + 1 2gh2
- − S∂xY = 0,
∂tY = 0. The equation ∂tY = 0 induces a stationary wave associated to the source term; we also note that q is a Riemann invariant for this wave. To approximate solutions of ∂tW + ∂xF(W) = S(W), we thus use the approximate Riemann solver displayed on the right (assuming λL < 0 < λR).
WL WR λL λR W ∗
L
W ∗
R
11 / 46
SLIDE 16
Development of high-order well-balanced schemes for geophysical flows A well-balanced scheme Structure of the scheme
Modification of the HLL scheme
We now wish to apply the Harten-Lax consistency condition to ∂tW + ∂xF(W) = S(W). Recall this condition: 1 ∆x ∆x/2
−∆x/2
- W(∆t, x; WL, WR)dx =
1 ∆x ∆x/2
−∆x/2
WR(∆t, x; WL, WR)dx, first step: compute 1 ∆x ∆x/2
−∆x/2
- W(∆t, x)dx (straightforward)
1 ∆x ∆x/2
−∆x/2
- W(∆t, x)dx = WL+WR
2 − λR ∆t ∆x(WR−W ∗
R) + λL
∆t ∆x(WL−W ∗
L)
second step: compute 1 ∆x ∆x/2
−∆x/2
WR(∆t, x)dx
12 / 46
SLIDE 17
Development of high-order well-balanced schemes for geophysical flows A well-balanced scheme Structure of the scheme
Modification of the HLL scheme
1 ∆t 1 ∆x ∆t ∆x/2
−∆x/2
(∂tWR + ∂xF(WR) ) dx dt = 0
13 / 46
SLIDE 18
Development of high-order well-balanced schemes for geophysical flows A well-balanced scheme Structure of the scheme
Modification of the HLL scheme
1 ∆t 1 ∆x ∆t ∆x/2
−∆x/2
(∂tWR + ∂xF(WR) ) dx dt = 0 0 = 1 ∆t 1 ∆x ∆x/2
−∆x/2
WR(∆t, x)dx − ∆x/2
−∆x/2
WR(0, x)dx
- +
1 ∆t 1 ∆x ∆t F(WR)
- t, −∆x
2
- dt −
∆t F(WR)
- t, ∆x
2
- dt
- WL
WR λL λR W ∗
L
W ∗
R
−∆x/2 ∆x/2 t x ∆t
13 / 46
SLIDE 19
Development of high-order well-balanced schemes for geophysical flows A well-balanced scheme Structure of the scheme
Modification of the HLL scheme
1 ∆t 1 ∆x ∆t ∆x/2
−∆x/2
(∂tWR + ∂xF(WR) ) dx dt = 0 0 = 1 ∆t 1 ∆x ∆x/2
−∆x/2
WR(∆t, x)dx − ∆x/2
−∆x/2
WR(0, x)dx
- +
1 ∆t 1 ∆x ∆t F(WR)
- t, −∆x
2
- dt −
∆t F(WR)
- t, ∆x
2
- dt
- 1
∆x ∆x/2
−∆x/2
WR(∆t, x)dx = WL + WR 2 − ∆t ∆x(F(WR) − F(WL))
13 / 46
SLIDE 20
Development of high-order well-balanced schemes for geophysical flows A well-balanced scheme Structure of the scheme
Modification of the HLL scheme
1 ∆t 1 ∆x ∆t ∆x/2
−∆x/2
(∂tWR + ∂xF(WR) − S(WR)) dx dt = 0
13 / 46
SLIDE 21
Development of high-order well-balanced schemes for geophysical flows A well-balanced scheme Structure of the scheme
Modification of the HLL scheme
1 ∆t 1 ∆x ∆t ∆x/2
−∆x/2
(∂tWR + ∂xF(WR) − S(WR)) dx dt = 0 0 = 1 ∆t 1 ∆x ∆x/2
−∆x/2
WR(∆t, x)dx − ∆x/2
−∆x/2
WR(0, x)dx
- +
1 ∆t 1 ∆x ∆t F(WR)
- t, −∆x
2
- dt −
∆t F(WR)
- t, ∆x
2
- dt
- −
1 ∆t 1 ∆x ∆t ∆x/2
−∆x/2
S(WR)(t, x) dx dt
WL WR λL λR W ∗
L
W ∗
R
−∆x/2 ∆x/2 t x ∆t
13 / 46
SLIDE 22
Development of high-order well-balanced schemes for geophysical flows A well-balanced scheme Structure of the scheme
Modification of the HLL scheme
1 ∆t 1 ∆x ∆t ∆x/2
−∆x/2
(∂tWR + ∂xF(WR) − S(WR)) dx dt = 0 0 = 1 ∆t 1 ∆x ∆x/2
−∆x/2
WR(∆t, x)dx − ∆x/2
−∆x/2
WR(0, x)dx
- +
1 ∆t 1 ∆x ∆t F(WR)
- t, −∆x
2
- dt −
∆t F(WR)
- t, ∆x
2
- dt
- −
1 ∆t 1 ∆x ∆t ∆x/2
−∆x/2
S(WR)(t, x) dx dt 1 ∆x ∆x/2
−∆x/2
WR(∆t, x)dx ≃ WL + WR 2 − ∆t ∆x(F(WR) − F(WL)) + S∆t
13 / 46
SLIDE 23
Development of high-order well-balanced schemes for geophysical flows A well-balanced scheme Structure of the scheme
Modification of the HLL scheme
We have 4 unknowns to determine: W ∗
L =
h∗
L
q∗
L
- and W ∗
R =
h∗
R
q∗
R
- .
q is a 0-Riemann invariant we take q∗
L = q∗ R = q∗ (relation 1)
Harten-Lax consistency gives us the following two relations: λRh∗
R − λLh∗ L = (λR − λL)hHLL (relation 2)
q∗ = qHLL + S∆x λR − λL (relation 3) next step: obtain a fourth relation
14 / 46
SLIDE 24
Development of high-order well-balanced schemes for geophysical flows A well-balanced scheme The full scheme for a general source term
Obtaining an additional relation
Assume that WL and WR define a steady state, i.e. that they satisfy the following discrete version of the steady relation ∂xF(W) = S(W) (where [X] = XR − XL): 1 ∆x
- q2
1 h
- + g
2
- h2
= S. For the steady state to be preserved, it is sufficient to have h∗
L = hL, h∗ R = hR
and q∗ = q0.
WL WR
Assuming a steady state, we show that q∗ = q0, as follows: q∗ = qHLL + S∆x λR − λL = q0 − 1 λR − λL
- q2
1 h
- + g
2
- h2
− S∆x
- = q0.
15 / 46
SLIDE 25
Development of high-order well-balanced schemes for geophysical flows A well-balanced scheme The full scheme for a general source term
Obtaining an additional relation
In order to determine an addition relation, we consider the discrete steady relation, satisfied when WL and WR define a steady state: q2 1 hR − 1 hL
- + g
2
- (hR)2 − (hL)2
= S∆x. To ensure that h∗
L = hL and h∗ R = hR, we impose that h∗ L and h∗ R
satisfy the above relation, as follows: q2 1 h∗
R
− 1 h∗
L
- + g
2
- (h∗
R)2 − (h∗ L)2
= S∆x.
16 / 46
SLIDE 26
Development of high-order well-balanced schemes for geophysical flows A well-balanced scheme The full scheme for a general source term
Determination of h∗
L and h∗ R
The intermediate water heights satisfy the following relation: −q2 h∗
R − h∗ L
h∗
Lh∗ R
- + g
2(h∗
L + h∗ R)(h∗ R − h∗ L) = S∆x.
Recall that q∗ is known and is equal to q0 for a steady state. Instead of the above relation, we choose the following linearization: −(q∗)2 hLhR (h∗
R − h∗ L) + g
2(hL + hR)(h∗
R − h∗ L) = S∆x,
which can be rewritten as follows: −(q∗)2 hLhR + g 2(hL + hR)
- α
(h∗
R − h∗ L) = S∆x.
17 / 46
SLIDE 27
Development of high-order well-balanced schemes for geophysical flows A well-balanced scheme The full scheme for a general source term
Determination of h∗
L and h∗ R
With the consistency relation between h∗
L and h∗ R, the intermediate
water heights satisfy the following linear system:
- α(h∗
R − h∗ L) = S∆x,
λRh∗
R − λLh∗ L = (λR − λL)hHLL.
Using both relations linking h∗
L and h∗ R, we obtain
h∗
L = hHLL −
λRS∆x α(λR − λL), h∗
R = hHLL −
λLS∆x α(λR − λL), where α = −(q∗)2 hLhR + g 2(hL + hR)
- with q∗ = qHLL +
S∆x λR − λL .
18 / 46
SLIDE 28
Development of high-order well-balanced schemes for geophysical flows A well-balanced scheme The full scheme for a general source term
Correction to ensure non-negative h∗
L and h∗ R
However, these expressions of h∗
L and h∗ R do not guarantee that the
intermediate heights are non-negative: instead, we use the following cutoff (see Audusse, Chalons, Ung (2014)): h∗
L = min
- hHLL −
λRS∆x α(λR − λL)
- +
,
- 1 − λR
λL
- hHLL
- ,
h∗
R = min
- hHLL −
λLS∆x α(λR − λL)
- +
,
- 1 − λL
λR
- hHLL
- .
Note that this cutoff does not interfere with: the consistency condition λRh∗
R − λLh∗ L = (λR − λL)hHLL;
the well-balance property, since it is not activated when WL and WR define a steady state.
19 / 46
SLIDE 29
Development of high-order well-balanced schemes for geophysical flows A well-balanced scheme The full scheme for a general source term
Summary
The two-state approximate Riemann solver with intermediate states W ∗
L =
h∗
L
q∗
- and W ∗
R =
h∗
R
q∗
- given by
q∗ = qHLL + S∆x λR − λL , h∗
L = min
- hHLL −
λRS∆x α(λR − λL)
- +
,
- 1 − λR
λL
- hHLL
- ,
h∗
R = min
- hHLL −
λLS∆x α(λR − λL)
- +
,
- 1 − λL
λR
- hHLL
- ,
is consistent, non-negativity-preserving and well-balanced. next step: determination of S according to the source term definition (topography and/or friction).
20 / 46
SLIDE 30
Development of high-order well-balanced schemes for geophysical flows A well-balanced scheme The cases of the topography and friction source terms
The topography source term
We now consider S(W) = St(W) = −gh∂xZ: the smooth steady states are governed by ∂x q2 h
- + g
2∂x
- h2
= −gh∂xZ, q2 2 ∂x 1 h2
- + g∂x(h + Z) = 0,
− − − − − − − →
discretization
q2 1 h
- + g
2
- h2
= St∆x, q2 2 1 h2
- + g[h + Z] = 0.
We can exhibit an expression of q2
0 and thus obtain
St = −g 2hLhR hL + hR [Z] ∆x + g 2∆x [h]3 hL + hR . However, when ZL = ZR, we have St = O(∆x), i.e. a loss of consistency with St (see for instance Berthon, Chalons (2016)).
21 / 46
SLIDE 31
Development of high-order well-balanced schemes for geophysical flows A well-balanced scheme The cases of the topography and friction source terms
The topography source term
Instead, we set, for some constant C > 0, St = −g 2hLhR hL + hR [Z] ∆x + g 2∆x [h]3
c
hL + hR , [h]c =
- hR − hL
if |hR − hL| ≤ C∆x, sgn(hR − hL) C∆x
- therwise.
Theorem: Well-balance for the topography source term If WL and WR define a smooth steady state, i.e. if they satisfy q2 2 1 h2
- + g[h + Z] = 0,
then we have W ∗
L = WL and W ∗ R = WR and the approximate
Riemann solver is well-balanced.
22 / 46
SLIDE 32
Development of high-order well-balanced schemes for geophysical flows A well-balanced scheme The cases of the topography and friction source terms
The friction source term
We consider, in this case, S(W) = Sf(W) = −kq|q|h−η. The average of Sf we choose is Sf = −k¯ q|¯ q|h−η, with ¯ q the harmonic mean of qL and qR (note that ¯ q = q0 at the equilibrium); h−η a well-chosen discretization of h−η, depending on hL and hR, and ensuring the well-balance property. We determine h−η using the same technique (with µ0 = sgn(q0)):
∂x q2 h
- + g
2∂x
- h2
= −kq0|q0|h−η, q2 ∂xhη−1 η − 1 − g ∂xhη+2 η + 2 = kq0|q0|, − − − − − − − →
discretization
q2 1 h
- + g
2
- h2
= −kµ0q2
0h−η∆x,
q2
- hη−1
η − 1 − g
- hη+2
η + 2 = kµ0q2
0∆x.
23 / 46
SLIDE 33
Development of high-order well-balanced schemes for geophysical flows A well-balanced scheme The cases of the topography and friction source terms
The friction source term
The expression for q2
0 we obtained is now used to get:
h−η = [h2] 2 η + 2 [hη+2] − µ0 k∆x 1 h
- + [h2]
2 [hη−1] η − 1 η + 2 [hη+2]
- ,
which gives Sf = −k¯ q|¯ q|h−η (h−η is consistent with h−η). Theorem: Well-balance for the friction source term If WL and WR define a smooth steady state, i.e. verify q2
- hη−1
η − 1 + g
- hη+2
η + 2 = −kq0|q0|∆x, then we have W ∗
L = WL and W ∗ R = WR and the approximate
Riemann solver is well-balanced.
24 / 46
SLIDE 34
Development of high-order well-balanced schemes for geophysical flows A well-balanced scheme The cases of the topography and friction source terms
Friction and topography source terms
With both source terms, the scheme preserves the following discretization of the steady relation ∂xF(W) = S(W): q2 1 h
- + g
2
- h2
= St∆x + Sf∆x. The intermediate states are therefore given by: q∗ = qHLL + (St + Sf)∆x λR − λL ; h∗
L = min
- hHLL − λR(St + Sf)∆x
α(λR − λL)
- +
,
- 1 − λR
λL
- hHLL
- ;
h∗
R = min
- hHLL − λL(St + Sf)∆x
α(λR − λL)
- +
,
- 1 − λL
λR
- hHLL
- .
25 / 46
SLIDE 35
Development of high-order well-balanced schemes for geophysical flows A well-balanced scheme The cases of the topography and friction source terms
The full Godunov-type scheme
x t tn+1 tn xi xi− 1
2
xi+ 1
2
W n
i
W R
i− 1
2
W L
i+ 1
2
λR
i− 1
2
λL
i+ 1
2
- W ∆(x, tn+1)
We define W n+1
i
= 1 ∆x xi+ 1
2
xi− 1
2
W ∆(x, tn+1)dx: then W n+1
i
= W n
i − ∆t
∆x
- λL
i+ 1
2
- W L
i+ 1
2 − W n
i
- − λR
i− 1
2
- W R
i− 1
2 − W n
i
- ,
which can be rewritten, after straightforward computations,
W n+1
i
= W n
i − ∆t
∆x
- F n
i+ 1
2 − F n
i− 1
2
- + ∆t
(St)n
i− 1
2+(St)n
i+ 1
2
2 + (Sf)n
i− 1
2+(Sf)n
i+ 1
2
2 .
26 / 46
SLIDE 36
Development of high-order well-balanced schemes for geophysical flows A well-balanced scheme The cases of the topography and friction source terms
Summary
We have presented a scheme that: is consistent with the shallow-water equations with friction and topography; is well-balanced for friction and topography steady states; preserves the non-negativity of the water height; is not able to correctly approximate wet/dry interfaces due to the stiffness of the friction: we require a semi-implicitation of the friction source term. next step: introduction of this semi-implicitation
27 / 46
SLIDE 37
Development of high-order well-balanced schemes for geophysical flows A well-balanced scheme Source terms contribution to the finite volume scheme
Semi-implicit finite volume scheme
We use a splitting method with an explicit treatment of the flux and the topography and an implicit treatment of the friction.
1 explicitly solve ∂tW + ∂xF(W) = St(W) to get
W
n+ 1
2
i
= W n
i − ∆t
∆x
- Fn
i+ 1
2 − Fn
i− 1
2
- + ∆t
- 1
2
- (St)n
i− 1
2 + (St)n
i+ 1
2
- 2 implicitly solve ∂tW = Sf(W) to get
hn+1
i
= h
n+ 1
2
i
IVP: ∂tq = −kq|q|(hn+1
i
)−η q(xi, tn) = q
n+ 1
2
i
qn+1
i
28 / 46
SLIDE 38
Development of high-order well-balanced schemes for geophysical flows A well-balanced scheme Source terms contribution to the finite volume scheme
Semi-implicit finite volume scheme
Solving the IVP yields: qn+1
i
= (hn+1
i
)ηq
n+ 1
2
i
(hn+1
i
)η + k ∆t
- q
n+ 1
2
i
- .
We use the following approximation of (hn+1
i
)η, which provides us with an expression of qn+1
i
that is equal to q0 at the equilibrium: (hη)n+1
i
= 2µ
n+ 1
2
i
µn
i
- h−ηn+1
i− 1
2
+
- h−ηn+1
i+ 1
2
+ k ∆t µ
n+ 1
2
i
qn
i .
semi-implicit treatment of the friction source term scheme able to model wet/dry transitions scheme still well-balanced and non-negativity-preserving
29 / 46
SLIDE 39
Development of high-order well-balanced schemes for geophysical flows 1D numerical experiments