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

Numerical Solution • Assume | v s − v f | 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 � q ℓ if x ≤ ¯ x , Riemann problem: ∂ t q + A ( q ) ∂ x q = 0 with I.C. q ( x, 0) = q r if x > ¯ x . Q n +1 i t n +1 Q n i → approximate solution on cell ( x i − 1 / 2 , x i +1 / 2 ) . Discontinuities at cell interfaces ⇒ Riemann problems. t n Q n i 1) At each cell interface x i +1 / 2 between Q n i and Q n i +1 → solve Riemann problem with data Q n i and Q n i +1 . i → Q n +1 2) Use solution of local Riemann problems to update solution Q n . i Riemann Solver: → Set of waves W k and speeds s k representing the (approximate) Riemann solution structure. k W k • ∆q ≡ Q i +1 − Q i = � k s k W k • For conservative systems ∂ t q + ∂ x F ( q ) = 0 : F ( Q i +1 ) −F ( Q i ) = � Wave-Propagation Algorithm: Q n +1 = Q n i − ∆t ∆x ( A + ∆Q i − 1 / 2 + A − ∆Q i +1 / 2 ) , i k ( s k i +1 / 2 ) ± W k fluctuations: A ± ∆Q i +1 / 2 = � s + =max( s, 0) , s − =min( s, 0) . i +1 / 2 , III. Numerical Solution – p. 13/31

Numerical Solution Homogeneous System ( b ( x ) = 0 , D = 0 ) q = ( h s , h s v s , h f , h f v f ) T , ∂ t q + ∂ x f ( q ) + w ( q, ∂ x q ) = 0 , � T , s + g s + g 1 − γ f + g � h s v s , h s v 2 2 h 2 2 h s h f , h f v f , h f v 2 2 h 2 f ( q ) = f T . w ( q, ∂ x q ) = (0 , γgh s ∂ x h f , 0 , gh f ∂ x h s ) ⋆ Solid and fluid mass equations are conservative. ⋆ Mixture momentum equation is conservative: ∂ t m + ∂ x f m ( q ) = 0 , f m ( q ) = f (2) ( q ) + γ f (4) ( q ) + γ g h s h f . m = h s v s + γ h f v f , ∂h f ∂h s ⋆ Non-conservative products γgh s ∂x , gh f ∂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 ∂ t q + A ( q ) ∂ x q = 0 . At each local cell interface x i +1 / 2 between solution values Q i and Q i +1 solve a Riemann problem for a linearized system ∂ t q + ˆ A ( Q i , Q i +1 ) ∂ x q = 0 with initial data Q i and Q i +1 . The Roe matrix ˆ A ( Q i , Q i +1 ) is defined so as to guarantee conservation for the mass of each phase and for the momentum of the mixture: A ( p, :) ( Q i +1 − Q i ) , f ( p ) ( Q i +1 ) − f ( p ) ( Q i ) = ˆ p = 1 , 3 , A (2 , :) + γ ˆ f m ( Q i +1 ) − f m ( Q i ) = ( ˆ A (4 , :) )( Q i +1 − Q i ) . We take ˆ A = A (ˆ h s , ˆ h f , ˆ v s , ˆ v f ) , with the choice � � h θ,i v θ,i + h θ,i +1 v θ,i +1 h θ = h θ,i + h θ,i +1 ˆ ˆ v θ = , θ = s, f . and � � 2 h θ,i + h θ,i +1 Then: waves W k = α k ˆ r k , and speeds s k = ˆ r k , ∆q = � 4 k =1 α k ˆ λ k , k = 1 , . . . , 4 . r k , ˆ λ k } = eigenpairs of ˆ { ˆ A . III. Numerical Solution – p. 15/31

Wave Propagation Algorithms (LeVeque, 1997) – F-Wave Formulation Basic software: CLAWPACK . k W k ; Roe: W k = α k ˆ r k , s k = ˆ 1) Classical Riemann solver: ∆q = � λ k 2) F-wave Approach: For a conservative system ∂ t q + ∂ x F ( q ) = 0 , k Z k . decompose flux jump ∆ F ≡ F ( Q i +1 ) − F ( Q i ) = � Local Riemann solution approximated by f-waves Z k and associated speeds s k . Roe: Z k = ζ k ˆ r k , s k = ˆ r k , ˆ λ k } = eigenpairs of Roe matrix for F ′ ( q ) . λ k , { ˆ Fluctuations � � Z k A + ∆Q i +1 / 2 = Z k A − ∆Q i +1 / 2 = i +1 / 2 , i +1 / 2 . k : s k k : s k i +1 / 2 < 0 i +1 / 2 > 0 Algorithm Q n +1 i − ∆t = Q n ∆x ( A + ∆Q i − 1 / 2 + A − ∆Q i +1 / 2 ) i (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 . k W k ; Roe: W k = α k ˆ r k , s k = ˆ 1) Classical Riemann solver: ∆q = � λ k 2) F-wave Approach: For a conservative system ∂ t q + ∂ x F ( q ) = 0 , k Z k . decompose flux jump ∆ F ≡ F ( Q i +1 ) − F ( Q i ) = � Local Riemann solution approximated by f-waves Z k and associated speeds s k . Roe: Z k = ζ k ˆ r k , s k = ˆ r k , ˆ λ k } = eigenpairs of Roe matrix for F ′ ( q ) . λ k , { ˆ Fluctuations � � Z k A + ∆Q i +1 / 2 = Z k A − ∆Q i +1 / 2 = i +1 / 2 , i +1 / 2 . k : s k k : s k i +1 / 2 < 0 i +1 / 2 > 0 Algorithm (high-resolution) Q n +1 i − ∆t ∆x ( A + ∆Q i − 1 / 2 + A − ∆Q i +1 / 2 ) − ∆t = Q n ∆x ( F c i +1 / 2 − F c i − 1 / 2 ) i 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 Here non-conservative system ∂ t q + ∂ x f ( q ) + w ( q, ∂ x q ) = 0 , Difficulty: T . w ( q, ∂ x q ) = (0 , γgh s ∂ x h f , 0 , gh f ∂ x h s ) 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 Here non-conservative system ∂ t q + ∂ x f ( q ) + w ( q, ∂ x q ) = 0 , Difficulty: T . w ( q, ∂ x q ) = (0 , γgh s ∂ x h f , 0 , gh f ∂ x h s ) 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, ∂ x q ) consistent with Roe linearization and define approximate flux difference h f ∆h s ) T . ∆ ˜ F = ∆f + (0 , γ g ˆ h s ∆h f , 0 , g ˆ Then decompose ∆ ˜ F = � 4 k =1 ζ k ˆ r k , Z k = ζ k ˆ s k = ˆ r k , λ k , k = 1 , . . . , 4 . and set r k , ˆ λ k } = eigenpairs of ˆ { ˆ A . ∆ ˜ F = ˆ A ∆q Note: (classical Roe property). III. Numerical Solution – p. 17/31

Topography Source Terms Consider system with topography terms: q = ( h s , h s v s , h f , h f v f ) T , ∂ t q + A ( q ) ∂ x q = ψ b ( q ) , T , ψ b ( q ) = − (0 , gh s ∂ x b, 0 , gh f ∂ x b ) b = b ( x ) . Need well-balancing: efficient modeling of equilibrium and quasi-equilibrium states → A ( q ) ∂ x q ≈ ψ b ( q ) . Steady states at rest ( v s = v f = 0 ): h f h s + h f + b = const. = const. and h s That is h + b = const. ϕ = const. and 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: ∆ ˜ i +1 / 2 ∆x = � 4 F− Ψ b k =1 ζ k ˆ r k . Then, same algorithm with f-waves Z k = ζ k ˆ r k and speeds s k = ˆ λ 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. h f ∆b ) T , i +1 / 2 ∆x = − (0 , g ˆ h s ∆b, 0 , g ˆ Ψ b ∆b = b i +1 − b i . We take Then, if initially steady state ⇒ Z k ≡ 0 ⇒ updating formula gives Q n +1 = Q n i 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: ∂ t q + A ( q ) ∂ x q = ψ b ( q ) + ψ D ( q ) , T , F D = D ( h s + h f )( v f − v s ) . ψ D ( q ) = (0 , γ F D , 0 , − F D ) Drag function: D = D /ρ f , D = S 1 ( ϕ ; σ ) + S 2 ( ϕ ; σ ) | v f − v s | . σ = parameters, e.g. ρ s , ρ f , d s (grain diameter). Note: At rest ψ D ( q ) = 0 ⇒ no influence on balance conditions at rest. Fractional Step Method 1. Solve over ∆t the system ∂ t q + A ( q ) ∂ x q − ψ b ( q ) = 0 , as described. 2. Solve exactly over ∆t the system of ODEs ∂ t q = ψ 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. P( λ ) λ 1 , 4 External eigenvalues computed through Newton iteration with starting guess min( v f , v s ) − a for λ 1 and max( v f , v s ) + a for λ 4 . λ 1 λ 2 λ 3 λ 4 0 min(v f ,v s ) − a max(v f ,v s ) + a For λ ∈ [min( v f , v s ) − a, λ 1 ] : P ′ ( λ ) < 0 , P ′′ ( λ ) > 0 , For λ ∈ [ λ 4 , max( v f , v s ) + a ] : P ′ ( λ ) > 0 , P ′′ ( λ ) > 0 . a = sqrt(g h) • Known λ 1 , 4 : Vieta’s formulas to obtain the internal eigenvalues λ 2 , 3 . • Explicit expressions of right eigenvectors r k and left eigenvectors l k in terms of λ k , k = 1 , . . . , 4 . III. Numerical Solution – p. 21/31

Initial flow hump with higher fluid content Flow height at t = 0 1.1 1 0.9 −15 −10 −5 0 5 10 15 Solid volume fraction at t = 0 0.6 0.55 0.5 0.45 0.4 −15 −10 −5 0 5 10 15 Phase velocities at t =0 0.2 0.1 0 −0.1 solid −0.2 fluid −15 −10 −5 0 5 10 15 Grid cells = 1000 , 2 nd order (MC limiter), CFL = 0 . 9 . III. Numerical Solution – p. 22/31

Initial flow hump with higher fluid content Flow height at t = 0.5 1.1 1 0.9 −15 −10 −5 0 5 10 15 Solid volume fraction at t = 0.5 0.6 0.55 0.5 0.45 0.4 −15 −10 −5 0 5 10 15 Phase velocities at t =0.5 0.2 0.1 0 −0.1 solid −0.2 fluid −15 −10 −5 0 5 10 15 III. Numerical Solution – p. 22/31

Initial flow hump with higher fluid content Flow height at t = 1.5 1.1 1 0.9 −15 −10 −5 0 5 10 15 Solid volume fraction at t = 1.5 0.6 0.55 0.5 0.45 0.4 −15 −10 −5 0 5 10 15 Phase velocities at t =1.5 0.2 0.1 0 −0.1 solid −0.2 fluid −15 −10 −5 0 5 10 15 III. Numerical Solution – p. 22/31

Initial flow hump with higher fluid content Flow height at t = 2.5 1.1 1 0.9 −15 −10 −5 0 5 10 15 Solid volume fraction at t = 2.5 0.6 0.55 0.5 0.45 0.4 −15 −10 −5 0 5 10 15 Phase velocities at t =2.5 0.2 0.1 0 −0.1 solid −0.2 fluid −15 −10 −5 0 5 10 15 III. Numerical Solution – p. 22/31

Initial flow hump with higher fluid content Flow height at t = 3.5 1.1 1 0.9 −15 −10 −5 0 5 10 15 Solid volume fraction at t = 3.5 0.6 0.55 0.5 0.45 0.4 −15 −10 −5 0 5 10 15 Phase velocities at t =3.5 0.2 0.1 0 −0.1 solid −0.2 fluid −15 −10 −5 0 5 10 15 III. Numerical Solution – p. 22/31

Only initial variation of h Only initial variation of ϕ Flow depth at t = 3.5 Flow depth at t = 3.5 1.15 1.15 1.1 1.1 1.05 1.05 1 1 0.95 0.95 −15 −10 −5 0 5 10 15 −15 −10 −5 0 5 10 15 Solid volume fraction at t = 3.5 Solid volume fraction at t = 3.5 0.6 0.6 0.55 0.55 0.5 0.5 0.45 0.45 0.4 0.4 −15 −10 −5 0 5 10 15 −15 −10 −5 0 5 10 15 Phase velocities at t =3.5 Phase velocities at t =3.5 0.2 0.2 0.1 0.1 0 0 v s v s −0.1 −0.1 v f v f −0.2 −0.2 −15 −10 −5 0 5 10 15 −15 −10 −5 0 5 10 15 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] � 0 . 25(cos( π ( x − 0 . 5) / 0 . 1) + 1) if | x − 0 . 5 | < 0 . 1 , b ( x ) = 0 otherwise . h + b at t = 0 For − 0 . 6 < x < − 0 . 5 : 1 h ( x, 0) = h 0 + ˜ h and 0 ϕ ( x, 0) = ϕ 0 − ˜ ϕ . −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 φ at t = 0 h 0 = 1 , ϕ 0 = 0 . 6 , 0.6 ϕ = 10 − 3 . ˜ h = ˜ 0 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 Grid cells = 100 , 2 nd order, CFL = 0 . 9 . Reference curve: Grid cells = 1000 . III. Numerical Solution – p. 24/31

Perturbation of a steady state at rest h + b at t = 0 1.0005 1.0004 1.0003 1.0002 1.0001 1 0.9999 0.9998 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 φ at t = 0 0.6004 0.6002 0.6 0.5998 0.5996 0.5994 0.5992 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 III. Numerical Solution – p. 25/31

Perturbation of a steady state at rest h + b at t = 0.25 1.0005 1.0004 1.0003 1.0002 1.0001 1 0.9999 0.9998 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 φ at t = 0.25 0.6004 0.6002 0.6 0.5998 0.5996 0.5994 0.5992 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 III. Numerical Solution – p. 25/31

Perturbation of a steady state at rest h + b at t = 0.5 1.0005 1.0004 1.0003 1.0002 1.0001 1 0.9999 0.9998 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 φ at t = 0.5 0.6004 0.6002 0.6 0.5998 0.5996 0.5994 0.5992 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 III. Numerical Solution – p. 25/31

Perturbation of a steady state at rest h + b at t = 0.75 1.0005 1.0004 1.0003 1.0002 1.0001 1 0.9999 0.9998 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 φ at t = 0.75 0.6004 0.6002 0.6 0.5998 0.5996 0.5994 0.5992 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 III. Numerical Solution – p. 25/31

Perturbation of a steady state at rest h + b at t = 1 1.0005 1.0004 1.0003 1.0002 1.0001 1 0.9999 0.9998 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 φ at t = 1 0.6004 0.6002 0.6 0.5998 0.5996 0.5994 0.5992 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 III. Numerical Solution – p. 25/31

Perturbation of a steady state at rest h + b at t = 1.25 1.0005 1.0004 1.0003 1.0002 1.0001 1 0.9999 0.9998 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 φ at t = 1.25 0.6004 0.6002 0.6 0.5998 0.5996 0.5994 0.5992 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 III. Numerical Solution – p. 25/31

Perturbation of a steady state at rest h + b at t = 1.5 1.0005 1.0004 1.0003 1.0002 1.0001 1 0.9999 0.9998 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 φ at t = 1.5 0.6004 0.6002 0.6 0.5998 0.5996 0.5994 0.5992 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 III. Numerical Solution – p. 25/31

Perturbation of a steady state at rest h + b at t = 1.75 1.0005 1.0004 1.0003 1.0002 1.0001 1 0.9999 0.9998 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 φ at t = 1.75 0.6004 0.6002 0.6 0.5998 0.5996 0.5994 0.5992 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 III. Numerical Solution – p. 25/31

Perturbation of a steady state at rest h + b at t = 2 1.0005 1.0004 1.0003 1.0002 1.0001 1 0.9999 0.9998 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 φ at t = 2 0.6004 0.6002 0.6 0.5998 0.5996 0.5994 0.5992 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 III. Numerical Solution – p. 25/31

Perturbation of a steady state at rest h + b at t = 2.25 1.0005 1.0004 1.0003 1.0002 1.0001 1 0.9999 0.9998 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 φ at t = 2.25 0.6004 0.6002 0.6 0.5998 0.5996 0.5994 0.5992 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 III. Numerical Solution – p. 25/31

Perturbation of a steady state at rest h + b at t = 2.5 1.0005 1.0004 1.0003 1.0002 1.0001 1 0.9999 0.9998 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 φ at t = 2.5 0.6004 0.6002 0.6 0.5998 0.5996 0.5994 0.5992 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 III. Numerical Solution – p. 25/31

Perturbation of a steady state at rest h + b at t = 2.75 1.0005 1.0004 1.0003 1.0002 1.0001 1 0.9999 0.9998 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 φ at t = 2.75 0.6004 0.6002 0.6 0.5998 0.5996 0.5994 0.5992 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 III. Numerical Solution – p. 25/31

Perturbation of a steady state at rest h + b at t = 3 1.0005 1.0004 1.0003 1.0002 1.0001 1 0.9999 0.9998 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 φ at t = 3 0.6004 0.6002 0.6 0.5998 0.5996 0.5994 0.5992 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 III. Numerical Solution – p. 25/31

Perturbation of a steady state at rest h + b at t = 3.25 1.0005 1.0004 1.0003 1.0002 1.0001 1 0.9999 0.9998 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 φ at t = 3.25 0.6004 0.6002 0.6 0.5998 0.5996 0.5994 0.5992 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 III. Numerical Solution – p. 25/31

Perturbation of a steady state at rest h + b at t = 3.5 1.0005 1.0004 1.0003 1.0002 1.0001 1 0.9999 0.9998 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 φ at t = 3.5 0.6004 0.6002 0.6 0.5998 0.5996 0.5994 0.5992 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 III. Numerical Solution – p. 25/31

Perturbation of a steady state at rest h + b at t = 3.75 1.0005 1.0004 1.0003 1.0002 1.0001 1 0.9999 0.9998 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 φ at t = 3.75 0.6004 0.6002 0.6 0.5998 0.5996 0.5994 0.5992 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 III. Numerical Solution – p. 25/31

Perturbation of a steady state at rest h + b at t = 4 1.0005 1.0004 1.0003 1.0002 1.0001 1 0.9999 0.9998 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 φ at t = 4 0.6004 0.6002 0.6 0.5998 0.5996 0.5994 0.5992 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 III. Numerical Solution – p. 25/31

Perturbation of a steady state at rest h + b at t = 4.25 1.0005 1.0004 1.0003 1.0002 1.0001 1 0.9999 0.9998 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 φ at t = 4.25 0.6004 0.6002 0.6 0.5998 0.5996 0.5994 0.5992 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 III. Numerical Solution – p. 25/31

Perturbation of a steady state at rest h + b at t = 4.5 1.0005 1.0004 1.0003 1.0002 1.0001 1 0.9999 0.9998 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 φ at t = 4.5 0.6004 0.6002 0.6 0.5998 0.5996 0.5994 0.5992 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 III. Numerical Solution – p. 25/31

Perturbation of a steady state at rest h + b at t = 4.75 1.0005 1.0004 1.0003 1.0002 1.0001 1 0.9999 0.9998 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 φ at t = 4.75 0.6004 0.6002 0.6 0.5998 0.5996 0.5994 0.5992 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 III. Numerical Solution – p. 25/31

Perturbation of a steady state at rest h + b at t = 5 1.0005 1.0004 1.0003 1.0002 1.0001 1 0.9999 0.9998 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 φ at t = 5 0.6004 0.6002 0.6 0.5998 0.5996 0.5994 0.5992 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 III. Numerical Solution – p. 25/31

Perturbation of a steady state at rest h + b at t = 5.25 1.0005 1.0004 1.0003 1.0002 1.0001 1 0.9999 0.9998 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 φ at t = 5.25 0.6004 0.6002 0.6 0.5998 0.5996 0.5994 0.5992 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 III. Numerical Solution – p. 25/31

Perturbation of a steady state at rest h + b at t = 5.5 1.0005 1.0004 1.0003 1.0002 1.0001 1 0.9999 0.9998 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 φ at t = 5.5 0.6004 0.6002 0.6 0.5998 0.5996 0.5994 0.5992 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 III. Numerical Solution – p. 25/31

Perturbation of a steady state at rest h + b at t = 5.75 1.0005 1.0004 1.0003 1.0002 1.0001 1 0.9999 0.9998 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 φ at t = 5.75 0.6004 0.6002 0.6 0.5998 0.5996 0.5994 0.5992 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 III. Numerical Solution – p. 25/31

Perturbation of a steady state at rest h + b at t = 6 1.0005 1.0004 1.0003 1.0002 1.0001 1 0.9999 0.9998 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 φ at t = 6 0.6004 0.6002 0.6 0.5998 0.5996 0.5994 0.5992 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 III. Numerical Solution – p. 25/31

Perturbation of a steady state at rest h + b at t = 6.5 1.0005 1.0004 1.0003 1.0002 1.0001 1 0.9999 0.9998 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 φ at t = 6.5 0.6004 0.6002 0.6 0.5998 0.5996 0.5994 0.5992 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 III. Numerical Solution – p. 25/31

Perturbation of a steady state at rest h + b at t = 7 1.0005 1.0004 1.0003 1.0002 1.0001 1 0.9999 0.9998 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 φ at t = 7 0.6004 0.6002 0.6 0.5998 0.5996 0.5994 0.5992 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 III. Numerical Solution – p. 25/31

Perturbation of a steady state at rest h + b at t = 10 1.0005 1.0004 1.0003 1.0002 1.0001 1 0.9999 0.9998 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 φ at t = 10 0.6004 0.6002 0.6 0.5998 0.5996 0.5994 0.5992 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 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 v s = v f ≡ v : g ( h + b ) + 1 2 v 2 = const. ϕ = const. , hv = const. , Test 1: Convergence to a steady subcritical flow over a bump (as single-phase s.w.) h + b 8 0 . 2 − 0 . 05( x − 10) 2 2.5 if 8 <x< 12 , < b ( x )= 2 0 otherwise . : 1.5 h = 2 , ϕ = const., v s = v f =0 . I.C. 1 ( hv ) in = 4 . 42 , h out = 2 . B.C. 0.5 computed exact 0 100 grid cells; solution at t = 50 . −5 0 5 10 15 20 25 Now: take initial disturbance of ϕ . ϕ = 10 − 3 , ϕ ( x, 0) = ϕ 0 + ˜ ϕ , ϕ 0 = 0 . 6 , ˜ for − 3 . 5 ≤ x ≤ − 2 . 5 . Grid cells = 150 , Reference curve: 1500 cells. 2 nd order; CFL = 0 . 9 . III. Numerical Solution – p. 26/31

Perturbation of a steady flow in motion h − h steady at t = 0 −4 x 10 4 2 0 −2 −4 −5 0 5 10 15 20 25 φ at t = 0 0.601 0.6005 0.6 0.5995 −5 0 5 10 15 20 25 III. Numerical Solution – p. 27/31

Perturbation of a steady flow in motion h − h steady at t = 1 −4 x 10 4 2 0 −2 −4 −5 0 5 10 15 20 25 φ at t = 1 0.601 0.6005 0.6 0.5995 −5 0 5 10 15 20 25 III. Numerical Solution – p. 27/31

Perturbation of a steady flow in motion h − h steady at t = 2 −4 x 10 4 2 0 −2 −4 −5 0 5 10 15 20 25 φ at t = 2 0.601 0.6005 0.6 0.5995 −5 0 5 10 15 20 25 III. Numerical Solution – p. 27/31

Perturbation of a steady flow in motion h − h steady at t = 3 −4 x 10 4 2 0 −2 −4 −5 0 5 10 15 20 25 φ at t = 3 0.601 0.6005 0.6 0.5995 −5 0 5 10 15 20 25 III. Numerical Solution – p. 27/31

Perturbation of a steady flow in motion h − h steady at t = 4 −4 x 10 4 2 0 −2 −4 −5 0 5 10 15 20 25 φ at t = 4 0.601 0.6005 0.6 0.5995 −5 0 5 10 15 20 25 III. Numerical Solution – p. 27/31

Perturbation of a steady flow in motion h − h steady at t = 5 −4 x 10 4 2 0 −2 −4 −5 0 5 10 15 20 25 φ at t = 5 0.601 0.6005 0.6 0.5995 −5 0 5 10 15 20 25 III. Numerical Solution – p. 27/31

Perturbation of a steady flow in motion h − h steady at t = 6 −4 x 10 4 2 0 −2 −4 −5 0 5 10 15 20 25 φ at t = 6 0.601 0.6005 0.6 0.5995 −5 0 5 10 15 20 25 III. Numerical Solution – p. 27/31

Perturbation of a steady flow in motion h − h steady at t = 7 −4 x 10 4 2 0 −2 −4 −5 0 5 10 15 20 25 φ at t = 7 0.601 0.6005 0.6 0.5995 −5 0 5 10 15 20 25 III. Numerical Solution – p. 27/31

Perturbation of a steady flow in motion h − h steady at t = 8 −4 x 10 4 2 0 −2 −4 −5 0 5 10 15 20 25 φ at t = 8 0.601 0.6005 0.6 0.5995 −5 0 5 10 15 20 25 III. Numerical Solution – p. 27/31

Perturbation of a steady flow in motion h − h steady at t = 9 −4 x 10 4 2 0 −2 −4 −5 0 5 10 15 20 25 φ at t = 9 0.601 0.6005 0.6 0.5995 −5 0 5 10 15 20 25 III. Numerical Solution – p. 27/31

Perturbation of a steady flow in motion h − h steady at t = 10 −4 x 10 4 2 0 −2 −4 −5 0 5 10 15 20 25 φ at t = 10 0.601 0.6005 0.6 0.5995 −5 0 5 10 15 20 25 III. Numerical Solution – p. 27/31

Perturbation of a steady flow in motion h − h steady at t = 11 −4 x 10 4 2 0 −2 −4 −5 0 5 10 15 20 25 φ at t = 11 0.601 0.6005 0.6 0.5995 −5 0 5 10 15 20 25 III. Numerical Solution – p. 27/31

Perturbation of a steady flow in motion h − h steady at t = 12 −4 x 10 4 2 0 −2 −4 −5 0 5 10 15 20 25 φ at t = 12 0.601 0.6005 0.6 0.5995 −5 0 5 10 15 20 25 III. Numerical Solution – p. 27/31

Perturbation of a steady flow in motion h − h steady at t = 13 −4 x 10 4 2 0 −2 −4 −5 0 5 10 15 20 25 φ at t = 13 0.601 0.6005 0.6 0.5995 −5 0 5 10 15 20 25 III. Numerical Solution – p. 27/31

Perturbation of a steady flow in motion h − h steady at t = 14 −4 x 10 4 2 0 −2 −4 −5 0 5 10 15 20 25 φ at t = 14 0.601 0.6005 0.6 0.5995 −5 0 5 10 15 20 25 III. Numerical Solution – p. 27/31

Perturbation of a steady flow in motion h − h steady at t = 15 −4 x 10 4 2 0 −2 −4 −5 0 5 10 15 20 25 φ at t = 15 0.601 0.6005 0.6 0.5995 −5 0 5 10 15 20 25 III. Numerical Solution – p. 27/31

Perturbation of a steady flow in motion h − h steady at t = 16 −4 x 10 4 2 0 −2 −4 −5 0 5 10 15 20 25 φ at t = 16 0.601 0.6005 0.6 0.5995 −5 0 5 10 15 20 25 III. Numerical Solution – p. 27/31

Perturbation of a steady flow in motion h − h steady at t = 17 −4 x 10 4 2 0 −2 −4 −5 0 5 10 15 20 25 φ at t = 17 0.601 0.6005 0.6 0.5995 −5 0 5 10 15 20 25 III. Numerical Solution – p. 27/31

Perturbation of a steady flow in motion h − h steady at t = 18 −4 x 10 4 2 0 −2 −4 −5 0 5 10 15 20 25 φ at t = 18 0.601 0.6005 0.6 0.5995 −5 0 5 10 15 20 25 III. Numerical Solution – p. 27/31

Perturbation of a steady flow in motion h − h steady at t = 19 −4 x 10 4 2 0 −2 −4 −5 0 5 10 15 20 25 φ at t = 19 0.601 0.6005 0.6 0.5995 −5 0 5 10 15 20 25 III. Numerical Solution – p. 27/31

Perturbation of a steady flow in motion h − h steady at t = 20 −4 x 10 4 2 0 −2 −4 −5 0 5 10 15 20 25 φ at t = 20 0.601 0.6005 0.6 0.5995 −5 0 5 10 15 20 25 III. Numerical Solution – p. 27/31

Perturbation of a steady flow in motion h − h steady at t = 21 −4 x 10 4 2 0 −2 −4 −5 0 5 10 15 20 25 φ at t = 21 0.601 0.6005 0.6 0.5995 −5 0 5 10 15 20 25 III. Numerical Solution – p. 27/31

Perturbation of a steady flow in motion h − h steady at t = 22 −4 x 10 4 2 0 −2 −4 −5 0 5 10 15 20 25 φ at t = 22 0.601 0.6005 0.6 0.5995 −5 0 5 10 15 20 25 III. Numerical Solution – p. 27/31

Perturbation of a steady flow in motion h − h steady at t = 23 −4 x 10 4 2 0 −2 −4 −5 0 5 10 15 20 25 φ at t = 23 0.601 0.6005 0.6 0.5995 −5 0 5 10 15 20 25 III. Numerical Solution – p. 27/31

Perturbation of a steady flow in motion h − h steady at t = 38 −4 x 10 4 2 0 −2 −4 −5 0 5 10 15 20 25 φ at t = 38 0.601 0.6005 0.6 0.5995 −5 0 5 10 15 20 25 III. Numerical Solution – p. 27/31

Numerical experiments with drag Flow hump with higher fluid content No drag Drag included Flow depth at t = 0.5 Flow depth at t = 0.5 1.15 1.15 1.1 1.1 1.05 1.05 1 1 0.95 0.95 −15 −10 −5 0 5 10 15 −15 −10 −5 0 5 10 15 Solid volume fraction at t = 0.5 Solid volume fraction at t = 0.5 0.6 0.6 0.55 0.55 0.5 0.5 0.45 0.45 0.4 0.4 −15 −10 −5 0 5 10 15 −15 −10 −5 0 5 10 15 Phase velocities at t =0.5 Phase velocities at t =0.5 0.2 0.2 0.1 0.1 0 0 v s v s −0.1 −0.1 v f v f −0.2 −0.2 −15 −10 −5 0 5 10 15 −15 −10 −5 0 5 10 15 III. Numerical Solution – p. 28/31

Numerical experiments with drag Flow hump with higher fluid content No drag Drag included Flow depth at t = 1.5 Flow depth at t = 1.5 1.15 1.15 1.1 1.1 1.05 1.05 1 1 0.95 0.95 −15 −10 −5 0 5 10 15 −15 −10 −5 0 5 10 15 Solid volume fraction at t = 1.5 Solid volume fraction at t = 1.5 0.6 0.6 0.55 0.55 0.5 0.5 0.45 0.45 0.4 0.4 −15 −10 −5 0 5 10 15 −15 −10 −5 0 5 10 15 Phase velocities at t =1.5 Phase velocities at t =1.5 0.2 0.2 0.1 0.1 0 0 v s v s −0.1 −0.1 v f v f −0.2 −0.2 −15 −10 −5 0 5 10 15 −15 −10 −5 0 5 10 15 III. Numerical Solution – p. 28/31

Numerical experiments with drag Flow hump with higher fluid content No drag Drag included Flow depth at t = 2.5 Flow depth at t = 2.5 1.15 1.15 1.1 1.1 1.05 1.05 1 1 0.95 0.95 −15 −10 −5 0 5 10 15 −15 −10 −5 0 5 10 15 Solid volume fraction at t = 2.5 Solid volume fraction at t = 2.5 0.6 0.6 0.55 0.55 0.5 0.5 0.45 0.45 0.4 0.4 −15 −10 −5 0 5 10 15 −15 −10 −5 0 5 10 15 Phase velocities at t =2.5 Phase velocities at t =2.5 0.2 0.2 0.1 0.1 0 0 v s v s −0.1 −0.1 v f v f −0.2 −0.2 −15 −10 −5 0 5 10 15 −15 −10 −5 0 5 10 15 III. Numerical Solution – p. 28/31

Numerical experiments with drag Flow hump with higher fluid content No drag Drag included Flow depth at t = 3.5 Flow depth at t = 3.5 1.15 1.15 1.1 1.1 1.05 1.05 1 1 0.95 0.95 −15 −10 −5 0 5 10 15 −15 −10 −5 0 5 10 15 Solid volume fraction at t = 3.5 Solid volume fraction at t = 3.5 0.6 0.6 0.55 0.55 0.5 0.5 0.45 0.45 0.4 0.4 −15 −10 −5 0 5 10 15 −15 −10 −5 0 5 10 15 Phase velocities at t =3.5 Phase velocities at t =3.5 0.2 0.2 0.1 0.1 0 0 v s v s −0.1 −0.1 v f v f −0.2 −0.2 −15 −10 −5 0 5 10 15 −15 −10 −5 0 5 10 15 III. Numerical Solution – p. 28/31

Numerical Experiment: Dam-Break Problem Initially: discontinuity between two constant states with flow at rest ( v s = v f =0 ). Left: h ℓ = 3 , ϕ ℓ = 0 . 7 ; Right: h r = 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 ( v s = v f =0 ). Left: h ℓ = 3 , ϕ ℓ = 0 . 7 ; Right: h r = 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 ( v s = v f =0 ). Left: h ℓ = 3 , ϕ ℓ = 0 . 7 ; Right: h r = 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 ( v s = v f =0 ). Left: h ℓ = 3 , ϕ ℓ = 0 . 7 ; Right: h r = 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 ( v s = v f =0 ). Left: h ℓ = 3 , ϕ ℓ = 0 . 7 ; Right: h r = 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 v s = v f = v eq in fractional step. � v eq = h s v s + γh f v f t 0 = limit for t → ∞ of solution of ∂ t q = ψ D ( q ) . � h s + γh f � III. Numerical Solution – p. 29/31

Dam-Break Problem h ℓ = 3 , ϕ ℓ = 0 . 7 ; h r = 2 , ϕ r = 0 . 4 . 1. No drag contribution Grid cells = 1000 Flow depth at t = 0.5 3 2.5 2 −5 −4 −3 −2 −1 0 1 2 3 4 5 Solid volume fraction at t = 0.5 0.7 0.6 0.5 0.4 −5 −4 −3 −2 −1 0 1 2 3 4 5 Phase velocities at t =0.5 v s 1.5 v f 1 0.5 0 −5 −4 −3 −2 −1 0 1 2 3 4 5 − − − Exact solution reduced model (instantaneous velocity equilibrium) III. Numerical Solution – p. 30/31

Dam-Break Problem h ℓ = 3 , ϕ ℓ = 0 . 7 ; h r = 2 , ϕ r = 0 . 4 . 2. Drag effects included Grid cells = 1000 Flow depth at t = 0.5 3 2.5 2 −5 −4 −3 −2 −1 0 1 2 3 4 5 Solid volume fraction at t = 0.5 0.7 0.6 0.5 0.4 −5 −4 −3 −2 −1 0 1 2 3 4 5 Phase velocities at t =0.5 v s 1.5 v f 1 0.5 0 −5 −4 −3 −2 −1 0 1 2 3 4 5 − − − Exact solution reduced model (instantaneous velocity equilibrium) · · · No drag III. Numerical Solution – p. 30/31

Dam-Break Problem h ℓ = 3 , ϕ ℓ = 0 . 7 ; h r = 2 , ϕ r = 0 . 4 . 2. Drag effects included Infinitely large drag v s = v f = v eq in fractional step Grid cells = 1000 Flow depth at t = 0.5 Flow depth at t = 0.5 3 3 2.5 2.5 2 2 −5 −4 −3 −2 −1 0 1 2 3 4 5 −5 −4 −3 −2 −1 0 1 2 3 4 5 Solid volume fraction at t = 0.5 Solid volume fraction at t = 0.5 0.7 0.7 0.6 0.6 0.5 0.5 0.4 0.4 −5 −4 −3 −2 −1 0 1 2 3 4 5 −5 −4 −3 −2 −1 0 1 2 3 4 5 Phase velocities at t =0.5 Phase velocities at t =0.5 v s v s 1.5 1.5 v f v f 1 1 0.5 0.5 0 0 −5 −4 −3 −2 −1 0 1 2 3 4 5 −5 −4 −3 −2 −1 0 1 2 3 4 5 − − − Exact solution reduced model (instantaneous velocity equilibrium) · · · No drag h s v s + γh f v f v eq = = equilibrium velocity. h s + γh f 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 h s , h f ≥ 0 ⇔ h ≥ 0 , ϕ ∈ [0 , 1] . Further work: friction terms, 2D model, complex topography... IV. Summary and Future Work – p. 31/31