Structure Preserving Numerical Methods for Hyperbolic Systems of - PowerPoint PPT Presentation

Structure Preserving Numerical Methods for Hyperbolic Systems of Conservation and Balance Laws Alina Chertock North Carolina State University chertock@math.ncsu.edu joint work with S. Cui, M. Herty, A. Kurganov, X. Liu, S.N. ¨ Ozcan and E. Tadmor

Systems of Balance Laws U t + f ( U ) x + g ( U ) y = 1 U t + f ( U ) x + g ( U ) y = S ( U ) ε S ( U ) Examples: Examples: • low Mach number compressible • Gas dynamics with pipe-wall flows friction • low Froude number shallow • Euler equations with water flows gravity/friction • diffusive relaxation in kinetic • shallow water equations with models Coriolis forces Applications: Applications: • various two-phase flows such as • astrophysical and atmospheric bubbles in water phenomena in many fields including supernova explosions • unmostly incompressible flows with regions of high • (solar) climate modeling and compressibility such as weather forecasting underwater explosions • atmospheric flows 1

Systems of Balance Laws U t + f ( U ) x + g ( U ) y = 1 U t + f ( U ) x + g ( U ) y = S ( U ) or ε S ( U ) • Challenges: certain structural properties of these hyperbolic problems (conservation or balance law, equilibrium state, positivity, assymptotic regimes, etc.) are essential in many applications; • Goal: to design numerical methods that are not only consistent with the given PDEs, but – preserve the structural properties at the discrete level – well-balanced numerical methods – remain accurate and robust in certain asymptotic regimes of physical interest – asymptotic preserving numerical methods [P. LeFloch; 2014] 2

Well-Balanced (WB) Methods U t + f ( U ) x + g ( U ) y = S ( U ) • In many physical applications, solutions of the system are small perturbations of the steady states; • These perturbations may be smaller than the size of the truncation error on a coarse grid; • To overcome this difficulty, one can use very fine grid, but in many physically relevant situations, this may be unaffordable; Goal: • to design a well-balanced numerical method, that is, the method which is capable of exactly preserving some steady state solutions; • perturbations of these solutions will be resolved on a coarse grid in a non-oscillatory way. 3

Asymptotic Preserving (AP) Methods U t + f ( U ) x + g ( U ) y = 1 ε S ( U ) • Solutions of many hyperbolic systemes reveal a multiscale character and thus their numerical resolution presence some major difficulties; • Such problems are typically characterized by the occurence of a small parameter by 0 < ε ≪ 1 ; • The solutions show a nonuniform behavior as ε → 0 ; • the type of the limiting solution is different in nature from that of the solutions for finite values of ε > 0 . Goal: • asymptotic passage from one model to another should be preserved at the discrete level; • for a fixed mesh size and time step, AP method should automatically transform into a stable discretization of the limitting model as ε → 0 . 4

Finite-Volume Methods – 1-D � � = 1 U t + f ( U ) x = S ε S � k ≈ 1 n U ( y, t n ) dy : cell averages over C j := ( x j − 1 • U 2 , x j + 1 2 ) ∆ y C k • Semi-discrete FV method: 2 ( t ) − F j − 1 2 ( t ) F j + 1 d dt U j ( t ) = − + S j ∆ x 2 ( t ) : numerical fluxes F j + 1 S j : quadrature approximating the corresponding source terms • Central-Upwind (CU) Scheme: [Kurganov, Lin, Noelle, Petrova, Tadmor, et al.; 2000–2007] 5

� � � � { U j ( t ) } → � U E , W U ( · , t ) → ( t ) → 2 ( t ) → { U j ( t + ∆ t ) } F j + 1 j (Discontinuous) piecewise-linear reconstruction: � U ( y, t ) := U j ( t ) + ( U x ) j ( x − x j ) , x ∈ C j It is conservative, second-order accurate, and non-oscillatory provided the slopes, { ( U y ) k } , are computed by a nonlinear limiter Example — Generalized Minmod Limiter � � θ U j − U j − 1 , U j +1 − U j − 1 , θ U j +1 − U j ( U y ) j = minmod ∆ x 2∆ x ∆ x  where min j { z j } , if z j > 0 ∀ j,   minmod( z 1 , z 2 , ... ) := max j { z j } , if z j < 0 ∀ j,   0 , otherwise , and θ ∈ [1 , 2] is a constant 6

� � � � { U j ( t ) } → � U E , W U ( · , t ) → ( t ) → 2 ( t ) → { U j ( t + ∆ t ) } F j + 1 j U E j and U W are the point values at x j + 1 2 and x j − 1 2 : j � U ( y, t ) = U j + ( U x ) j ( x − x j ) , x ∈ C j j := U j + ∆ x k+1/2 U E 2 ( U x ) j := U j − ∆ x U W 2 ( U x ) j k j k−1/2 j j−1/2 j+1/2 7

� � � � { U j ( t ) } → � U E , W U ( · , t ) → ( t ) → 2 ( t ) → { U j ( t + ∆ t ) } F j + 1 j F j + 1 2 − F j − 1 d 2 dt U j = − + S j ∆ x where a + j ) − a − 2 f ( U E 2 f ( U W j +1 ) � � j + 1 j + 1 U W j +1 − U E 2 = + α j + 1 F j + 1 j a + 2 − a − 2 j + 1 j + 1 2 a + 2 a − j + 1 j + 1 2 α j + 1 2 = a + 2 − a − j + 1 j + 1 2 � � � � a + a − λ ( U E j ) , λ ( U W λ ( U E j ) , λ ( U W 2 = max j +1 ) , 0 , 2 = min j +1 ) , 0 j + 1 j + 1 2-D extension is dimension-by-dimension 8

Non Well-Balanced Property – Example � ρ t + q x = 0 , q t + f 2 ( ρ, q ) x = − s ( ρ, q ) For steady-state solution: q = Const and ρ = ρ ( x ) Implementing the CU scheme results in  ❍❍❍❍❍❍❍❍❍❍❍❍❍❍❍❍❍ ✟ ✟✟✟✟✟✟✟✟✟✟✟✟✟✟✟✟✟ a + j − a − 2 q E 2 q W dρ j dt = − 1 j +1 j + 1 j + 1  2 ( ρ W j +1 − ρ E + α j + 1 j )  a + 2 − a − ∆ x j + 1 j + 1 2 ❍  ✟✟✟✟✟✟✟✟✟✟✟✟✟✟✟✟✟✟ ❍❍❍❍❍❍❍❍❍❍❍❍❍❍❍❍❍❍ a + j − 1 − a − 2 q E 2 q W j j − 1 j − 1  2 ( ρ W j − ρ E − + α j − 1 j − 1 )  � = 0 a + 2 − a − j − 1 j − 1 2 • The steady state would not be preserved at the discrete level; • This would also true for the first-order version of the scheme; • For smooth solutions, the balance error is expected to be of order (∆ x ) 2 , but a coarse grid solution may contain large spurious waves. 9

Well-Balanced Methods “Balance is not something you find, it’s something you create” 10

1-D 2 × 2 Systems of Balance Laws � ρ t + f 1 ( ρ, q ) x = 0 , q t + f 2 ( ρ, q ) x = − s ( ρ, q ) , Steady state solution : f 1 ( ρ, q ) x ≡ 0 , f 2 ( ρ, q ) x + s ( ρ, q ) ≡ 0 or K := f 1 ( ρ, q ) ≡ Const , � x ∀ x, t L := f 2 ( ρ, q ) + s ( ρ, q ) dξ ≡ Const Numerical Challenges : to exactly balance the flux and source terms, i.e., to exactly preserve the steady states. How to design a well-balanced scheme? 11

Well-Balanced Scheme � ρ t + f 1 ( ρ, q ) x = 0 , q t + f 2 ( ρ, q ) x = − s ( ρ, q ) • Incorporate the source term into the flux: � � x ρ t + f 1 ( ρ, q ) x = 0 , R := s ( ρ, q ) dξ q t + ( f 2 ( ρ, q ) x + R ) x = 0 , • Rewrite � ρ t + K x = 0 , q t + L x = 0 where K := f 1 ( ρ, q ) , L := f 2 ( ρ, q ) x + R • Define conservative variables U = ( ρ, q ) T equilibrium variables W := ( K, L ) T 12

Well-Balanced Scheme U t + f ( U ) x = 0 � � � � ρ K U = , f ( U ) = W := q L Semi-discrete FV method: 2 ( t ) − F j − 1 2 ( t ) F j + 1 d dt U j ( t ) = − ∆ x Two major modifications : • Well-balanced reconstruction – performed on the equilibrium rather than conservative variables : � � � � � � W E , W U E , W { U j ( t ) } → � U ( · , t ) → ( t ) → ( t ) → 2 ( t ) → { U j ( t +∆ t ) } F j + 1 j j • Well-balanced evolution 13

Well-Balanced Reconstruction Given : U j ( t ) = ( ρ j , q j ) T – cell averages ) T – point values, where Need : W E , W = ( K E , W , L E , W j j j � x K := f 1 ( ρ, q ) , L := f 2 ( ρ, q ) x + R, R := s ( ρ, q ) dξ x j � • Compute R j = s ( ρ, q ) dξ by the midpoint quadrature rule and using the following recursive relation: R j = 1 R 1 / 2 ≡ 0 , 2( R j − 1 2 + R j + 1 2 ) , R j + 1 2 = R ( x j + 1 2 ) = R j − 1 2 + ∆ x s ( x j ,ρ j ,q j ) • Compute the point values of K and L at x j from the cell averages, ρ j and q j : K j = f 1 ( ρ j ,q j ) , L j = f 2 ( ρ j ,q j ) + R j 14

Well-Balanced Reconstruction • Apply the minmod reconstruction procedure to { K j , L j } and obtain the point values at the cell interfaces: j = K j + ∆ x j = L j + ∆ x K E L E 2 ( K x ) j , 2 ( L x ) j , = K j − ∆ x j = L j − ∆ x K W L W 2 ( K x ) j , 2 ( L x ) j j • Finally, equipped with the values of K E , W , L E , W and R j ± 1 2 , solve j j K E j = f 1 ( ρ E j , q E L E j = f 2 ( ρ E j , q E j ) , j ) + R j + 1 2 , K W = f 1 ( ρ W j , q W L W j = f 2 ( ρ W j , q W j ) , j ) + R j − 1 j 2 for U E , W = ( ρ E , W , q E , W ) T . j j j 15

Well-Balanced Evolution 2 − F j − 1 F j + 1 d 2 dt U j = − ∆ x 1 where 0.8 0.6 H 0.4 a + j − a − 2 K E 2 K W 0.2 j +1 j + 1 j + 1 F (1) 2 = j + 1 a + 2 − a − 0 0 0.01 0.02 0.03 0.04 j + 1 j + 1 ψ 2 � | K j +1 − K j | � · | Ω | 2 ( ρ W j +1 − ρ E + α j + 1 j ) H K j , K j +1 } , ∆ x max j a + j − a − 2 L E 2 L W j +1 j + 1 j + 1 F (2) 2 = j + 1 a + 2 − a − j + 1 j + 1 2 � | L j +1 − L j | � | Ω | 2 ( q W j +1 − q E + α j + 1 j ) H · , ∆ x max j { L j , L j +1 } 16

Proof of the Well-Balanced Property Theorem . The central-upwind semi-discrete schemes coupled with the well-balanced reconstruction and evolution is well-balanced in the sense that it preserves the corresponding steady states exactly. 17