A Numerical Model for Two-Phase Shallow Granular Flows over Variable - - PowerPoint PPT Presentation

A Numerical Model for Two-Phase Shallow Granular Flows

• ver Variable Topography

Marica Pelanti Département de Mathématiques et Applications, ENS Paris, France Joint work with: François Bouchut (DMA-ENS) Anne Mangeney (IPGP)

GdR CHANT Workshop on Models and Numerical Methods for Granular Materials, ENPC, November 19-21, 2007

– p. 1/31

Two-Phase Granular Flow

Outline

• I. Objectives and Motivations
• II. Physical and Mathematical Model

⋆ Eigenvalues analysis

• III. Numerical Solution Method

⋆ Solver for the homogeneous system ⋆ Topography source terms, inter-phase drag source terms ⋆ Numerical experiments

• IV. Summary and Future Work

– p. 2/31

Objectives and Motivations

Goal: Development of a numerical two-phase shallow flow model for mixtures of solid grains and fluid over variable basal surface. Ultimate objective: Numerical simulation of geophysical flows such as avalanches and debris flows over natural terrain. (Project DMA-ENS & IPGP.)

• Gravitational geophysical flows typically involve both solid granular material

and interstitial fluid.

• Inter-phase forces influence flow mechanics: flow deformation, mobility,

run-out, deposit.

Source: USGS

• I. Objectives and Motivations

– p. 3/31

Thin Layer Granular Flow Models – State of the Art in brief Thin layer (shallow flow) models: H/L ≪ 1. H, L = characteristic flow depth and length. Continuum flow equations are scaled and depth-averaged.

• First 1D single-phase dry granular flow model: Savage and Hutter, 1989.

Extensive work on dry granular flows: 2D models, complex topography.

• I. Objectives and Motivations

– p. 4/31

Thin Layer Granular Flow Models – State of the Art in brief Thin layer (shallow flow) models: H/L ≪ 1. H, L = characteristic flow depth and length. Continuum flow equations are scaled and depth-averaged.

• First 1D single-phase dry granular flow model: Savage and Hutter, 1989.

Extensive work on dry granular flows: 2D models, complex topography.

• Solid-Fluid Mixture Model (Iverson, 1997; Iverson and Denlinger, 2001)

(Similar mixture theory approach: Pudasaini–Wang–Hutter, 2005).

Hp: 1) constant fluid volume fraction; 2) fluid velocity = solid velocity. System: mass and momentum equations for the mixture. No inherent pore fluid motion description ⇒ needs supplementary specification of pore fluid pressure evolution ⋆ 2D model, treats irregular topography; Numerical simulation of many laboratory experiments and debris-flow-flume tests.

• I. Objectives and Motivations

– p. 4/31

Thin Layer Granular Flow Models – State of the Art in brief Thin layer (shallow flow) models: H/L ≪ 1. H, L = characteristic flow depth and length. Continuum flow equations are scaled and depth-averaged.

• First 1D single-phase dry granular flow model: Savage and Hutter, 1989.

Extensive work on dry granular flows: 2D models, complex topography.

• Solid-Fluid Mixture Model (Iverson, 1997; Iverson and Denlinger, 2001)

(Similar mixture theory approach: Pudasaini–Wang–Hutter, 2005).

Hp: 1) constant fluid volume fraction; 2) fluid velocity = solid velocity. System: mass and momentum equations for the mixture. No inherent pore fluid motion description ⇒ needs supplementary specification of pore fluid pressure evolution ⋆ 2D model, treats irregular topography; Numerical simulation of many laboratory experiments and debris-flow-flume tests.

• A Two-Phase Model: Pitman and Le, 2005.

⋆ Retains mass and momentum equations for both solid and fluid phases. However: non-conservative mixture momentum eq.; no general topography. Numerical method only for reduced model that ignores fluid inertial terms.

• I. Objectives and Motivations

– p. 4/31

Two-Phase Shallow Granular Flow

Physical and Mathematical Model

Pitman–Le approach (2005) Consider thin layer of a mixture of solid grains and fluid flowing over a smooth basal surface. Two-phase flow equations [Anderson and Jackson, 1967] ∂t(ρsϕ) + ∇ · (ρsϕVs) = 0 , ρsϕ(∂tVs + (Vs · ∇)Vs) = ∇ · Ts + ϕ∇ · Tf + I + ρsϕg , ∂t(ρf(1 − ϕ)) + ∇ · (ρf(1 − ϕ)Vf) = 0 , ρf(1 − ϕ)(∂tVf + (Vf · ∇)Vf) = (1 − ϕ)∇ · Tf − I + ρf(1 − ϕ)g . ρs, ρf = solid and fluid specific densities; ϕ = solid volume fraction; Vs, Vf = velocities; Ts, Tf = stress tensors; g = gravity vector; I = non-buoyancy interaction forces (e.g. drag).

• II. Physical and Mathematical Model

– p. 5/31

Two-Phase Shallow Granular Flow Model Material Constitutive Behaviour

• Solid and fluid are incompressible: material densities ρs, ρf = constant.
• Fluid is inviscid; only fluid stress is a pressure.
• Solid modeled as Coulomb material (as Savage–Hutter).

Coulomb friction law for solid shear stresses: T xz

s

= −sgn(Vs)νT zz

s ,

ν ≥ 0 , Vs = solid sliding velocity. Earth-pressure relation for lateral normal stresses: T xx

s

= KT zz

s . (K = 1)

Boundary Conditions

• Free upper surface stress-free, and material surface for both phases.
• Both solid and fluid motion tangent to the basal surface (no-slip condition).

Moreover: Only non-buoyancy interaction force I is drag: I = D(Vf − Vs). Under the shallow flow hypothesis H/L ≪ 1 : Scale and depth-average the governing two-phase flow equations.

• II. Physical and Mathematical Model

– p. 6/31

Model Equations (1D and hp. small topography slopes)

Depth-Averaged Solid and Fluid Mass and Momentum Equations ∂ ∂t (ϕh) + ∂ ∂x (ϕhvs) = 0 , ∂ ∂t (ϕhvs) + ∂ ∂x

• ϕhv2

s + g 1−γ

2 ϕh2

• + γϕg

2 ∂h2 ∂x = −gϕh ∂b ∂x − sgn(vs)νb g(1 − γ)ϕh + γDh(vf − vs) , ∂ ∂t ((1 − ϕ)h) + ∂ ∂x ((1 − ϕ)hvf) = 0 , ∂ ∂t ((1 − ϕ)hvf) + ∂ ∂x

• (1 − ϕ)hv2

f

• + (1 − ϕ)g

2 ∂h2 ∂x = −g(1 − ϕ)h ∂b ∂x − Dh(vf − vs) .

h = flow depth; ϕ = depth-averaged solid volume fraction; vs, vf = averaged solid and fluid velocities; γ =

ρf ρs < 1,

ρs, ρf = material specific densities (constant); b(x) = bottom topography; g = gravity constant; νb = tan δbed, δbed = basal friction angle; D = average drag function (= ¯ D/ρf).

• II. Physical and Mathematical Model

– p. 7/31

Presented model vs. original Pitman–Le model Presented model: variant of the original Pitman–Le model. Different averaging approximation of fluid motion equation ⇒ Different fluid momentum balance: Pitman–Le:

∂ ∂t (hvf) + ∂ ∂x

• hv2

f

• + g

2 ∂h2 ∂x = 0.

Here:

∂ ∂t ((1 − ϕ)hvf) + ∂ ∂x

• (1 − ϕ)hv2

f

• + (1 − ϕ) g

2 ∂h2 ∂x = 0.

Difference: τ ≡ hvf (∂t(1 − ϕ) + vf∂x(1 − ϕ)) = (1 − ϕ)vf∂x(ϕh(vs − vf)). ⇒ Different mixture momentum balance. Here: conservative equation for the momentum of the mixture

∂ ∂t ` (ϕvs + γ(1 − ϕ)vf)h ´ + ∂ ∂x “` ϕv2

s + γ(1 − ϕ)v2 f

´ h + g 2 (ϕ + γ(1 − ϕ))h2” = 0 .

Consistent with conservative mixture momentum equation of two-phase flow system before averaging and expected physical behaviour.

• II. Physical and Mathematical Model

– p. 8/31

Two-Phase Shallow Granular Flow Equations. Formulation (hs, hsvs, hf, hfvf). Set hs = ϕ h, hf = (1 − ϕ) h, ϕ = solid volume fraction.

(Here no friction.)

∂hs ∂t + ∂ ∂x (hsvs) = 0 , ∂ ∂t (hsvs) + ∂ ∂x

• hsv2

s + g

2h2

s + g 1−γ

2 hshf

• + γghs

∂hf ∂x = −ghs ∂b ∂x + γF D, ∂hf ∂t + ∂ ∂x (hf vf) = 0 , ∂ ∂t (hfvf) + ∂ ∂x

• hfv2

f + g

2h2

f

• + g hf

∂hs ∂x = −ghf ∂b ∂x − F D.

Drag force F D = D(hs + hf)(vf − vs); γ = ρf/ρs .

Note: Similar to two-layer shallow flow model, except additional cross term

∂ ∂x

• g 1−γ

2 hshf

• in the solid momentum balance.
• II. Physical and Mathematical Model

– p. 9/31

Eigenvalues Analysis

Consider ∂tq + A(q)∂xq = 0 , q ∈ R4 , A ∈ R4×4 . In general: eigenvalues λk , k = 1, . . . , 4, cannot be expressed explicitly. Define: a = √gh

and

β =

• 1

2(1 − ϕ)(1 − γ) < 1 .

If vf = vs ≡ v then A has real distinct eigenvalues (ϕ = 1) λ1,4 = v ∓ a

and

λ2,3 = v ∓ aβ . We can show that: There are always two real external eigenvalues λ1,4 , and, moreover min(vf, vs) − a ≤ λ1 < Re(λ2) ≤ Re(λ3) < λ4 ≤ max(vf, vs) + a . Furthermore:

• If |vs − vf| < 2aβ or |vs − vf| > 2a then all the eigenvalues are real and

distinct (ϕ = 1) ⇒ the system is strictly hyperbolic.

• If 2aβ < |vs − vf| < 2a then the internal eigenvalues may be complex.

Hyperbolicity at least when the velocity difference |vs − vf| is sufficiently small.

• II. Physical and Mathematical Model

– p. 10/31

Eigenvectors

Right Eigenvectors q = (hs, hsvs, hf, hfvf)T. Assume hs, hf = 0. For k = 1, . . . , 4 : rk =    1 λk ξk ξkλk   , ξk = (λk − vs)2 − g

• hs + 1−γ

2 hf

• g 1+γ

2 hs

= ghf (λk − vf)2 − ghf . Note: Can show that 1st and 4th fields are genuinely nonlinear: ∇λk · rk = 0, ∀q. Left Eigenvectors, L = R−1 lk = nk P ′(λk) , nk = ( ϑs,k(λk − 2vs), ϑs,k, ϑf(λk − 2vf), ϑf) , P(λ) = characteristic polynomial, ϑs,k = (λk − vf)2 − ghf = g hf ξk

and

ϑf = g 1 + γ 2 hs .

• II. Physical and Mathematical Model

– p. 11/31

Numerical Solution

• Assume |vs − vf| small enough so that the system is strictly hyperbolic.

Class of methods used: Godunov-type Finite Volume Schemes (Schemes based on Riemann Solvers) Difficulties:

• Non-conservative system. Many well-established efficient finite volume

schemes: only for conservation laws. (New difficulty with respect to dry granular flow models and mixture models.)

• Topography source terms need to be discretized so that the method is

well-balanced = it preserves steady states and captures accurately perturbations. Well-known difficulty for systems with sources.

• Positivity preservation: computed values of flow depth and phase volume

fractions must be positive. Important to handle interfaces between flow fronts and dry bed zones (h = 0). → still to be addressed. Here we will consider regimes without dry bed areas.

• III. Numerical Solution

– p. 12/31

Godunov-Type Schemes

Riemann problem: ∂tq + A(q)∂xq = 0 with I.C. q(x, 0) =

• qℓ

if x ≤ ¯

x , qr

if x > ¯

x .

tn tn+1 Qn

i

Qn+1

i

Qn

i → approximate solution on cell (xi−1/2, xi+1/2).

Discontinuities at cell interfaces ⇒ Riemann problems.

1) At each cell interface xi+1/2 between Qn

i and Qn i+1 → solve Riemann

problem with data Qn

i and Qn i+1.

2) Use solution of local Riemann problems to update solution Qn

i → Qn+1 i

. Riemann Solver: → Set of waves Wk and speeds sk representing the (approximate) Riemann solution structure.

• ∆q ≡ Qi+1−Qi =

k Wk

• For conservative systems ∂tq + ∂xF(q) = 0: F(Qi+1)−F(Qi) =

k skWk

Wave-Propagation Algorithm: Qn+1

i

= Qn

i − ∆t ∆x(A+∆Qi−1/2 + A−∆Qi+1/2),

fluctuations: A±∆Qi+1/2 =

k(sk i+1/2)±Wk i+1/2 ,

s+=max(s, 0), s−=min(s, 0).

• III. Numerical Solution

– p. 13/31

Numerical Solution Homogeneous System

(b(x) = 0, D = 0)

∂tq + ∂xf(q) + w(q, ∂xq) = 0 , q = (hs, hsvs, hf, hfvf)T , f(q) =

• hsvs, hs v2

s + g 2 h2 s + g 1−γ 2 hs hf, hf vf, hf v2 f + g 2 h2 f

T , w(q, ∂xq) = (0, γghs ∂xhf, 0, ghf ∂xhs)

T .

⋆ Solid and fluid mass equations are conservative. ⋆ Mixture momentum equation is conservative: ∂tm + ∂xfm(q) = 0 , m = hsvs + γ hfvf , fm(q) = f (2)(q) + γ f (4)(q) + γ g hshf . ⋆ Non-conservative products γghs

∂hf ∂x , ghf ∂hs ∂x in the momentum balances

couple sets of equations of the two phases ⇒ avoid uncoupled schemes that may generate instabilities. We employ a Roe-type Riemann Solver.

• III. Numerical Solution

– p. 14/31

Roe-type Riemann Solver

Consider the quasi-linear form of the system ∂tq + A(q)∂xq = 0 . At each local cell interface xi+1/2 between solution values Qi and Qi+1 solve a Riemann problem for a linearized system ∂tq + ˆ A(Qi, Qi+1)∂xq = 0 with initial data Qi and Qi+1. The Roe matrix ˆ A(Qi, Qi+1) is defined so as to guarantee conservation for the mass of each phase and for the momentum of the mixture: f (p)(Qi+1) − f (p)(Qi) = ˆ A(p,:) (Qi+1 − Qi) , p = 1, 3 , fm(Qi+1) − fm(Qi) = ( ˆ A(2,:) + γ ˆ A(4,:))(Qi+1 − Qi) . We take ˆ A = A(ˆ hs, ˆ hf, ˆ vs, ˆ vf), with the choice ˆ hθ = hθ,i + hθ,i+1 2

and

ˆ vθ =

• hθ,i vθ,i +
• hθ,i+1 vθ,i+1
• hθ,i +
• hθ,i+1

, θ = s, f . Then: waves Wk = αk ˆ rk , ∆q = 4

k=1 αk ˆ

rk , and speeds sk = ˆ λk , k = 1, . . . , 4. {ˆ rk , ˆ λk} = eigenpairs of ˆ A.

• III. Numerical Solution

– p. 15/31

Wave Propagation Algorithms (LeVeque, 1997) – F-Wave Formulation

Basic software: CLAWPACK.

1) Classical Riemann solver: ∆q =

k Wk ; Roe: Wk =αkˆ

rk , sk = ˆ λk 2) F-wave Approach: For a conservative system ∂tq + ∂xF(q) = 0 , decompose flux jump ∆F ≡ F(Qi+1) − F(Qi) =

k Zk .

Local Riemann solution approximated by f-waves Zk and associated speeds sk. Roe: Zk = ζk ˆ rk , sk = ˆ λk , {ˆ rk , ˆ λk} = eigenpairs of Roe matrix for F′(q). Fluctuations A−∆Qi+1/2 =

• k:sk

i+1/2<0

Zk

i+1/2 ,

A+∆Qi+1/2 =

• k:sk

i+1/2>0

Zk

i+1/2 .

Algorithm Qn+1

i

= Qn

i − ∆t ∆x(A+∆Qi−1/2 + A−∆Qi+1/2)

(2) can be equivalent to (1), but useful framework to include source terms.

• III. Numerical Solution

– p. 16/31

Wave Propagation Algorithms (LeVeque, 1997) – F-Wave Formulation

Basic software: CLAWPACK.

1) Classical Riemann solver: ∆q =

k Wk ; Roe: Wk =αkˆ

rk , sk = ˆ λk 2) F-wave Approach: For a conservative system ∂tq + ∂xF(q) = 0 , decompose flux jump ∆F ≡ F(Qi+1) − F(Qi) =

k Zk .

Local Riemann solution approximated by f-waves Zk and associated speeds sk. Roe: Zk = ζk ˆ rk , sk = ˆ λk , {ˆ rk , ˆ λk} = eigenpairs of Roe matrix for F′(q). Fluctuations A−∆Qi+1/2 =

• k:sk

i+1/2<0

Zk

i+1/2 ,

A+∆Qi+1/2 =

• k:sk

i+1/2>0

Zk

i+1/2 .

Algorithm (high-resolution) Qn+1

i

= Qn

i − ∆t ∆x(A+∆Qi−1/2 + A−∆Qi+1/2)− ∆t ∆x(F c i+1/2 − F c i−1/2)

F c

i+1/2 = correction fluxes for second order accuracy

(2) can be equivalent to (1), but useful framework to include source terms.

• III. Numerical Solution

– p. 16/31

The Roe-Type Solver in the F-Wave Framework Difficulty: Here non-conservative system ∂tq + ∂xf(q) + w(q, ∂xq) = 0 , w(q, ∂xq) = (0, γghs ∂xhf, 0, ghf ∂xhs)

T .

We lack a flux function F to be used for f-wave decomposition.

• III. Numerical Solution

– p. 17/31

The Roe-Type Solver in the F-Wave Framework Difficulty: Here non-conservative system ∂tq + ∂xf(q) + w(q, ∂xq) = 0 , w(q, ∂xq) = (0, γghs ∂xhf, 0, ghf ∂xhs)

T .

We lack a flux function F to be used for f-wave decomposition. Nonetheless can still formulate our Roe-type method into the f-wave framework: Take local linearization of w(q, ∂xq) consistent with Roe linearization and define approximate flux difference ∆ ˜ F = ∆f + (0, γ g ˆ hs ∆hf, 0, g ˆ hf ∆hs)T . Then decompose ∆ ˜ F = 4

k=1 ζk ˆ

rk , and set Zk = ζk ˆ rk , sk = ˆ λk , k = 1, . . . , 4 . {ˆ rk , ˆ λk} = eigenpairs of ˆ A. Note: ∆ ˜ F = ˆ A ∆q (classical Roe property).

• III. Numerical Solution

– p. 17/31

Topography Source Terms

Consider system with topography terms: ∂tq + A(q)∂xq = ψ b(q) , q = (hs, hsvs, hf, hfvf)T , ψ b(q) = − (0, ghs∂xb, 0, ghf ∂xb)

T ,

b = b(x) . Need well-balancing: efficient modeling of equilibrium and quasi-equilibrium states → A(q)∂xq ≈ ψ b(q) . Steady states at rest (vs = vf = 0): hs + hf + b = const.

and

hf hs = const. That is h + b = const.

and

ϕ = const.

• III. Numerical Solution

– p. 18/31

Well-Balancing Topography Terms - F-Wave Method

(Bale–LeVeque–Mitran–Rossmanith, 2002)

Idea: Concentrate source term at interfaces → Ψ b

i+1/2 and incorporate topography

contribution Ψ b

i+1/2 ∆x into the Riemann solution. Now we decompose:

∆ ˜ F−Ψ b

i+1/2 ∆x = 4 k=1 ζk ˆ

rk . Then, same algorithm with f-waves Zk = ζk ˆ rk and speeds sk = ˆ λk. The interface source term Ψ b

i+1/2 must satisfy the discrete steady state condition

∆ ˜ F/∆x = Ψ b

i+1/2 ,

whenever initial data correspond to equilibrium at rest. We take Ψ b

i+1/2 ∆x = −(0, g ˆ

hs ∆b, 0, g ˆ hf ∆b)T , ∆b = bi+1 − bi . Then, if initially steady state ⇒ Zk ≡ 0 ⇒ updating formula gives Qn+1

i

= Qn

i

⇒ equilibrium is maintained. If solution close to a steady state, it is the deviation from equilibrium that is decomposed ⇒ perturbations are well modeled.

• III. Numerical Solution

– p. 19/31

Interphase Drag Terms

Consider system with drag source terms: ∂tq + A(q)∂xq = ψ b(q) + ψD(q) , ψD(q) = (0, γ F D, 0, −F D)

T ,

F D = D(hs + hf)(vf − vs) . Drag function: D = D/ρf , D = S1(ϕ; σ) + S2(ϕ; σ)|vf − vs|. σ = parameters, e.g. ρs, ρf, ds (grain diameter). Note: At rest ψD(q) = 0 ⇒ no influence on balance conditions at rest. Fractional Step Method

• 1. Solve over ∆t the system ∂tq + A(q)∂xq − ψb(q) = 0, as described.
• 2. Solve exactly over ∆t the system of ODEs ∂tq = ψD(q).

→ Efficient modeling of both fast and slow velocity relaxation processes.

• III. Numerical Solution

– p. 20/31

Eigenvalues Computation Explicit expression of the eigenvalues not available: need numerical computation.

λ1 λ2 λ3 λ4 min(vf,vs) − a max(vf,vs) + a a = sqrt(g h)

P(λ)

External eigenvalues λ1,4 computed through Newton iteration with starting guess min(vf, vs) − a for λ1 and max(vf, vs) + a for λ4.

For λ ∈ [min(vf, vs) − a, λ1] : P ′(λ) < 0, P ′′(λ) > 0, For λ ∈ [λ4, max(vf, vs) + a] : P ′(λ) > 0, P ′′(λ) > 0.

• Known λ1,4: Vieta’s formulas to obtain the internal eigenvalues λ2,3.
• Explicit expressions of right eigenvectors rk and left eigenvectors lk in

terms of λk, k = 1, . . . , 4.

• III. Numerical Solution

– p. 21/31

Initial flow hump with higher fluid content

−15 −10 −5 5 10 15 0.9 1 1.1

Flow height at t = 0

−15 −10 −5 5 10 15 0.4 0.45 0.5 0.55 0.6

Solid volume fraction at t = 0

−15 −10 −5 5 10 15 −0.2 −0.1 0.1 0.2

Phase velocities at t =0

solid fluid

Grid cells = 1000, 2nd order (MC limiter), CFL = 0.9. III. Numerical Solution

– p. 22/31

Initial flow hump with higher fluid content

−15 −10 −5 5 10 15 0.9 1 1.1

Flow height at t = 0.5

−15 −10 −5 5 10 15 0.4 0.45 0.5 0.55 0.6

Solid volume fraction at t = 0.5

−15 −10 −5 5 10 15 −0.2 −0.1 0.1 0.2

Phase velocities at t =0.5

solid fluid

• III. Numerical Solution

– p. 22/31

Initial flow hump with higher fluid content

−15 −10 −5 5 10 15 0.9 1 1.1

Flow height at t = 1.5

−15 −10 −5 5 10 15 0.4 0.45 0.5 0.55 0.6

Solid volume fraction at t = 1.5

−15 −10 −5 5 10 15 −0.2 −0.1 0.1 0.2

Phase velocities at t =1.5

solid fluid

• III. Numerical Solution

– p. 22/31

Initial flow hump with higher fluid content

−15 −10 −5 5 10 15 0.9 1 1.1

Flow height at t = 2.5

−15 −10 −5 5 10 15 0.4 0.45 0.5 0.55 0.6

Solid volume fraction at t = 2.5

−15 −10 −5 5 10 15 −0.2 −0.1 0.1 0.2

Phase velocities at t =2.5

solid fluid

• III. Numerical Solution

– p. 22/31

Initial flow hump with higher fluid content

−15 −10 −5 5 10 15 0.9 1 1.1

Flow height at t = 3.5

−15 −10 −5 5 10 15 0.4 0.45 0.5 0.55 0.6

Solid volume fraction at t = 3.5

−15 −10 −5 5 10 15 −0.2 −0.1 0.1 0.2

Phase velocities at t =3.5

solid fluid

• III. Numerical Solution

– p. 22/31

Only initial variation of h Only initial variation of ϕ

−15 −10 −5 5 10 15 0.95 1 1.05 1.1 1.15

Flow depth at t = 3.5

−15 −10 −5 5 10 15 0.4 0.45 0.5 0.55 0.6

Solid volume fraction at t = 3.5

−15 −10 −5 5 10 15 −0.2 −0.1 0.1 0.2

Phase velocities at t =3.5

vs vf −15 −10 −5 5 10 15 0.95 1 1.05 1.1 1.15

Flow depth at t = 3.5

−15 −10 −5 5 10 15 0.4 0.45 0.5 0.55 0.6

Solid volume fraction at t = 3.5

−15 −10 −5 5 10 15 −0.2 −0.1 0.1 0.2

Phase velocities at t =3.5

vs vf

• III. Numerical Solution

– p. 23/31

Numerical Test: Perturbation of a steady state at rest Extension of LeVeque’s classical test [JCP, vol. 146, 1998] b(x) =

• 0.25(cos(π(x − 0.5)/0.1) + 1)

if |x − 0.5| < 0.1 ,

• therwise.

−0.8 −0.6 −0.4 −0.2 0.2 0.4 0.6 0.8 1 1

h + b at t = 0

−0.8 −0.6 −0.4 −0.2 0.2 0.4 0.6 0.8 1 0.6

φ at t = 0

For −0.6 < x < −0.5 : h(x, 0) = h0 + ˜ h and ϕ(x, 0) = ϕ0 − ˜ ϕ . h0 = 1 , ϕ0 = 0.6 , ˜ h = ˜ ϕ = 10−3 . Grid cells = 100, 2nd order, CFL = 0.9. Reference curve: Grid cells = 1000.

• III. Numerical Solution

– p. 24/31

Perturbation of a steady state at rest

−0.8 −0.6 −0.4 −0.2 0.2 0.4 0.6 0.8 1 0.9998 0.9999 1 1.0001 1.0002 1.0003 1.0004 1.0005

h + b at t = 0

−0.8 −0.6 −0.4 −0.2 0.2 0.4 0.6 0.8 1 0.5992 0.5994 0.5996 0.5998 0.6 0.6002 0.6004

φ at t = 0

• III. Numerical Solution

– p. 25/31

Perturbation of a steady state at rest

−0.8 −0.6 −0.4 −0.2 0.2 0.4 0.6 0.8 1 0.9998 0.9999 1 1.0001 1.0002 1.0003 1.0004 1.0005

h + b at t = 0.25

−0.8 −0.6 −0.4 −0.2 0.2 0.4 0.6 0.8 1 0.5992 0.5994 0.5996 0.5998 0.6 0.6002 0.6004

φ at t = 0.25

• III. Numerical Solution

– p. 25/31

Perturbation of a steady state at rest

−0.8 −0.6 −0.4 −0.2 0.2 0.4 0.6 0.8 1 0.9998 0.9999 1 1.0001 1.0002 1.0003 1.0004 1.0005

h + b at t = 0.5

−0.8 −0.6 −0.4 −0.2 0.2 0.4 0.6 0.8 1 0.5992 0.5994 0.5996 0.5998 0.6 0.6002 0.6004

φ at t = 0.5

• III. Numerical Solution

– p. 25/31

Perturbation of a steady state at rest

−0.8 −0.6 −0.4 −0.2 0.2 0.4 0.6 0.8 1 0.9998 0.9999 1 1.0001 1.0002 1.0003 1.0004 1.0005

h + b at t = 0.75

−0.8 −0.6 −0.4 −0.2 0.2 0.4 0.6 0.8 1 0.5992 0.5994 0.5996 0.5998 0.6 0.6002 0.6004

φ at t = 0.75

• III. Numerical Solution

– p. 25/31

Perturbation of a steady state at rest

−0.8 −0.6 −0.4 −0.2 0.2 0.4 0.6 0.8 1 0.9998 0.9999 1 1.0001 1.0002 1.0003 1.0004 1.0005

h + b at t = 1

−0.8 −0.6 −0.4 −0.2 0.2 0.4 0.6 0.8 1 0.5992 0.5994 0.5996 0.5998 0.6 0.6002 0.6004

φ at t = 1

• III. Numerical Solution

– p. 25/31

Perturbation of a steady state at rest

−0.8 −0.6 −0.4 −0.2 0.2 0.4 0.6 0.8 1 0.9998 0.9999 1 1.0001 1.0002 1.0003 1.0004 1.0005

h + b at t = 1.25

−0.8 −0.6 −0.4 −0.2 0.2 0.4 0.6 0.8 1 0.5992 0.5994 0.5996 0.5998 0.6 0.6002 0.6004

φ at t = 1.25

• III. Numerical Solution

– p. 25/31

Perturbation of a steady state at rest

−0.8 −0.6 −0.4 −0.2 0.2 0.4 0.6 0.8 1 0.9998 0.9999 1 1.0001 1.0002 1.0003 1.0004 1.0005

h + b at t = 1.5

−0.8 −0.6 −0.4 −0.2 0.2 0.4 0.6 0.8 1 0.5992 0.5994 0.5996 0.5998 0.6 0.6002 0.6004

φ at t = 1.5

• III. Numerical Solution

– p. 25/31

Perturbation of a steady state at rest

−0.8 −0.6 −0.4 −0.2 0.2 0.4 0.6 0.8 1 0.9998 0.9999 1 1.0001 1.0002 1.0003 1.0004 1.0005

h + b at t = 1.75

−0.8 −0.6 −0.4 −0.2 0.2 0.4 0.6 0.8 1 0.5992 0.5994 0.5996 0.5998 0.6 0.6002 0.6004

φ at t = 1.75

• III. Numerical Solution

– p. 25/31

Perturbation of a steady state at rest

−0.8 −0.6 −0.4 −0.2 0.2 0.4 0.6 0.8 1 0.9998 0.9999 1 1.0001 1.0002 1.0003 1.0004 1.0005

h + b at t = 2

−0.8 −0.6 −0.4 −0.2 0.2 0.4 0.6 0.8 1 0.5992 0.5994 0.5996 0.5998 0.6 0.6002 0.6004

φ at t = 2

• III. Numerical Solution

– p. 25/31

Perturbation of a steady state at rest

−0.8 −0.6 −0.4 −0.2 0.2 0.4 0.6 0.8 1 0.9998 0.9999 1 1.0001 1.0002 1.0003 1.0004 1.0005

h + b at t = 2.25

−0.8 −0.6 −0.4 −0.2 0.2 0.4 0.6 0.8 1 0.5992 0.5994 0.5996 0.5998 0.6 0.6002 0.6004

φ at t = 2.25

• III. Numerical Solution

– p. 25/31

Perturbation of a steady state at rest

−0.8 −0.6 −0.4 −0.2 0.2 0.4 0.6 0.8 1 0.9998 0.9999 1 1.0001 1.0002 1.0003 1.0004 1.0005

h + b at t = 2.5

−0.8 −0.6 −0.4 −0.2 0.2 0.4 0.6 0.8 1 0.5992 0.5994 0.5996 0.5998 0.6 0.6002 0.6004

φ at t = 2.5

• III. Numerical Solution

– p. 25/31

Perturbation of a steady state at rest

−0.8 −0.6 −0.4 −0.2 0.2 0.4 0.6 0.8 1 0.9998 0.9999 1 1.0001 1.0002 1.0003 1.0004 1.0005

h + b at t = 2.75

−0.8 −0.6 −0.4 −0.2 0.2 0.4 0.6 0.8 1 0.5992 0.5994 0.5996 0.5998 0.6 0.6002 0.6004

φ at t = 2.75

• III. Numerical Solution

– p. 25/31

Perturbation of a steady state at rest

−0.8 −0.6 −0.4 −0.2 0.2 0.4 0.6 0.8 1 0.9998 0.9999 1 1.0001 1.0002 1.0003 1.0004 1.0005

h + b at t = 3

−0.8 −0.6 −0.4 −0.2 0.2 0.4 0.6 0.8 1 0.5992 0.5994 0.5996 0.5998 0.6 0.6002 0.6004

φ at t = 3

• III. Numerical Solution

– p. 25/31

Perturbation of a steady state at rest

−0.8 −0.6 −0.4 −0.2 0.2 0.4 0.6 0.8 1 0.9998 0.9999 1 1.0001 1.0002 1.0003 1.0004 1.0005

h + b at t = 3.25

−0.8 −0.6 −0.4 −0.2 0.2 0.4 0.6 0.8 1 0.5992 0.5994 0.5996 0.5998 0.6 0.6002 0.6004

φ at t = 3.25

• III. Numerical Solution

– p. 25/31

Perturbation of a steady state at rest

−0.8 −0.6 −0.4 −0.2 0.2 0.4 0.6 0.8 1 0.9998 0.9999 1 1.0001 1.0002 1.0003 1.0004 1.0005

h + b at t = 3.5

−0.8 −0.6 −0.4 −0.2 0.2 0.4 0.6 0.8 1 0.5992 0.5994 0.5996 0.5998 0.6 0.6002 0.6004

φ at t = 3.5

• III. Numerical Solution

– p. 25/31

Perturbation of a steady state at rest

−0.8 −0.6 −0.4 −0.2 0.2 0.4 0.6 0.8 1 0.9998 0.9999 1 1.0001 1.0002 1.0003 1.0004 1.0005

h + b at t = 3.75

−0.8 −0.6 −0.4 −0.2 0.2 0.4 0.6 0.8 1 0.5992 0.5994 0.5996 0.5998 0.6 0.6002 0.6004

φ at t = 3.75

• III. Numerical Solution

– p. 25/31

Perturbation of a steady state at rest

−0.8 −0.6 −0.4 −0.2 0.2 0.4 0.6 0.8 1 0.9998 0.9999 1 1.0001 1.0002 1.0003 1.0004 1.0005

h + b at t = 4

−0.8 −0.6 −0.4 −0.2 0.2 0.4 0.6 0.8 1 0.5992 0.5994 0.5996 0.5998 0.6 0.6002 0.6004

φ at t = 4

• III. Numerical Solution

– p. 25/31

Perturbation of a steady state at rest

−0.8 −0.6 −0.4 −0.2 0.2 0.4 0.6 0.8 1 0.9998 0.9999 1 1.0001 1.0002 1.0003 1.0004 1.0005

h + b at t = 4.25

−0.8 −0.6 −0.4 −0.2 0.2 0.4 0.6 0.8 1 0.5992 0.5994 0.5996 0.5998 0.6 0.6002 0.6004

φ at t = 4.25

• III. Numerical Solution

– p. 25/31

Perturbation of a steady state at rest

−0.8 −0.6 −0.4 −0.2 0.2 0.4 0.6 0.8 1 0.9998 0.9999 1 1.0001 1.0002 1.0003 1.0004 1.0005

h + b at t = 4.5

−0.8 −0.6 −0.4 −0.2 0.2 0.4 0.6 0.8 1 0.5992 0.5994 0.5996 0.5998 0.6 0.6002 0.6004

φ at t = 4.5

• III. Numerical Solution

– p. 25/31

Perturbation of a steady state at rest

−0.8 −0.6 −0.4 −0.2 0.2 0.4 0.6 0.8 1 0.9998 0.9999 1 1.0001 1.0002 1.0003 1.0004 1.0005

h + b at t = 4.75

−0.8 −0.6 −0.4 −0.2 0.2 0.4 0.6 0.8 1 0.5992 0.5994 0.5996 0.5998 0.6 0.6002 0.6004

φ at t = 4.75

• III. Numerical Solution

– p. 25/31

Perturbation of a steady state at rest

−0.8 −0.6 −0.4 −0.2 0.2 0.4 0.6 0.8 1 0.9998 0.9999 1 1.0001 1.0002 1.0003 1.0004 1.0005

h + b at t = 5

−0.8 −0.6 −0.4 −0.2 0.2 0.4 0.6 0.8 1 0.5992 0.5994 0.5996 0.5998 0.6 0.6002 0.6004

φ at t = 5

• III. Numerical Solution

– p. 25/31

Perturbation of a steady state at rest

−0.8 −0.6 −0.4 −0.2 0.2 0.4 0.6 0.8 1 0.9998 0.9999 1 1.0001 1.0002 1.0003 1.0004 1.0005

h + b at t = 5.25

−0.8 −0.6 −0.4 −0.2 0.2 0.4 0.6 0.8 1 0.5992 0.5994 0.5996 0.5998 0.6 0.6002 0.6004

φ at t = 5.25

• III. Numerical Solution

– p. 25/31

Perturbation of a steady state at rest

−0.8 −0.6 −0.4 −0.2 0.2 0.4 0.6 0.8 1 0.9998 0.9999 1 1.0001 1.0002 1.0003 1.0004 1.0005

h + b at t = 5.5

−0.8 −0.6 −0.4 −0.2 0.2 0.4 0.6 0.8 1 0.5992 0.5994 0.5996 0.5998 0.6 0.6002 0.6004

φ at t = 5.5

• III. Numerical Solution

– p. 25/31

Perturbation of a steady state at rest

−0.8 −0.6 −0.4 −0.2 0.2 0.4 0.6 0.8 1 0.9998 0.9999 1 1.0001 1.0002 1.0003 1.0004 1.0005

h + b at t = 5.75

−0.8 −0.6 −0.4 −0.2 0.2 0.4 0.6 0.8 1 0.5992 0.5994 0.5996 0.5998 0.6 0.6002 0.6004

φ at t = 5.75

• III. Numerical Solution

– p. 25/31

Perturbation of a steady state at rest

−0.8 −0.6 −0.4 −0.2 0.2 0.4 0.6 0.8 1 0.9998 0.9999 1 1.0001 1.0002 1.0003 1.0004 1.0005

h + b at t = 6

−0.8 −0.6 −0.4 −0.2 0.2 0.4 0.6 0.8 1 0.5992 0.5994 0.5996 0.5998 0.6 0.6002 0.6004

φ at t = 6

• III. Numerical Solution

– p. 25/31

Perturbation of a steady state at rest

−0.8 −0.6 −0.4 −0.2 0.2 0.4 0.6 0.8 1 0.9998 0.9999 1 1.0001 1.0002 1.0003 1.0004 1.0005

h + b at t = 6.5

−0.8 −0.6 −0.4 −0.2 0.2 0.4 0.6 0.8 1 0.5992 0.5994 0.5996 0.5998 0.6 0.6002 0.6004

φ at t = 6.5

• III. Numerical Solution

– p. 25/31

Perturbation of a steady state at rest

−0.8 −0.6 −0.4 −0.2 0.2 0.4 0.6 0.8 1 0.9998 0.9999 1 1.0001 1.0002 1.0003 1.0004 1.0005

h + b at t = 7

−0.8 −0.6 −0.4 −0.2 0.2 0.4 0.6 0.8 1 0.5992 0.5994 0.5996 0.5998 0.6 0.6002 0.6004

φ at t = 7

• III. Numerical Solution

– p. 25/31

Perturbation of a steady state at rest

−0.8 −0.6 −0.4 −0.2 0.2 0.4 0.6 0.8 1 0.9998 0.9999 1 1.0001 1.0002 1.0003 1.0004 1.0005

h + b at t = 10

−0.8 −0.6 −0.4 −0.2 0.2 0.4 0.6 0.8 1 0.5992 0.5994 0.5996 0.5998 0.6 0.6002 0.6004

φ at t = 10

• III. Numerical Solution

– p. 25/31

Numerical Test: Perturbation of a steady flow moving over a bump Steady state conditions for a moving flow with vs = vf ≡ v: ϕ = const. , hv = const. , g(h + b) + 1 2v2 = const. Test 1: Convergence to a steady subcritical flow over a bump (as single-phase s.w.)

b(x)= 8 < : 0.2−0.05(x−10)2 if 8<x<12,

• therwise.

I.C. h = 2 , ϕ = const., vs =vf =0. B.C. (hv)in = 4.42, hout = 2 . 100 grid cells; solution at t = 50.

−5 5 10 15 20 25 0.5 1 1.5 2 2.5

h + b

computed exact

Now: take initial disturbance of ϕ. ϕ(x, 0) = ϕ0 + ˜ ϕ , ϕ0 = 0.6 , ˜ ϕ = 10−3 , for − 3.5 ≤ x ≤ −2.5 .

Grid cells = 150, Reference curve: 1500 cells. 2nd order; CFL = 0.9.

• III. Numerical Solution

– p. 26/31

Perturbation of a steady flow in motion

−5 5 10 15 20 25 −4 −2 2 4 x 10

−4

h − hsteady at t = 0

−5 5 10 15 20 25 0.5995 0.6 0.6005 0.601

φ at t = 0

• III. Numerical Solution

– p. 27/31

Perturbation of a steady flow in motion

−5 5 10 15 20 25 −4 −2 2 4 x 10

−4

h − hsteady at t = 1

−5 5 10 15 20 25 0.5995 0.6 0.6005 0.601

φ at t = 1

• III. Numerical Solution

– p. 27/31

Perturbation of a steady flow in motion

−5 5 10 15 20 25 −4 −2 2 4 x 10

−4

h − hsteady at t = 2

−5 5 10 15 20 25 0.5995 0.6 0.6005 0.601

φ at t = 2

• III. Numerical Solution

– p. 27/31

Perturbation of a steady flow in motion

−5 5 10 15 20 25 −4 −2 2 4 x 10

−4

h − hsteady at t = 3

−5 5 10 15 20 25 0.5995 0.6 0.6005 0.601

φ at t = 3

• III. Numerical Solution

– p. 27/31

Perturbation of a steady flow in motion

−5 5 10 15 20 25 −4 −2 2 4 x 10

−4

h − hsteady at t = 4

−5 5 10 15 20 25 0.5995 0.6 0.6005 0.601

φ at t = 4

• III. Numerical Solution

– p. 27/31

Perturbation of a steady flow in motion

−5 5 10 15 20 25 −4 −2 2 4 x 10

−4

h − hsteady at t = 5

−5 5 10 15 20 25 0.5995 0.6 0.6005 0.601

φ at t = 5

• III. Numerical Solution

– p. 27/31

Perturbation of a steady flow in motion

−5 5 10 15 20 25 −4 −2 2 4 x 10

−4

h − hsteady at t = 6

−5 5 10 15 20 25 0.5995 0.6 0.6005 0.601

φ at t = 6

• III. Numerical Solution

– p. 27/31

Perturbation of a steady flow in motion

−5 5 10 15 20 25 −4 −2 2 4 x 10

−4

h − hsteady at t = 7

−5 5 10 15 20 25 0.5995 0.6 0.6005 0.601

φ at t = 7

• III. Numerical Solution

– p. 27/31

Perturbation of a steady flow in motion

−5 5 10 15 20 25 −4 −2 2 4 x 10

−4

h − hsteady at t = 8

−5 5 10 15 20 25 0.5995 0.6 0.6005 0.601

φ at t = 8

• III. Numerical Solution

– p. 27/31

Perturbation of a steady flow in motion

−5 5 10 15 20 25 −4 −2 2 4 x 10

−4

h − hsteady at t = 9

−5 5 10 15 20 25 0.5995 0.6 0.6005 0.601

φ at t = 9

• III. Numerical Solution

– p. 27/31

Perturbation of a steady flow in motion

−5 5 10 15 20 25 −4 −2 2 4 x 10

−4

h − hsteady at t = 10

−5 5 10 15 20 25 0.5995 0.6 0.6005 0.601

φ at t = 10

• III. Numerical Solution

– p. 27/31

Perturbation of a steady flow in motion

−5 5 10 15 20 25 −4 −2 2 4 x 10

−4

h − hsteady at t = 11

−5 5 10 15 20 25 0.5995 0.6 0.6005 0.601

φ at t = 11

• III. Numerical Solution

– p. 27/31

Perturbation of a steady flow in motion

−5 5 10 15 20 25 −4 −2 2 4 x 10

−4

h − hsteady at t = 12

−5 5 10 15 20 25 0.5995 0.6 0.6005 0.601

φ at t = 12

• III. Numerical Solution

– p. 27/31

Perturbation of a steady flow in motion

−5 5 10 15 20 25 −4 −2 2 4 x 10

−4

h − hsteady at t = 13

−5 5 10 15 20 25 0.5995 0.6 0.6005 0.601

φ at t = 13

• III. Numerical Solution

– p. 27/31

Perturbation of a steady flow in motion

−5 5 10 15 20 25 −4 −2 2 4 x 10

−4

h − hsteady at t = 14

−5 5 10 15 20 25 0.5995 0.6 0.6005 0.601

φ at t = 14

• III. Numerical Solution

– p. 27/31

Perturbation of a steady flow in motion

−5 5 10 15 20 25 −4 −2 2 4 x 10

−4

h − hsteady at t = 15

−5 5 10 15 20 25 0.5995 0.6 0.6005 0.601

φ at t = 15

• III. Numerical Solution

– p. 27/31

Perturbation of a steady flow in motion

−5 5 10 15 20 25 −4 −2 2 4 x 10

−4

h − hsteady at t = 16

−5 5 10 15 20 25 0.5995 0.6 0.6005 0.601

φ at t = 16

• III. Numerical Solution

– p. 27/31

Perturbation of a steady flow in motion

−5 5 10 15 20 25 −4 −2 2 4 x 10

−4

h − hsteady at t = 17

−5 5 10 15 20 25 0.5995 0.6 0.6005 0.601

φ at t = 17

• III. Numerical Solution

– p. 27/31

Perturbation of a steady flow in motion

−5 5 10 15 20 25 −4 −2 2 4 x 10

−4

h − hsteady at t = 18

−5 5 10 15 20 25 0.5995 0.6 0.6005 0.601

φ at t = 18

• III. Numerical Solution

– p. 27/31

Perturbation of a steady flow in motion

−5 5 10 15 20 25 −4 −2 2 4 x 10

−4

h − hsteady at t = 19

−5 5 10 15 20 25 0.5995 0.6 0.6005 0.601

φ at t = 19

• III. Numerical Solution

– p. 27/31

Perturbation of a steady flow in motion

−5 5 10 15 20 25 −4 −2 2 4 x 10

−4

h − hsteady at t = 20

−5 5 10 15 20 25 0.5995 0.6 0.6005 0.601

φ at t = 20

• III. Numerical Solution

– p. 27/31

Perturbation of a steady flow in motion

−5 5 10 15 20 25 −4 −2 2 4 x 10

−4

h − hsteady at t = 21

−5 5 10 15 20 25 0.5995 0.6 0.6005 0.601

φ at t = 21

• III. Numerical Solution

– p. 27/31

Perturbation of a steady flow in motion

−5 5 10 15 20 25 −4 −2 2 4 x 10

−4

h − hsteady at t = 22

−5 5 10 15 20 25 0.5995 0.6 0.6005 0.601

φ at t = 22

• III. Numerical Solution

– p. 27/31

Perturbation of a steady flow in motion

−5 5 10 15 20 25 −4 −2 2 4 x 10

−4

h − hsteady at t = 23

−5 5 10 15 20 25 0.5995 0.6 0.6005 0.601

φ at t = 23

• III. Numerical Solution

– p. 27/31

Perturbation of a steady flow in motion

−5 5 10 15 20 25 −4 −2 2 4 x 10

−4

h − hsteady at t = 38

−5 5 10 15 20 25 0.5995 0.6 0.6005 0.601

φ at t = 38

• III. Numerical Solution

– p. 27/31

Numerical experiments with drag Flow hump with higher fluid content

No drag Drag included

−15 −10 −5 5 10 15 0.95 1 1.05 1.1 1.15

Flow depth at t = 0.5

−15 −10 −5 5 10 15 0.4 0.45 0.5 0.55 0.6

Solid volume fraction at t = 0.5

−15 −10 −5 5 10 15 −0.2 −0.1 0.1 0.2

Phase velocities at t =0.5

vs vf −15 −10 −5 5 10 15 0.95 1 1.05 1.1 1.15

Flow depth at t = 0.5

−15 −10 −5 5 10 15 0.4 0.45 0.5 0.55 0.6

Solid volume fraction at t = 0.5

−15 −10 −5 5 10 15 −0.2 −0.1 0.1 0.2

Phase velocities at t =0.5

vs vf

• III. Numerical Solution

– p. 28/31

Numerical experiments with drag Flow hump with higher fluid content

No drag Drag included

−15 −10 −5 5 10 15 0.95 1 1.05 1.1 1.15

Flow depth at t = 1.5

−15 −10 −5 5 10 15 0.4 0.45 0.5 0.55 0.6

Solid volume fraction at t = 1.5

−15 −10 −5 5 10 15 −0.2 −0.1 0.1 0.2

Phase velocities at t =1.5

vs vf −15 −10 −5 5 10 15 0.95 1 1.05 1.1 1.15

Flow depth at t = 1.5

−15 −10 −5 5 10 15 0.4 0.45 0.5 0.55 0.6

Solid volume fraction at t = 1.5

−15 −10 −5 5 10 15 −0.2 −0.1 0.1 0.2

Phase velocities at t =1.5

vs vf

• III. Numerical Solution

– p. 28/31

Numerical experiments with drag Flow hump with higher fluid content

No drag Drag included

−15 −10 −5 5 10 15 0.95 1 1.05 1.1 1.15

Flow depth at t = 2.5

−15 −10 −5 5 10 15 0.4 0.45 0.5 0.55 0.6

Solid volume fraction at t = 2.5

−15 −10 −5 5 10 15 −0.2 −0.1 0.1 0.2

Phase velocities at t =2.5

vs vf −15 −10 −5 5 10 15 0.95 1 1.05 1.1 1.15

Flow depth at t = 2.5

−15 −10 −5 5 10 15 0.4 0.45 0.5 0.55 0.6

Solid volume fraction at t = 2.5

−15 −10 −5 5 10 15 −0.2 −0.1 0.1 0.2

Phase velocities at t =2.5

vs vf

• III. Numerical Solution

– p. 28/31

Numerical experiments with drag Flow hump with higher fluid content

No drag Drag included

−15 −10 −5 5 10 15 0.95 1 1.05 1.1 1.15

Flow depth at t = 3.5

−15 −10 −5 5 10 15 0.4 0.45 0.5 0.55 0.6

Solid volume fraction at t = 3.5

−15 −10 −5 5 10 15 −0.2 −0.1 0.1 0.2

Phase velocities at t =3.5

vs vf −15 −10 −5 5 10 15 0.95 1 1.05 1.1 1.15

Flow depth at t = 3.5

−15 −10 −5 5 10 15 0.4 0.45 0.5 0.55 0.6

Solid volume fraction at t = 3.5

−15 −10 −5 5 10 15 −0.2 −0.1 0.1 0.2

Phase velocities at t =3.5

vs vf

• III. Numerical Solution

– p. 28/31

Numerical Experiment: Dam-Break Problem Initially: discontinuity between two constant states with flow at rest (vs =vf =0). Left: hℓ = 3, ϕℓ = 0.7 ; Right: hr = 2, ϕr = 0.4 .

• III. Numerical Solution

– p. 29/31

Numerical Experiment: Dam-Break Problem Initially: discontinuity between two constant states with flow at rest (vs =vf =0). Left: hℓ = 3, ϕℓ = 0.7 ; Right: hr = 2, ϕr = 0.4 . Compare:

• 1. Solution of two-phase model with no drag.
• III. Numerical Solution

– p. 29/31

Numerical Experiment: Dam-Break Problem Initially: discontinuity between two constant states with flow at rest (vs =vf =0). Left: hℓ = 3, ϕℓ = 0.7 ; Right: hr = 2, ϕr = 0.4 . Compare:

• 1. Solution of two-phase model with no drag.
• 2. Solution of two-phase model with drag effects.
• III. Numerical Solution

– p. 29/31

Numerical Experiment: Dam-Break Problem Initially: discontinuity between two constant states with flow at rest (vs =vf =0). Left: hℓ = 3, ϕℓ = 0.7 ; Right: hr = 2, ϕr = 0.4 . Compare:

• 1. Solution of two-phase model with no drag.
• 2. Solution of two-phase model with drag effects.
• 3. Solution of reduced model derived theoretically from two-phase model by

assuming drag strong enough to drive instantaneously phase velocities to equilibrium. ⇒ Hyperbolic system of three equations: ⋆ Mass and momentum conservation for the mixture + advection for ϕ. Riemann problems can be solved exactly.

• III. Numerical Solution

– p. 29/31

Numerical Experiment: Dam-Break Problem Initially: discontinuity between two constant states with flow at rest (vs =vf =0). Left: hℓ = 3, ϕℓ = 0.7 ; Right: hr = 2, ϕr = 0.4 . Compare:

• 1. Solution of two-phase model with no drag.
• 2. Solution of two-phase model with drag effects.
• 3. Solution of reduced model derived theoretically from two-phase model by

assuming drag strong enough to drive instantaneously phase velocities to equilibrium. ⇒ Hyperbolic system of three equations: ⋆ Mass and momentum conservation for the mixture + advection for ϕ. Riemann problems can be solved exactly.

• 4. Solution of two-phase model in the limit of infinitely large drag:

Impose numerically instantaneous velocity equilibrium by setting vs = vf = veq in fractional step. veq = hsvs+γhf vf

hs+γhf

• t0 = limit for t → ∞ of solution of ∂tq = ψD(q).
• III. Numerical Solution

– p. 29/31

Dam-Break Problem

hℓ = 3, ϕℓ = 0.7 ; hr = 2, ϕr = 0.4 .

• 1. No drag contribution

Grid cells = 1000

−5 −4 −3 −2 −1 1 2 3 4 5 2 2.5 3

Flow depth at t = 0.5

−5 −4 −3 −2 −1 1 2 3 4 5 0.4 0.5 0.6 0.7

Solid volume fraction at t = 0.5

−5 −4 −3 −2 −1 1 2 3 4 5 0.5 1 1.5

Phase velocities at t =0.5 vs vf

− − − Exact solution reduced model (instantaneous velocity equilibrium)

• III. Numerical Solution

– p. 30/31

Dam-Break Problem

hℓ = 3, ϕℓ = 0.7 ; hr = 2, ϕr = 0.4 .

• 2. Drag effects included

Grid cells = 1000

−5 −4 −3 −2 −1 1 2 3 4 5 2 2.5 3

Flow depth at t = 0.5

−5 −4 −3 −2 −1 1 2 3 4 5 0.4 0.5 0.6 0.7

Solid volume fraction at t = 0.5

−5 −4 −3 −2 −1 1 2 3 4 5 0.5 1 1.5

Phase velocities at t =0.5 vs vf

− − − Exact solution reduced model (instantaneous velocity equilibrium) · · · No drag

• III. Numerical Solution

– p. 30/31

Dam-Break Problem

hℓ = 3, ϕℓ = 0.7 ; hr = 2, ϕr = 0.4 .

• 2. Drag effects included

Infinitely large drag

Grid cells = 1000

vs = vf = veq in fractional step

−5 −4 −3 −2 −1 1 2 3 4 5 2 2.5 3

Flow depth at t = 0.5

−5 −4 −3 −2 −1 1 2 3 4 5 0.4 0.5 0.6 0.7

Solid volume fraction at t = 0.5

−5 −4 −3 −2 −1 1 2 3 4 5 0.5 1 1.5

Phase velocities at t =0.5 vs vf

−5 −4 −3 −2 −1 1 2 3 4 5 2 2.5 3

Flow depth at t = 0.5

−5 −4 −3 −2 −1 1 2 3 4 5 0.4 0.5 0.6 0.7

Solid volume fraction at t = 0.5

−5 −4 −3 −2 −1 1 2 3 4 5 0.5 1 1.5

Phase velocities at t =0.5 vs vf

− − − Exact solution reduced model (instantaneous velocity equilibrium) · · · No drag veq =

hsvs+γhf vf hs+γhf

= equilibrium velocity.

• III. Numerical Solution

– p. 30/31

Summary A mathematical and numerical two-phase shallow flow model has been presented for grain/fluid mixtures over variable topography. Numerical solution technique: Finite Volume Method based on a Roe-type Riemann Solver, which includes treatment of topography and inter-phase drag terms. This is only a very first step towards the modeling of realistic geophysical flows.

• IV. Summary and Future Work

– p. 31/31

Summary A mathematical and numerical two-phase shallow flow model has been presented for grain/fluid mixtures over variable topography. Numerical solution technique: Finite Volume Method based on a Roe-type Riemann Solver, which includes treatment of topography and inter-phase drag terms. This is only a very first step towards the modeling of realistic geophysical flows. Current Work Major issue: positivity preservation of flow depth and phase volume fractions, to handle dry bed states (h = 0) and/or vanishing of one phase (ϕ = 0, ϕ = 1). Need to guarantee hs, hf ≥ 0 ⇔ h ≥ 0, ϕ ∈ [0, 1]. Further work: friction terms, 2D model, complex topography...

• IV. Summary and Future Work

– p. 31/31

Summary A mathematical and numerical two-phase shallow flow model has been presented for grain/fluid mixtures over variable topography. Numerical solution technique: Finite Volume Method based on a Roe-type Riemann Solver, which includes treatment of topography and inter-phase drag terms. This is only a very first step towards the modeling of realistic geophysical flows. Current Work Major issue: positivity preservation of flow depth and phase volume fractions, to handle dry bed states (h = 0) and/or vanishing of one phase (ϕ = 0, ϕ = 1). Need to guarantee hs, hf ≥ 0 ⇔ h ≥ 0, ϕ ∈ [0, 1]. Further work: friction terms, 2D model, complex topography... ...Thank you for your attention!

• IV. Summary and Future Work

– p. 31/31