Numerical Hydrodynamics: A Primer Wilhelm Kley Institut f ur - PowerPoint PPT Presentation

Numerical Hydrodynamics: A Primer Wilhelm Kley Institut f¨ ur Astronomie & Astrophysik & Kepler Center for Astro and Particle Physics T¨ ubingen Bad Honnef 2016

Hydrodynamics: Hydrodynamic Equations The Euler-Equations in conservative form read ∂ρ ∂ t + ∇· ( ρ� u ) = 0 (1) ∂ ( ρ� u ) −∇ p + ρ� + ∇· ( ρ� u ⊗ � u ) = k (2) ∂ t ∂ ( ρǫ ) + ∇· ( ρǫ� − p ∇· � u ) = u (3) ∂ t u : Velocity, � � k : external forces, ǫ specific internal energy The equations describe the conservation of mass, momentum and energy. For completion we need an equation of state (eos): p = ( γ − 1 ) ρǫ (4) Using this and eq. (3), we can rewrite the energy equation as an equation for the pressure ∂ p ∂ t + ∇· ( p � u ) = − ( γ − 1 ) p ∇ · � u (5) W. Kley Numerical Hydrodynamics: A Primer Bad Honnef 2016 2

Hydrodynamics: Reformulating Expanding the divergences on the left side and use for the momentum and energy equation the continuity equation ∂ρ ∂ t + ( � − ρ ∇ · � u · ∇ ) ρ = u (6) ∂� u − 1 ρ ∇ p + � ∂ t + ( � u · ∇ ) � u = k (7) ∂ p ∂ t + ( � − γ p ∇· � u · ∇ ) p = u (8) Since all quantities depend on space ( � r ) and time ( t ), for example ρ ( � r , t ) , we can use for the left side the total time derivative (Lagrange-Formulation). For example, for the density one obtains D ρ Dt = ∂ρ ∂ t + ( � u · ∇ ) ρ = − ρ ∇ · � u . (9) The Operator Dt = ∂ D ∂ t + � u · ∇ (10) is called material derivative (equivalent to the total time derivative, d / dt ). W. Kley Numerical Hydrodynamics: A Primer Bad Honnef 2016 3

Hydrodynamics: Lagrange-Formulation Use now the material derivative D ρ − ρ ∇ · � = u (11) Dt D � u − 1 ρ ∇ p + � = k (12) Dt Dp − γ p ∇· � = u (13) Dt These equations describe the change of the quantities in the comoving frame = Lagrange-Formulation. For the Euler-Formulation, one analysed the changes at a specific, fixed point in space ! The Lagrange-Formulation can be used conveniently for 1D-problems, for example the radial stellar oscillations, using comoving mass-shells. For the Euler-Formulation a fixed grid is used. W. Kley Numerical Hydrodynamics: A Primer Bad Honnef 2016 4

Numerical Hydrodynamics: The problem Consider the evolution of the full time-dependent hydrodynamic equations. The non-linear partial differential equations of hydrodynamics will be solved numerically Continuum ⇒ Discretisation W. Kley Numerical Hydrodynamics: A Primer Bad Honnef 2016 5

Numerical Hydrodynamics: Method of solution Particle methods (Lagrange) Grid-based methods (Euler) moving Grid/Particle fixed Grid - flow moved grid - matter flows through grid ρ d � u dt = −∇ p � ∂� � u ∂ t + � u · ∇ � ρ u = −∇ p Well known method: Smoothed Particle Hydrodynamics , SPH Methods: • finite differences non-conservative • Control Volume conservative • Riemann-solver wave properties ’smeared out particles’ • Problem: Discontinuities good for free boundaries, self-gravity W. Kley Numerical Hydrodynamics: A Primer Bad Honnef 2016 6

Numerical Hydrodynamics: consider: 1D Euler equations describe conservation of mass, momentum and energy ∂ρ ∂ t + ∂ρ u = 0 (14) ∂ x ∂ρ u + ∂ρ uu − ∂ p = (15) ∂ t ∂ x ∂ x ∂ρǫ ∂ t + ∂ρǫ u − p ∂ u = (16) ∂ x ∂ x ρ : density u : velocity p : pressure ǫ : internal specific energy (Energy/Mass) with the equation of state p = ( γ − 1 ) ρǫ (17) γ : adiabatic exponent partial differential equation in space and time → need discretisation in space and time. W. Kley Numerical Hydrodynamics: A Primer Bad Honnef 2016 7

Numerical Hydrodynamics: Discretisation ψ Consider funktion: ψ ( x , t ) discretisation in space cover with a grid ∆ x = x max − x min N X j j−1 j+1 ψ n j cell average of ψ ( x , t ) at the grdipoint x j at time t n � ( j + 1 / 2 )∆ x 1 ψ n j = ψ ( x j , t n ) ≈ ψ ( x , n ∆ t ) dx ∆ x ( j − 1 / 2 )∆ x ψ n j is piecewise constant. j spatial index, n time step. W. Kley Numerical Hydrodynamics: A Primer Bad Honnef 2016 8

Numerical Hydrodynamics: Integration in time Consider general equation ∂ψ ∂ t = L ( ψ ( x , t )) (18) with a (spatial) differential operator L . typical Discretisation (1. order in time), at time: t = t n = n ∆ t = ψ n + 1 − ψ n ∂ψ ∂ t ≈ ψ ( t + ∆ t ) − ψ ( t ) = L ( ψ n ) (19) ∆ t ∆ t now at a special location, the grid point x j (with moving terms) ψ n + 1 = ψ n j + ∆ tL ( ψ n k ) (20) j L ( ψ n k ) : discretised differential operator L (here explicit) - k in L ( ψ k ): set of spatial indices: - typical for 2. order: k ∈ { j − 2 , j − 1 , j , j + 1 , j + 2 } (need information from left and right, 5 point ’Stencil’) W. Kley Numerical Hydrodynamics: A Primer Bad Honnef 2016 9

Numerical Hydrodynamics: Operator-Splitting ∂� A ∂ t = L 1 ( � A ) + L 2 ( � A ) (21) L i ( � A ) , i = 1 , 2 individul (Differential-)operators applied to the quantities � A = ( ρ, u , ǫ ) . Here, for 1D ideal hydrodynamics L 1 : Advection L 2 : pressure, or external forces To solve the full equation the solution is split in several substeps A n + ∆ tL 1 ( � � � A 1 A n ) = A n + 1 = � A 1 + ∆ tL 2 ( � � � A 2 A 1 ) = (22) L i is the differential operator to L i . W. Kley Numerical Hydrodynamics: A Primer Bad Honnef 2016 10

Numerical Hydrodynamics: advection-step ∂ρ − ∂ρ u = ∂ t ∂ x ∂ ( ρ u ) − ∂ ( ρ uu ) = ∂ t ∂ x ∂ ( ρǫ ) − ∂ ( ρǫ u ) = ∂ t ∂ x In explicit conservation form ∂ t + ∂� ∂� f ( � u u ) = 0 (23) ∂ x u = ( u 1 , u 2 , u 3 ) and � for � f = ( f 1 , f 2 , f 3 ) we have: u = ( ρ, ρ u , ρǫ ) and � � f = ( ρ u , ρ uu , ρǫ u ) . ρ n → ρ 1 = ρ n + 1 , u n → u 1 , ǫ n → ǫ 1 This step yields: W. Kley Numerical Hydrodynamics: A Primer Bad Honnef 2016 11

Numerical Hydrodynamics: Force terms Momentum equation ∂ u ∂ t = − 1 ∂ p (24) ρ ∂ x � � p j − p j − 1 1 u n + 1 = u j − ∆ t for j = 2 , N (25) j ρ n + 1 ∆ x ¯ j energy equation ∂ǫ ∂ t = − p ∂ u (26) ρ ∂ x � � p j u j + 1 − u j ǫ n + 1 = ǫ j − ∆ t for j = 1 , N (27) j ρ n + 1 ∆ x j on the right hand side we use the actual values for u , ǫ and p , i.e. here u 1 , p 1 , ǫ 1 . u 1 → u n + 1 , ǫ 1 → ǫ n + 1 This step yields: W. Kley Numerical Hydrodynamics: A Primer Bad Honnef 2016 12

Numerical Hydrodynamics: Model equation for advection The continuity equation was ∂ρ ∂ t + ∂ρ u ∂ x = 0 (28) Here, F m = ρ u is the mass flow Using the notation ρ → ψ and u → a = const . we obtain the Linear Advection Equation ∂ψ ∂ t + a ∂ψ ∂ x = 0 . (29) with a constant velocity a , the solution is a wave traveling to the right using ψ ( x , t = 0 ) = f ( x ) we get ψ ( x , t ) = f ( x − at ) Here f ( x ) is the initial condition at time t = 0, that is shifted by the advection with a constant velocity a to the right. The numerics should maintain this property as accurately as possible. W. Kley Numerical Hydrodynamics: A Primer Bad Honnef 2016 13

Numerical Hydrodynamics: Linear Advection FTCS: Forward Time Centered Space algorithm ∂ψ ∂ t + a ∂ψ ∂ x = 0 (30) ψ ψ ψ j j−1 j+1 Specify the grid : x x x j−1 j j+1 and write � ψ n + 1 n − ψ n ∂ψ � j j = (31) � ∂ t ∆ t � j � n ψ n j + 1 − ψ n ∂ψ � j − 1 = (32) � ∂ x 2 ∆ x � j it follows j − a ∆ t � � ψ n + 1 = ψ n ψ n j + 1 − ψ n (33) j − 1 j 2 ∆ x The method looks well motivated: but it is unstable for all ∆ t ! W. Kley Numerical Hydrodynamics: A Primer Bad Honnef 2016 14

Numerical Hydrodynamics: Upwind-Method I ∂ψ ∂ t + ∂ a ψ = 0 (34) ∂ x or α ∆ ψ ∂ψ ∂ t + a ∂ψ a t ∂ x = 0 (35) a : constant (velocity) > 0 ψ ( x , t ) arbitrary transport quantity change of ψ in grid cell j in out ψ n + 1 ∆ x = ψ n j ∆ x +∆ t ( F in − F out ) (36) j X The flux F in is for constant ψ j j j−1 j+1 purple regions will be trans- a ψ n F in = (37) j − 1 ported into the next neighbour a ψ n F out = (38) j cell Upwind-Method Information comes from upstream W. Kley Numerical Hydrodynamics: A Primer Bad Honnef 2016 15

Numerical Hydrodynamics: Upwind-Method II α Extension for non-constant states ψ a ∆ t � � x j − 1 / 2 − a ∆ t ψ (x) F in = a ψ I (39) I 2 ψ I ( x ) interpolation polynomial Here linear interpolation (straight line) in out This yields � + 1 � X ψ n F in = a 2 ( 1 − σ )∆ ψ j − 1 (40) j j−1 j+1 j − 1 � �� � j + x − x j 1 st Order ψ I ( x ) = ψ n ∆ x ∆ ψ j (41) � �� � 2 nd Order ∆ ψ j undivided differences with σ = a ∆ t / ∆ x 2nd order upwind � ∆ ψ j ≈ ∂ψ � ∆ x � ψ I ( x ) is evaluated in the middle ∂ x � x j of the purple area. W. Kley Numerical Hydrodynamics: A Primer Bad Honnef 2016 16