Derivation of a bedload transport model with viscous effects E. - PowerPoint PPT Presentation

Introduction Derivation of a bedload transport model with viscous effects E. Audusse, L. Boittin, M. Parisot, J. Sainte-Marie Project-team ANGE, Inria; CEREMA; LJLL, UPMC Universit´ e Paris VI; UMR CNRS 7958 May 30, 2017 L.Boittin ( ANGE, Inria ) 30.05.17 1 / 21

Introduction Why should we simulate sediment transport? Predict the evolution of the river topography Estimate sediment accumulation at the bottom of dams, in harbours... Estimate the stability of structures such as canals and bridges with scour Sediment accumulation Scouring L.Boittin ( ANGE, Inria ) 30.05.17 2 / 21

Introduction Different modes of sediment transport Source: http://theses.univ-lyon2.fr/documents/getpart.php?id=lyon2.2008.pintomartins d&part=154405 L.Boittin ( ANGE, Inria ) 30.05.17 3 / 21

Introduction Classical bed load transport model: Exner equation Clear water → Saint-Venant equations Sediment layer → Exner equation Based on mass conservation: ∂ t h 2 + ∇ · q 2 = 0 h 2 : thickness of sediment layer Solid discharge q 2 given by empirical relationships Einstein (1942), Meyer-Peter and M¨ uller (1948), Nielsen (1992) q 2 depends on the water depth h 1 and velocity u 1 : no intrinsic mechanism in the sediment L.Boittin ( ANGE, Inria ) 30.05.17 4 / 21

Introduction Goals: derive and simulate a new bedload transport model Formal derivation of a model, from the Navier-Stokes equations Coupled model with an energy equation E.D. Fern´ andez-Nieto et al. “Formal deduction of the Saint-Venant-Exner model including arbitrarily sloping sediment beds and associated energy”. In: Mathematical Modelling and Numerical Analysis (2016) Try to include the sediment rheology in the model L.Boittin ( ANGE, Inria ) 30.05.17 5 / 21

Model derivation kinematic condition η H 1 , U 1 , W 1 h 1 , ¯ u 1 no penetration κ ζ ζ kinematic condition H 2 , U 2 , W 2 h 2 , ¯ u 2 B no penetration κ B Saint-Venant Navier-Stokes equations, · thin layer approximation equations sediment transport · high viscosity, high friction model in the sediment layer L.Boittin ( ANGE, Inria ) 30.05.17 6 / 21

Model derivation Inviscid model with high friction Assume the following scaling: H L = ε W = ε U , 1 Re 1 = 1 Fr = 1 , ε , Re 2 = 1 , K ζ = ε, K B = 1 . Assume that the space variations of the µ k are of the order ε 2 . Then, for ε small enough, the system � ∂ t h 1 + ∇ x · ( h 1 ¯ u 1 ) = 0 , h 2 2 Fr 2 Id) = − h 1 Fr 2 ∇ x ( h 2 + B ) − 1 ∂ t ( h 1 ¯ u 1 ) + ∇ x · ( h 1 ¯ u 1 ⊗ ¯ u 1 + r κ ζ (¯ u 1 − ε ˜ u 2 ) , 1 � ∂ t h 2 + ε ∇ x · ( h 2 ˜ u 2 ) = 0 , u 2 ) = − h 2 κ B (˜ Fr 2 ∇ x ( h 2 + rh 1 + B ) + κ ζ (¯ u 1 − ε ˜ u 2 ) , with κ B (˜ u 2 ) = ˜ κ B ˜ u 2 is derived from the bilayer Navier-Stokes system with the following modeling errors: |H 2 − h 2 | = O ( ε 2 ) , |H 1 − h 1 | = O ( ε ) , u 2 | = O ( ε 2 ) . |U 1 − ¯ u 1 | = O ( ε ) , |U 2 − ε ˜ L.Boittin ( ANGE, Inria ) 30.05.17 7 / 21

Model derivation Proof Momentum eq. and boundary conditions: ∂ z ( µ 2 ∂ z U 2 ) = O ( ε 2 ) ,   µ 2 ∂ z U 2 = O ( ε 2 ) ,  ⇒ µ 2 ∂ z U 2 | B = O ( ε 2 ) at z = ζ µ 2 ∂ z U 2 = εκ B U 2 + O ( ε 2 ) , at z = B Imposes that U 2 = ¯ U 2 + O ( ε 2 ) = ε ˜ U 2 + O ( ε 2 ). Fine, but similar to an Exner model, and no rheology... L.Boittin ( ANGE, Inria ) 30.05.17 8 / 21

Model derivation Viscous model Assume the following scaling: H / L = ε W = ε U , 1 Re 1 = 1 Fr = 1 , ε , Re 2 = ε, K ζ = ε, K B = 1 . Assume that the space variations of the µ k are of the order ε 2 . Then, for ε small enough, the system � ∂ t h 1 + ∇ x · ( h 1 ¯ u 1 ) = 0 , h 2 2 Fr 2 Id) = − h 1 Fr 2 ∇ x ( h 2 + B ) − 1 ∂ t ( h 1 ¯ u 1 ) + ∇ x · ( h 1 ¯ u 1 ⊗ ¯ u 1 + r κ ζ (¯ u 1 − ε ˜ u 2 ) , 1 � ∂ t h 2 + ε ∇ x · ( h 2 ˜ u 2 ) = 0 , u 2 ) = − h 2 κ B (˜ Fr 2 ∇ x ( h 2 + rh 1 + B ) + κ ζ (¯ u 1 − ε ˜ u 2 ) , with κ B (˜ u 2 ) = ˜ κ B ˜ u 2 − ∇ x · ( µ 2 h 2 D x ˜ u 2 ) is derived from the bilayer Navier-Stokes system with the following modeling errors: |H 2 − h 2 | = O ( ε 2 ) , |H 1 − h 1 | = O ( ε ) , u 2 | = O ( ε 2 ) . |U 1 − ¯ u 1 | = O ( ε ) , |U 2 − ε ˜ L.Boittin ( ANGE, Inria ) 30.05.17 9 / 21

Model derivation Proof Momentum eq. and boundary conditions: ∂ z ( µ 2 ∂ z U 2 ) = O ( ε 2 ) ,    ⇒ U 2 = ¯ U 2 ( x , t ) + O ( ε 2 ) µ 2 ∂ z U 2 = O ( ε 2 ) , at z = ζ µ 2 ∂ z U 2 = O ( ε 2 ) , at z = B Vertically integrated horizontal momentum equation: U 2 ) + H 2 ∂ t ( H 2 ¯ + ∇ x · ( H 2 ¯ U 2 ⊗ ¯ U 2 ) Fr 2 ∇ x ( rh 1 + H 2 + B ) κ B ¯ = − ˜ U 2 + 1 ε ∇ x · ( µ 2 h 2 D x ¯ U 2 ) − κ ζ ( ¯ U 2 − ¯ u 1 ) + O ( ε ) , ε Main order terms: κ B ¯ U 2 − ∇ x · ( µ 2 h 2 D x ¯ ˜ U 2 ) = O ( ε ) . Imposes that ¯ U 2 = O ( ε ) = ε ˜ U 2 + O ( ε 2 ). L.Boittin ( ANGE, Inria ) 30.05.17 10 / 21

Model derivation Threshold for incipient motion Classical laws for sediment transport used in hydraulic engineering have a threshold for incipient motion (eg. Meyer-Peter and M¨ uller). Critical shear stress: τ c u 1 − h 2 Effective shear stress: τ eff = κ ζ ¯ Fr 2 ∇ x ( rh 1 + h 2 + B ) Velocity equation: κ B (˜ u 2 ) = τ eff u 2 || α + τ c ) τ eff Take κ B ( · ) such that κ B (¯ u 2 ) = (˜ κ || ¯ || τ eff || − ∇ x · ( µ 2 h 2 ∇ x ¯ u 2 ), with � κ if || τ eff || ≥ τ c ˜ ˜ κ = + ∞ if || τ eff || < τ c L.Boittin ( ANGE, Inria ) 30.05.17 11 / 21

Model analysis Model analysis Dissipative energy balance for the bilayer model For smooth enough solutions: u 1 + ε ∂ t ( K 1 + E ) + ∇ x · ( K 1 u 1 + h 1 φ 1 ¯ r h 2 φ 2 ˜ u 2 ) = − ε u 2 ) − κ ζ u 2 | 2 , r ˜ u 2 · κ B (˜ r | ¯ u 1 − ε ˜ with K 1 = 1 u 1 | 2 : kinematic energy of the water 2 h 1 | ¯ � h 1 � h 2 1 + h 2 � � �� E = h 1 2 + h 2 + B 2 + B : potential energy Fr 2 r 1 1 φ 1 = Fr 2 ( h 2 + h 1 + B ) , φ 2 = rFr 2 ( h 2 + rh 1 + B ): potentials Sediment layer only, without forcing Positivity Maximum principle for smooth solutions L.Boittin ( ANGE, Inria ) 30.05.17 12 / 21

Numerical tests A first idea for the numerical scheme Simplified model: Sediment layer alone No viscosity ⇒ nonlinear heat equation No topography ∂ t h 2 − ε κ B Fr 2 ∇ x · ( h 2 2 ∇ x h 2 ) = 0 . Explicit scheme: parabolic CFL condition ∆ t ≤ C (∆ x ) 2 Try an implicit scheme L.Boittin ( ANGE, Inria ) 30.05.17 13 / 21

Numerical tests 3 different schemes Implicit (linearized) finite volume schemes, staggered grid  h n +1 i +1 / 2 u n +1 i − 1 / 2 u n +1 = h n i − ∆ t ∆ x ( h n i +1 / 2 − h n i − 1 / 2 )  i   u n +1 / 2 ν i +1 ( u n +1 i +3 / 2 − u n +1 i ( u n +1 i +1 / 2 − u n +1 (∆ x ) 2 ( h n i +1 / 2 ) − h n i +1 / 2 − i − 1 / 2 )) h n +1 i +1 − h n +1  = − gh n  i ,  i +1 / 2 ∆ x This scheme dissipates the discrete energy. Centered scheme h n i + h n Take h n i +1 / 2 = i +1 2 Upwind with respect to ∇ h Take h n i +1 / 2 = max( h n i , h n i +1 ) Upwind with respect to u � i if u n +1 h n i +1 / 2 ≥ 0 Take h n i +1 / 2 = . A fixed point is needed. i +1 if u n +1 h n i +1 / 2 < 0 L.Boittin ( ANGE, Inria ) 30.05.17 14 / 21

Numerical tests Comparison of the schemes Regular Irregular 1.05 1 0.95 0.9 0.85 h 0.8 0.75 0.7 0.65 0.6 0.55 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 x Initial conditions L.Boittin ( ANGE, Inria ) 30.05.17 15 / 21

Numerical tests Comparison of the schemes Regular Irregular 0.96 1.7 centered centered upwind h upwind h 0.94 1.6 upwind u upwind u 0.92 1.5 0.9 1.4 0.88 1.3 h T /h 0 h T /h 0 0.86 1.2 0.84 1.1 0.82 1 0.8 0.9 0.78 0.8 0.76 0.7 10 -6 10 -5 10 -4 10 -3 10 -2 10 -1 10 -6 10 -5 10 -4 10 -3 10 -2 10 -1 x= t x= t Infinity norm at final time L.Boittin ( ANGE, Inria ) 30.05.17 15 / 21

Numerical tests Comparison of the schemes Regular Irregular 0.996 1 centered centered upwind h upwind h upwind u upwind u 0.994 0.995 0.992 0.99 normalized energy 0.99 normalized energy 0.985 0.988 0.98 0.986 0.975 0.984 0.97 0.982 0.965 0.98 0.96 10 -6 10 -5 10 -4 10 -3 10 -2 10 -1 10 -6 10 -5 10 -4 10 -3 10 -2 10 -1 x= t x= t Energy at final time L.Boittin ( ANGE, Inria ) 30.05.17 15 / 21

Numerical tests Problems with the centered scheme Energy dissipation, but oscillations! 0.85 1 centered centered upwind upwind h 0.8 upwind u 0.995 0.75 0.7 0.99 normalized energy 0.65 h 0.985 0.6 0.55 0.98 0.5 0.975 0.45 0.4 0.97 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.2 0.4 0.6 0.8 1 1.2 x t Solution at final time, starting from Energy dissipation for the three schemes smooth initial condition L.Boittin ( ANGE, Inria ) 30.05.17 16 / 21

Numerical tests Simulation of the forced model No topography: B ( x ) = 0 . Constant free surface: η = rh 1 + h 2 =constant The two upwind schemes behave differently 0.9 0.8 0.7 h 0.6 0.5 0.4 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 x Water velocity u 1 = 10, density ratio r = 0 . 6 L.Boittin ( ANGE, Inria ) 30.05.17 17 / 21