ME433 (2020SS) Project I – 1D Unsteady Compressible Navier-Stokes Due 5pm on April 27, 2020 1 Shock Tube or Riemann Problem The shock tube problem constitutes a particularly useful test case, as it simulates a shock wave, a contact discontinuity and an expansion wave region. Initially, a diaphragm in a long 1D tube separates two initial gas states (both air) at different pressures and densities, but both stationary (u = 0). At t = 0, the diaphragm breaks down suddenly, creating three types of waves. (1) A shock wave moving to the right with speed W , across which all variables jump in value, (2) an expansion wave region moving to the left, across which all variables change in value continuously, and (3) a contact discontinuity (moving with the flow to the right), across which u and P are continuous, but other variables experience a jump. These waves lead to 4 regions in the tube. 14 u2 = u3 W Region 1: Undisturbed portion of driven section, with original values of u1(= 0), P1, ρ1, T1 Region 2: Processed by shock wave propagating through it, P2 > P1, u2 > 0, ρ2 > ρ1. Region 3:Processed by expansion waves (isentropic waves), P3 = P2, u3 = u2, s3 6= s2; thus T3 6= T2, ρ3 6= ρ2. Region 4: Undisturbed protion of driver section, with original P4, ρ4, T4, u4(= 0). For further descriptions and derivation of analytical solution for inviscid flow, see Charles, Hirsch,1990, Numerical computation of internal and external flows, Chapter 16. 1 2 Governing equations ∂U ∂t + ∂F ∂x = J, where U = ρρu ρe , F = ρuρu2 + P (ρe)u , J = 0+∂τxx ∂x −P ∂u ∂x + τxx ∂u ∂x + ∂qx ∂x , Also, equation of state P = ρRT. Heat transferred due to thermal conduction qx = −K∂T ∂x , K = µcp Pr Viscous stress τxx = λ ∂u ∂x + 2µ ∂u ∂x , λ = −2 3 µ, µ = µref ( T Tref )3/2 Tref + 110 T + 110 , where µref = 1.846 is the µ value at reference temperature Tref = 298K. The internal energy is e = ∫ cvdt ≈ cvT ; cp ≡ cp(T ) ≈ const; cv ≡ cv(T ) ≈ const, where µo and To are the initial values of µ and T at t = 0. If the flow is inviscid, qx = τxx = 0. Initial condition: At t = 0, u = 0 throughout the tube. The values of P , ρ, and T = P/ρR jump across the diaphragm, { x < xo, P = P4, ρ = ρ4 x > xo, P = P1, ρ = ρ1 Boundary conditions: Adiabatic walls at x = 0, 10 u = 0, u1 = 0, uN = 0 ∂T ∂x = 0, T1 = 4 3 T2 − 1 3 T3, TN = 4 3 TN−1 − 1 3 TN−2 ∂P ∂x = 0, P1 = 4 3 P2 − 1 3 P3, PN = 4 3 PN−1 − 1 3 PN−2 3 Solution method Use MacCormack scheme. The following scheme yields 2nd-order accuracy in space and time, Predictor: U i = U ti + ∆t ( J ti − F ti − F ti−1 ∆x ) + Sti , Decode u, P , ρ, e, T , set B.C. Construct J, F , S. Corrector: U t+∆ti = 1 2 [ U ti + U i + ∆t ( J i − F i+1 − F i ∆x )] + Si Decode ut+∆t, P t+∆t, ρt+∆t, et+∆t, T t+∆t, set B.C. 2 At boundary, use the appropriate one-sided scheme to calculate ∂F/∂x. For 4th-order accuracy in space, use the following schemes, Predictor ( ∂F ∂x ) i = 7Fi − 8Fi−1 + Fi−2 6∆x Corrector ( ∂F ∂x ) i = −7F i + 8F i+1 − F i+2 6∆x Use central difference for viscous terms (those involving τ , K), put them all in J . Sti is the artificial viscosity, Sti = Cx |P ti+1 − 2P ti + P ti−1| P ti+1 + 2P t i + P t i−1 ( U ti+1 − 2U ti + U ti−1 ) Si = Cx |P i+1 − 2P i + P i−1| P i+1 + 2P i + P i−1 ( U i+1 − 2U i + U i−1 ) 0.01 ≤ Cx ≤ 0.3. Try Cx ≈ 0.1 first. 4 Parameters Domain size Lx = 10m. The diaphragm locates at xo = 5m. Air is used. Assume constant values of cp, cv, γ, Pr equal to the values at T = 298K. In SI units, cp = 1005, cv = 718, γ = cp/cv = 1.4, Pr = 0.707, R = 286.9. Note that µ and K vary with T in the simulation. P4 = 10 5Pa, P1 = 10 4Pa ρ4 = 1kg/m 3, ρ1 = 0.125kg/m 3 PR = P2/P1 = 3.031 Tasks: 1. Compare inviscid (analytical) results (µ = 0, K = 0) with the viscous (numerical) results. 2. Effects of grid resolution: Nx = 21, 101, 1001. 3. Compare 2nd-order (in space) with 4th-order scheme. 5 Report requirement 1. Calculate P , ρ, u, T and plot them vs. x at different times. 2. Plot/Compare numerical (viscous) and analytical (inviscid) results in the same figure. 3. Attach the main part of your program to the report. 4. Submit an electronic copy of your report before/on due date. 5. Report must have: Introduction, Formulation/Numerical Procedure, Results/Discussions, Conclusions, Appendix (code, …) Note: Undergraduate students can work in pairs. Sign both names. 3 6 Reference: Analytical solution for inviscid flow Regions 1 and 4: values of P , T , ρ, u etc. are the same as the initial conditions. Region 2: PR = P2/P1 can be calculated (given here as a known value). Obtained the values of other variables based on PR: PR = P2 P1 → P2 ρ2 ρ1 = 1 + αPR α + PR , α = γ + 1 γ − 1 , γ = cp cv , P2 = ρ2RT2 → ρ2, T2 u2 − u1 a1 = PR − 1 (1 + αPR)1/2 1√ γ(γ−1) 2 , a1 = √ γRT1 → u2 Region 3: P3 = P2, u3 = u2, ρ3 = ( P3 P4 ργ4 )1/γ → P3, ρ3, u3 a4 = a3 + γ − 1 2 u3, a3 = √ γRT3 → T3 Shock wave location (separation of Regions 1 and 2): W − u1 a1 = (PR − 1)a1 γ(u2 − u1) → Shock speed: W xshock = Wt → Location of shock: xshock Expansion wave region (−a4 ≤ x/t ≤ u3 − a3): all variables continuously change values—they are func- tions of x. u = 2 γ + 1 ( a4 + x t ) → u(x) a a4 = ( T T4 )1/2 = 1− γ − 1 2 u a4 → T (x) ρ ρ4 = [ 1− γ − 1 2 ( u a4 )]2/(γ−1) , P = ρRT → ρ(x), P (x) 4
