Central-Upwind Scheme for the Shallow Water System with Horizontal Temperature Gradients Alina Chertock North Carolina State University chertock@math.ncsu.edu joint work with Alexander Kurganov, Tulane University Yu Liu, Tulane University
Saint-Venant System of Shallow Water h t + ( hu ) x = 0 hu 2 + g � 2 h 2 � ( hu ) t + x = − ghB x z w=B+h h(x,t) B(x) 1
The 1-D Ripa System h t + ( hu ) x = 0 hu 2 + g � � 2 h 2 θ ( hu ) t + x = − ghθB x ( hθ ) t + ( uhθ ) x = 0 • h : water height • u : fluid velocity • θ : potential temperature. Specifically, θ is the reduced gravity g ∆Θ / Θ ref computed as the potential temperature difference ∆Θ from some reference value Θ ref . • B : bottom topography • g : gravitational constant If θ ≡ const , then the Ripa system reduces to the Saint-Venant system of shallow water equations 2
The Ripa System • Introduced in [Ripa [(1993,1995), Dellar (2003)] for modeling ocean currents. • The derivation of the system is based on considering multilayered ocean models, and vertically integrating the density, horizontal pressure gradient and velocity fields in each layer. • The model incorporates the horizontal temperature gradients, which results in the variations in the fluid density within each layer. 3
The Ripa System h t + ( hu ) x = 0 hu 2 + g � � 2 h 2 θ x = − ghθB x ( hu ) t + ( hθ ) t + ( uhθ ) x = 0 admits the energy (entropy) inequality � hu 2 + gh 2 θ � hu 2 � � �� + gh 2 θ + ghθB ≤ 0 2 + ghθB + u 2 2 t x • The eigenvalues are u ± √ ghθ, u, 0 • There is a nonlinear resonance when u ± c = 0 (wave speeds from different families of waves coincide) • There are no Riemann invariants for the Ripa system and therefore it is very hard to design upwind schemes since they are based on (approximate) Riemann problem solvers 4
Balance Law U t + f ( U ) x = S ( U ) f := ( hu, hu 2 + g 2 h 2 θ, uhθ ) T , U := ( h, hu, hθ ) T , S := (0 , − ghB x , 0) T Semi-discrete central-upwind scheme: 2 − H j − 1 H j + 1 d 2 dt U j = − + S j , ∆ x j ≈ 1 � n U ( x, t n ) dx, C j := ( x j − 1 2 , x j + 1 2 ) U ∆ x C j Numerical Challenges: • well-ballanced • positivity preserving 5
Steady States ( hu ) x = 0 ( hu ) x = 0 ⇐ ⇒ � u 2 hu 2 + g � = g � � 2 h 2 θ x = − ghθB x 2 + gθ ( h + B ) 2 hθ x x The system cannot be integrated, but admits several particular steady-state solutions, two of them are the following “lake at rest” ones: 1. θ ≡ constant , w := h + B ≡ constant , u ≡ 0 corresponds to flat water surface under the constant temperature p := g 2 h 2 θ ≡ constant , 2. B ≡ constant , u ≡ 0 corresponds to the contact wave across which h and θ jump while u and p remain constant Goal: to derive a well-balanced scheme which preserves both steady states! 6
Well-Balanced Scheme 1. θ ≡ constant , w := h + B ≡ constant , u ≡ 0 Well-balanced scheme should exactly balance the flux and source terms so that the steady state is preserved – the same approach as in the case of the central-upwind scheme for the Saint-Venant system (Kurganov & Petrova, 2007) p := g 2 h 2 θ ≡ constant , 2. B ≡ constant , u ≡ 0 In a well-balanced scheme, the pressure should remain oscillation-free across the temperature jump and thus the steady state will be exactly preserved – the same approach as in the case of the interface tracking method for compressible multifluids (Chertock, Karni & Kurganov, 2008) 7
Small Perturbation of Steady-State – Numerical Example 0 . 85(cos(10 π ( x + 0 . 9)) + 1) , − 1 . 0 ≤ x ≤ − 0 . 8 , B ( x ) = 1 . 25(cos(10 π ( x − 0 . 4)) + 1) , 0 . 3 ≤ x ≤ 0 . 5 , 0 , otherwise . It is easy to see that � (6 , 0 , 4) T , x < 0 ( w s , u s , θ s ) T ( x ) = (4 , 0 , 9) T , x > 0 is a piecewise constant steady-state solution, which is in fact a combination of two “lake at rest” steady states of type I connected through the temperature jump, which corresponds to a steady state of type II. 8
2.5 6 2 5.5 1.5 w B 5 1 4.5 0.5 4 0 −2 −1 0 1 2 −2 −1 0 1 2 t=0.1 t=0.4 6 6 5.5 5.5 w 5 w 5 4.5 4.5 4 4 −2 −1 0 1 2 −2 −1 0 1 2 9
Pressure Oscillations – Numerical Example • The initial condition is √ � (2 2 , 4 , 1) , x < 0 ( h ( x, 0) , u ( x, 0) , θ ( x, 0)) = (1 , 4 , 8) , x > 0 • Notice that p = 4 g for all x , thus initially there is no pressure jump • g = 1 • The bottom topography is flat 10
t=0 t=0 5 5 4.5 4.5 4 4 p p 3.5 3.5 3 3 −1000 −500 0 500 1000 −1000 −500 0 500 1000 t=10 t=10 4.4 4.4 4.3 4.3 p−zoom p−zoom 4.2 4.2 4.1 4.1 4 4 0 20 40 60 80 0 20 40 60 80 t=50 t=50 4.3 4.3 4.25 4.25 4.2 4.2 p−zoom p−zoom 4.15 4.15 4.1 4.1 4.05 4.05 4 4 3.95 3.95 100 150 200 250 300 350 400 100 150 200 250 300 350 400 11
Switching to Equilibrium Variable h t + ( hu ) x = 0 hu 2 + g � � 2 h 2 θ x = − ghB x ( hu ) t + ( hθ ) t + ( uhθ ) x = 0 � ( h, hu, hθ ) → ( w := h + B, hu, hθ ) w t + ( hu ) x = 0 � ( hu ) 2 w − B + g � 2 θ ( w − B ) 2 ( hu ) t + = − gθ ( w − B ) B x x ( hθ ) t + ( huθ ) x = 0 12
Semi-Discrete Central-Upwind Scheme Central-upwind schemes were developed for multidimensional hyperbolic systems of conservation laws in 2000–2007 by Kurganov, Lin, Noelle, Petrova, Tadmor, ... Central-upwind schemes are Godunov-type finite-volume projection- evolution methods: • at each time level a solution is globally approximated by a piecewise polynomial function, • which is then evolved to the new time level using the integral form of the system of hyperolic conservation/balance laws. 13
Semi-Discrete Central-Upwind Scheme 2 − H j − 1 H j + 1 d q := ( w, hu, hθ ) T 2 dt q j = − + S j , ∆ x a + q + a + q − − a − 2 a − � � � � 2 , B j + 1 2 , B j + 1 2 f 2 f j + 1 j + 1 j + 1 j + 1 j + 1 j + 1 � � 2 2 q + 2 − q − 2 2 = + H j + 1 j + 1 j + 1 a + a + 2 − a − 2 − a − 2 j + 1 j + 1 j + 1 j + 1 2 2 • q ± 2 : right/left point values at x j + 1 2 of a piecewise polynomial j + 1 reconstruction • a ± 2 : local right-/left-sided speeds j + 1 • B j + 1 2 = B ( x j + 1 2 ) 14
Reconstruction of Equilibrium Variables • To preserve the first steady state, we reconstruct the equilibrium variables ( θ, hu, w ) and obtain their point values at x j + 1 2 : 2 = θ j + ∆ x 2 = θ j +1 − ∆ x θ + θ − 2 ( θ x ) j , 2 ( θ x ) j j + 1 j + 1 2 = ( hu ) j + ∆ x 2 = ( hu ) j +1 − ∆ x ( hu ) − ( hu ) + 2 (( hu ) x ) j , 2 (( hu ) x ) j j + 1 j + 1 2 = w j + ∆ x 2 = w j +1 − ∆ x w − w + 2 ( w x ) j , 2 ( w x ) j j + 1 j + 1 • The point values of h , u and hθ are then computed as follows: ( hu ) ± j + 1 h ± 2 = w ± u ± ( hθ ) ± 2 = h ± 2 θ ± 2 2 − B j + 1 2 , 2 = , j + 1 j + 1 j + 1 j + 1 j + 1 j + 1 h ± 2 j + 1 2 15
Preservation of Positivity h ± 2 = w ± 2 − B j + 1 j + 1 j + 1 2 Step 1: Piecewise linear reconstruction of the bottom B j+1/2 B (x) B B j+1 j x j−1/2 x j+1/2 x j+3/2 16
Step 2: Positivity preserving reconstruction of w h ± 2 = w ± 2 − B j + 1 j + 1 j + 1 2 B j−1/2 w j − + w w j−1/2 j+1/2 B j B j+1/2 x x x j−1/2 j+1/2 j 17
+ w B j−1/2 B j−1/2 = j−1/2 w j w j − + w w j−1/2 j+1/2 − w j+1/2 B j B j B j+1/2 B j+1/2 x x x x x x j−1/2 j+1/2 j−1/2 j+1/2 j j if w − then take w − 2 , w + 2 = 2 w j − B j + 1 2 < B j + 1 2 = B j + 1 j + 1 j + 1 j − 1 2 2 if w + 2 , w + then take w − 2 < B j − 1 2 = 2 w j − B j − 1 2 = B j − 1 j − 1 j + 1 j − 1 2 2 18
We have proved that if an SSP ODE solver is used, then ¯ h n +1 = α − 2 h − 2 + α + 2 h + 2 + α − 2 h − 2 + α + 2 h + j − 1 j − 1 j − 1 j − 1 j + 1 j + 1 j + 1 j + 1 j 2 and n +1 = β − 2 + β + 2 h + 2 θ + 2 + β + 2 h + 2 θ + 2 h − 2 θ − 2 + β − 2 h − 2 θ − hθ j j − 1 j − 1 j − 1 j − 1 j − 1 j − 1 j + 1 j + 1 j + 1 j + 1 j + 1 j + 1 2 where the coefficients α ± 2 > 0 and β ± 2 > 0 provided an appropriate CFL j ± 1 j ± 1 condition is satisfied. This guarantees positivity of both h and θ = hθ h . 19
Approximation of the Source Term 2 ≡ ˆ Substitute the “lake at rest” values θ ± θ, ( hu ) ± 2 ≡ 0 , w ± 2 ≡ ˆ w into j + 1 j + 1 j + 1 the scheme ⇒ the numerical fluxes H j + 1 2 reduce to: 0 , g 2 , H (3) � T � T � � H (1) 2 , H (2) ˆ 2 ) 2 , 0 = θ ( ˆ w − B j + 1 j + 1 j + 1 2 Thus H (2) 2 − H (2) = − g ˆ θ j + 1 j − 1 � 2 ) 2 − ( ˆ 2 ) 2 � 2 − w − B j + 1 w − B j − 1 ( ˆ ∆ x 2∆ x � B j + 1 = g ˆ 2 − B j − 1 θ � 2 w − B j + 1 w − B j − 1 ˆ 2 + ˆ 2 ∆ x 2 The well-balanced quadrature: xj +1 2 = − g � (2) θ ( w − B ) B x dx S j ∆ x xj − 1 2 �� B j +1 2 − B j − 1 ≈ − g � � � � + θ + w + θ − w − 2 − B j +1 − B j − 1 j +1 j +1 j − 1 j − 1 2 ∆ x 2 2 2 2 2 2 20
Recommend
More recommend
