Simple Eulerian Methods for Compressible Fluids in Domains with - PowerPoint PPT Presentation

Simple Eulerian Methods for Compressible Fluids in Domains with Moving Boundaries Alina Chertock North Carolina State University chertock@math.ncsu.edu joint work with Armando Coco, Yuanzhen Cheng, Smadar Karni, Alexander Kurganov and Giovanni Russo supported in part by

Compressible Euler Equations       ρ ρu ρv ρu 2 + p       ρu ρuv       + + = 0 ρv 2 + p       ρv ρuv E u ( E + p ) v ( E + p ) t x y γ − 1 + ρ p 2( u 2 + v 2 ) E = • ρ is the fluid density; • u, v are the velocities; • E is the total energy; • p is the pressure; � γp • c = is the speed of sound. ρ 1

Two Types of Problems • Multicomponent Flow: We assume that γ is a piecewise constant function propagation with the fluid velocity according to γ t + uγ x + vγ y = 0 . • Problems with Changing Geometries: – The right boundary is a moving – Flow around an (oscillating) solid wall with a prescribed solid circle. equation of motion x = x B ( t ) I II Moving Boundary x B 2

Conservative Methods     ρ ρu ρu 2 + p     ρu     w t + F ( w ) x = 0 ⇔ + = 0     E u ( E + p ) ργ uργ t x Godunov-type discretization: � � j − ∆ t w n +1 w n ¯ = ¯ H j +1 / 2 − H j − 1 / 2 , j ∆ x w n w n H j +1 / 2 = H ( . . . , ¯ j , ¯ j +1 , . . . ) , H ( . . . , w , w , . . . ) = F ( w ) . The major difficulty: to ensure that the pressure equilibrium between fluid components is not disturbed by the numerical scheme. 3

What is going wrong? (S. Karni) 4

Multicomponent Flow     ρ ρu ρu 2 + p     ρu w t + f ( w ) x = 0 ⇐ ⇒ + = 0 E u ( E + p ) t x � � E − 1 2 ρu 2 EOS: p = ( γ − 1) − γp ∞ w = ( ρ, ρu, E ) T We assume that γ is a piecewise constant function propagation with the fluid velocity according to γ t + uγ x = 0 . 5

Multi-Fluid Example pressure CONSERVATIVE METHOD EXACT SOLUTION 1.12 1.1 1.08 1.06 1.04 1.02 1 0.98 0 0.2 0.4 0.6 0.8 1 � (1 . 0 , 1 . 0 , 1 . 0 , 1 . 6) T , if x < 0 . 25 , ( ρ, u, p, γ ) T = (0 . 1 , 1 . 0 , 1 . 0 , 1 . 4) T , if x > 0 . 25 . Different EOS... 6

What is going wrong? t n + 1 J-1 J J+1 t n x J � 1/2 x J + 1/2 material interface • Fluxes are computed using the information in the mixed cell. • No valid EOS in mixed cells. 7

Some Models/Schemes • Conserve not only total mass and momentum but also total energy. – Volume-of-Fluid Method (Noh & Woodward; Collela, Glaz & Ferguson). – Internal Energy Correction Algorithm (Jenny, Mueller & Thomann). – A Fluid-Mixture Type Algorithm (Shyue). • Conserve only total mass and momentum but not total energy. – Pressure Evolution Method (Karni). – Ghost-Fluid Method for the Poor (Abgrall & Karni). • Nonconservative – Ghost-Fluid Method (Fedkiw, Aslam, Merriman & Osher). • – Conservative locally moving mesh method for multifluid flows (A.C. & A. Kurganov) – A second-order finite-difference method for compressible fluids in domains with moving boundaries , (A. C., A. Coco, A. Kurganov and G. Russo) – Interface tracking method for compressible multifluids , (A.C., S. Karni & A. Kurganov) 8

Semi-Discrete Central-Upwind Scheme – 1-D Case     ρ ρu ρu 2 + p     w t + f ( w ) x = 0 ⇐ ⇒ ρu + = 0 E u ( E + p ) t x w := ( ρ, ρu, E ) T Computational cells: C j := [ x j − 1 2 , x j + 1 2 ] � w j ( t ) := 1 The cell averages: w ( x, t ) dx ∆ x C j Time evolution: 2 ( t ) − H j − 1 2 ( t ) H j + 1 d dt w j ( t ) = − ∆ x � � : numerical fluxes H j + 1 2 9

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

� � � � w ± { w j ( t ) } → � w ( · , t ) → 2 ( t ) → 2 ( t ) → { w j ( t + ∆ t ) } H j + 1 j + 1 w ± 2 ( t ) are the point values of j + 1 w ( x, t ) = w j + ( w x ) j ( x − x j ) , � x ∈ C j at x j + 1 2 : 2 + 0 , t ) = w j +1 − ∆ x w + 2 := � w ( x j + 1 2 ( w x ) j +1 j + 1 2 − 0 , t ) = w j + ∆ x w − 2 := � w ( x j − 1 2 ( w x ) j j + 1 11

� � � � w ± { w j ( t ) } → � w ( · , t ) → j ( t ) → 2 ( t ) → { w j ( t + ∆ t ) } H j + 1 2 ( t ) = H ( w − 2 ( t ) , w + 2 ( t )) H j + 1 j + 1 j + 1 Example — Central-Upwind Flux: [Kurganov, Lin, Noelle, Petrova, Tadmor, et al.; 2000–2007] a + 2 f ( w + a + 2 f ( w − 2 ) − a − 2 a − 2 ) � � j + 1 j + 1 j + 1 j + 1 j + 1 j + 1 w + 2 − w − 2 2 = + H j + 1 j + 1 j + 1 a + 2 − a − a + 2 − a − 2 j + 1 j + 1 j + 1 j + 1 2 2 � � a ± 2 ( t ) = a ± w + 2 , w − are one-sided local speeds, estimated using j + 1 j + 1 j + 1 j + 1 2 2 the largest and the smallest eigenvalues of the Jacobian ∂ f ∂ w . 2 ( t ) − H j − 1 2 ( t ) H j + 1 d dt w j ( t ) = − ∆ x 12

Boundary Treatment in 1-D • Assume that at some time level t the wall is located in cell J • All computational cells are divided into 3 groups: – Internal cells with reliable data – Boundary cells with unreliable data – External cells with no data Moving Boundary x J � 5/2 x J � 3/2 x J � 1/2 x B x J+1/2 x J+3/2 13

Boundary Treatment in 1-D • The cell averages w J ( t ) will not be evolved • To evolve the cell averages w J − 1 ( t ) we need: – the numerical flux H J − 1 2 ( t ) – the point values w + 2 ( t ) J − 1 w J � 1 Moving Boundary + w J � 1/2 x B x J � 5/2 x J � 3/2 x J � 1/2 x J+1/2 x J+3/2 14

Solid Wall Extrapolation The solid wall extrapolation of the solution values in cell ( J − 1) , � � E J − 1 − ρ J − 1 u 2 u J − 1 = ( ρu ) J − 1 J − 1 ρ J − 1 , , p J − 1 = ( γ − 1) ρ J − 1 2 To obtain the ”missing” ghost values: � � � � ρ d 2 x B ( t ) ∂p x = x B ( t ) = γp � � � � x = x B ( t ) = − p x x = x B ( t ) , � � � � dt 2 ∂ρ ρ x = x B ( t ) w J � 1 w gh + w J � 1/2 x B x J � 5/2 x J � 3/2 x J � 1/2 x J+1/2 x J+3/2 15

Solid Wall Extrapolation • Second-Order Approximation: · d 2 x B ( t ) ρ J − 1 + ρ gh = − p gh − p J − 1 dt 2 2 ∆ x p gh − p J − 1 = γ · p gh + p J − 1 ρ gh − ρ J − 1 ρ gh + ρ J − 1 • First-Order Approximation: ρ gh := ρ J − 1 , u gh := 2 ˙ x B ( t ) − u J − 1 , p gh := p J − 1 x B ( t ) : velocity of the wall at time t ˙ 16

Phase Space Interpolation Once w + if u ∗ > 0 then are available, we J − 1 � w ∗ , 2 if u ∗ − c ∗ < 0 compute the slopes ( w x ) J − 1 w + 2 = J − 1 w J − 1 , otherwise using else � w ∗ , 2 , w J − 1 and w + w − if u ∗ + c ∗ > 0 w + J − 3 J − 1 2 = 2 J − 1 w gh , otherwise w � J � 1 ρ w J � 3/2 * u Moving Boundary * p ρ ρ * J−1 gh + w J � 1/2 u u gh J−1 p p gh J−1 x J � 5/2 x J � 3/2 x J � 1/2 x B x J+1/2 x J+3/2 17

Solution of the RP – A Very Simple Structure • If u gh < u J − 1 , the solution consists of three constant states, connected by two shock waves. The intermediate state is equal to: p ∗ = p J − 1   � � � 2 x B ( t ) − u J − 1 ) 2 + ( γ + 1) ρ J − 1 ( ˙ 4 c J − 1  1 +  1 + 4 ( γ + 1) ( ˙ x B ( t ) − u J − 1 ) ( γ − 1) p J − 1 + ( γ + 1) p ∗ u ∗ = ˙ x B ( t ) , ρ ∗ = ρ J − 1 ( γ + 1) p J − 1 + ( γ − 1) p ∗ • If u gh > u J − 1 , the solution consists of three constant states, connected by two rarefaction waves.The intermediate state is given by: � � 2 γ/ ( γ − 1) 1 − ( γ − 1) ( ˙ x B ( t ) − u J − 1 ) p ∗ = p J − 1 2 c J − 1 � p ∗ � 1 /γ u ∗ = ˙ x B ( t ) , ρ ∗ = ρ J − 1 p J − 1 18

Time Evolution 2 ( t ) − H j − 1 2 ( t ) H j + 1 d dt w j ( t ) = − : w j ( t ) → w j ( t + ∆ t ) ∆ x • If x B ( t + ∆ t ) ∈ C J , the evolution step is completed. • If x B ( t + ∆ t ) ∈ C J − 1 , then – Cell J − 1 becomes the new boundary cell – Cell J turns into an external one – The computed values w J − 1 ( t + ∆ t ) will not be used any more J − 1 t W J � 1 t+ � t t x x J � 3/2 x J � 1/2 x J+1/2 19

• If x B ( t + ∆ t ) ∈ C J +1 , then – Cell ( J + 1) becomes a new boundary cell – Cell J turns into an internal one – We use the phase space interpolation between ( ρ J − 1 , u J − 1 , p J − 1 ) and ( ρ gh , u gh , p gh ) to obtain ( ρ ∗ , u ∗ , p ∗ ) and then w J ( t + ∆ t ) := w ∗ ∈ t W J t+ � t t x x J � 1/2 x J+1/2 x J+3/2 20