Lecture 4-5: Nozzle Design: Method of Characteristics

16.512, Rocket Propulsion Prof. Manuel Martinez-Sanchez Lecture 4-5: Nozzle Design: Method of Characteristics The Method of Characteristics (Ideal Gas) (Ref. Phillip Thompson Compressive Fluid Dynamics, McGraw Hill, 1972, Ch. 9) 2-D or axisymmetric Homentropic as well as isentropic ur → ∇×u = o plus ∇ ur (ρ u ) = 0 0 ur ur 1 ⎛ u2 ⎞ ur ur 1 and u ∇u + ∇p = o ⇒ ∇ ⎜ ⎟ + ω × u + ∇p = 0 ρ ρ ⎝ 2 ⎠ Use intrinsic co-ordinates r Eq. of motion along s : u r Eq. of motion along n : u2 1 ∂p ∂u2 2 = =− R ∂n ρ ∂n Now 1 ∂ϑ ≡ R ∂s u2 ∂ϑ ∂u = −u ∂s ∂n ∂u 1 ∂p + =o ∂s ρ ∂s ∂u ∂ϑ +u =o ∂n ∂s (good no p in it) ur (can also get this from ∇ × u = o ) Continuity: 16.512, Rocket Propulsion Prof. M. Martinez-Sanchez ∂ ( ρ u 2πr δn ) = o ∂s (1) (2) 1 ∂ρ 1 ∂u 1 ∂r 1 ∂δn + + + =o ρ ∂s u ∂s r ∂s ∂n ∂s ( δs ) sin ϑ = ∂r 1 ∂ϑ ( ∂s ) sin ϑ ∂s r ∂n ∂δn ⎛ ∂ϑ ⎞ δn ⎟ δs ds = ⎜ − ∂s ⎝ ∂n ⎠ 1 ∂ρ 1 ∂u ∂ϑ sin ϑ + − =− r ρ ∂s u ∂s ∂n Homentropic: dp = c 2 d ρ ⎛ 2 ⎛ ∂p ⎞ ⎞ ⎜⎜ c = ⎜ ⎟ ⎟⎟ ⎝ ∂ρ ⎠s ⎠ ⎝ so and, from s eq. of motion, − 1 ∂p 1 ∂u ∂ϑ sin ϑ + − =− 2 r ρ c ∂s u ∂s ∂n ∂p 1 ∂p ∂u : so you can eliminate = −u ∂s ρ ∂s ∂s u ∂u 1 ∂u ∂ϑ sin ϑ + − =− r c 2 ∂s u ∂s ∂n 1⎛ u2 ⎞ ∂u ∂ϑ sin ϑ 1 − + + = ⎜ ⎟ r u⎝ c 2 ⎠ ∂s ∂n 2 M -1 M 2 − 1 ∂u ∂ϑ sin ϑ + = u ∂s ∂n r Introduce the Mach angle μ = sin−1 tan μ Then (2) 1 = tan−1 M 1 2 M −1 (4) tan μ = M 2 − 1 ∂u ∂ϑ + =o u ∂n ∂s M 2 − 1 ∂u ∂ϑ tan μ sin ϑ + tan μ = u ∂s ∂n r And (4) Introduce the Prandtl-Meyer function ω (M) by dω = M2 − 1 du u (to be integrated later) 1 M2 − 1 tan μ ∂ω ∂ϑ + =o ∂n ∂s then ∂ω ∂ϑ tan μ sin ϑ + tan μ = ∂s ∂n r add and subtract to obtain the “Characteristics form” (single differential operator per equation) ∂ ⎞ tan μ sin ϑ ⎛ ∂ ⎜ ∂s + tan μ ∂n ⎟ ( ω + ϑ ) = r ⎝ ⎠ ⇒ ∂ ⎞ tan μ sin ϑ ⎛ ∂ ⎜ ∂s − tan μ ∂n ⎟ ( ϑ − ω) r ⎝ ⎠ ∂ ∂ ⎞ sin μ sin ϑ ⎛ ⎜ cos μ ∂s + sin μ ∂n ⎟ ( ω + ϑ ) = r ⎝ ⎠ ∂ ∂ ⎞ sin μ sin ϑ ⎛ ⎜ cos μ ∂s − sin μ ∂n ⎟ ( ϑ − ω) = − r ⎝ ⎠ r + ⎧cos μ ⎫ In (s, n) coordinates, 1 = ⎨ ⎬ ⎩ sin μ ⎭ cos μ r+ ∂ ∂ ∂ + sin μ =1 ∇ = ∂s ∂n ∂m+ cos μ r− ∂ ∂ ∂ − sin μ =1 ∇ = ∂s ∂n ∂m− r − ⎧ cos μ ⎫ 1 =⎨ ⎬ ⎩− sin μ ⎭ so (m+, m- are lengths along characteristics) ⎛ m− inclined ϑ + μ ⎞ ⎜⎜ + ⎟⎟ ⎝ m inclined ϑ − μ ⎠ sin μ sin ϑ ⎫ ⎧ ∂ ⎪ ∂m+ ( ϑ + ω) = + ⎪ r ⎪⎪ ⎪⎪ ⎨ ⎬ ⎪ ∂ sin μ sin ϑ ⎪ ⎪ ⎪ ϑ − ω) = − − ( r ⎩⎪ ∂m ⎭⎪ For 2-D, r → ∞ ⇒ ϑ + ω =const. along m+ ≡ I+ (inclined ϑ − μ ) ϑ − ω =const. along m- ≡ I- (inclined ϑ + μ ) 2-D Simple Regions Consider a uniform region; the flow from it enters some disturbed region, like a wall turning. One of the m families originate in the uniform region (m+ in example) and carries a constant invariant, e.g. ϑ + ω = ϑo + ωo everywhere downstream. Along one of the other characteristics (m- here), we carry a constant (which varies from ch. To ch. of that family); we can evaluate it at the wall, for instance ϑ − ω = ϑw − ωw along each m- line and ωw = ϑo + ωo − ϑw so ϑ − ω = 2ϑw − ( ϑo + ωo ) then 2ϑ = ϑo + ωo + ϑw − ωw = 2ϑw along each m− → 2ω = ϑo + ωo − ( ϑw − ωw ) along each m − → ϑ = ϑw ω = ϑo + ωo − ϑw = ωw so ϑ and ω (and hence M, μ) are constant along each m-, and these m- lines are straight (const. ϑ + μ ). This is a Simple Region (one of the invariants is constant). The turning of the flow is dictated by how ϑ changes on a m+ line, as different m- lines are crossed. Since on a m+ we have ϑ + ω = ϑo + ωo , changes of ϑ are equal and opposite to those in ω (ω increases as M does in the expansion, so ϑ decreases (increase negatively) at the same rate.) So, we can interpret ω as the magnitude of the isentropic flow turning in a simple region, i.e., when nothing varies along one characteristic family. Calculation of ω (M) dω = but M2 − 1 du = u ⎛ dM 1 dT ⎞ M2 − 1 ⎜ + 2 T ⎟⎠ ⎝ M dω = γ −1 2 ⎛ M M2 − 1 ⎜ 2 ⎜1 − γ −1 2 M ⎜⎜ M 1+ ⎝ 2 ⎞ ⎟ ⎟ dM = ⎟⎟ ⎠ M2 − 1 M dM γ −1 2 M 1+ 2 Integrates (with ω =0 at M=1) to ω = K tan−1 For M → ∞ M2 − 1 − tan−1 M 2 − 1 K ω→K γ ω (M → ∞ ) ⎛ ⎜⎜ K = ⎝ γ +1⎞ ⎟ γ − 1 ⎟⎠ ⎞ π π π π⎛ γ +1 − = ( K − 1) = ⎜ − 1⎟ ⎜ ⎟ 2 2 2 2 ⎝ γ −1 ⎠ 1.2 208o 1.25 180o 1.4 130.5o 5/3 90o (Max. turning from M=1) So rocket exhaust ( γ 1.2 − 1.3 ) can turn backwards at a sonic nozzle exit to vacuum) (but very long density, long mfp, so continuum approach fails at some point in the expansion, molecules then continue in straight line). Actually, one starts from some high Me>1, so the actual turning is through ω∞ − ω (Me ) , not ω∞ , even in a vacuum: Example of Application: Ideal 2-D plug Nozzle at Design Condition Along m− , ϑ − ω = ϑa − ωa Simple region (Prandtl-Meyer fan centered at lip L) In particular, at inlet M=1, so ϑthroat = −ωa tan μ = From geometry, μ − ϑ + ϕ = Sub. ϑ = ω − ωa ϑt = −ωa ϕ= π + K tan−1 2 r dϕ = dr 1 M2 − 1 π ϑt 2 μ − ω + ωa + ϕ = ϕ= π +ω−μ 2 M2 − 1 − tan−1 M 2 − 1 − tan−1 K − π 2 π + ωa 2 1 2 M −1 ⎛ ⎜⎜ K = ⎝ γ +1⎞ ⎟ γ − 1 ⎠⎟ M2 − 1 K ϕ = K tan−1 tan ⇒ ϕ = K M2 − 1 1 = K K tan μ r dϕ = dr 1 K tan dr =K r ϕ K ϕ⎞ ⎛ ϕ d ⎜ cos ⎟ K ⎠ K dϕ = −K 2 ⎝ ϕ ϕ cos cos K K sin r = For ϕ = 0,r = ht (throat height) In particular, at end of expansion ra = ht K2 ϕa ⎞ ⎛ ⎜ cos K ⎟ ⎝ ⎠ and r = const. K2 ϕ⎞ ⎛ ⎜ cos K ⎟ ⎝ ⎠ ht K2 ϕ⎞ ⎛ ⎜ cos K ⎟ ⎝ ⎠ M = Ma → ϕa = K tan−1 Ma2 − 1 K ha = ra sin μa = , then ra Ma xa = ra cos μa = ra 1 − 1 Ma2 Numerical Application γ Po γ − 1 2 ⎞ γ −1 ⎛ = 100 = ⎜ 1 + ⇒ Ma = Ma ⎟ Pa 2 ⎝ ⎠ Take γ = 1.3 , Then Ma = γ −1 ⎡ ⎤ 2 ⎢⎛ Po ⎞ γ ⎥ ⎜ ⎟ − 1⎥ γ − 1 ⎢⎝ Pa ⎠ ⎥⎦ ⎣⎢ 0.3 ⎤ 2 ⎡ 1.3 ⎢100 − 1⎥ = 3.554 0.3 ⎣ ⎦ ⎛ M2 − 1 ⎞ −1 M2 − 1 ωa = K tan−1 ⎜ ⎟ − tan K ⎝ ⎠ ωa = 67.35o ⇒ K= γ +1 = γ −1 2.3 = 2.769 0.3 ϑt = −67.35o ⎛ M2 − 1 ⎞ o Also ρa = K tan−1 ⎜ ⎟ = 141.01 ⎝ K ⎠ and so ra 1 1 = = K2 2.7692 ht ρa ⎞ 141.01 ⎞ ⎛ ⎛ ⎜ cos 2.769 ⎟ ⎜ cos K ⎟ ⎝ ⎠ ⎝ ⎠ From the geometry, ha=ra sin μa = ra = 34.41 ht ra Ma ha 1 = 34.41 ht 3.554 ha =9.684 ht 2 and xa ⎛ 1 ⎞ = ra cos μa = 34.41 1 − ⎜ ⎟ ht ⎝ 3.554 ⎠ xa =33.02 ht Very long and pointy, should be truncated. Non-Simple Regions. When the characteristics of both families intersect some upstream disturbance, they affect each other’s invariant, and characteristics are no longer straight, and no longer carry constant flow properties (except for their own invariant) 0-1-4 is a simple region 0-4-41 is uniform 1-4-10 is non-simple 4-41--11 is a simple region Then we need to calculate in a step-by-step manner, carrying to each point the two invariants I+, I- from neighboring upstream points, along the m+, m- lines from them to us. After this is done, we know the new segments of m+, m- from our point (slopes ϑ + μ , ϑ − μ ), so we can extend the grid as we go. Notice the flowfield properties can be found first everywhere and only then we need to come back and place the points geometrically. Example: Design a 2-D ideal nozzle to expand from near sonic conditions (M0 = 1.1) to Me = 3. Use only 4 characteristics. Use a corner expansion as a starter, γ =1.25 K = ωo = 3 tan−1 1.21 − 1 − tan−1 1.21 − 1 = 1.435o ; ϑo = 0 3 ωe = 3 tan−1 9 −1 − tan−1 9 − 1 = 59.413o ; ϑe = 0 3 2.25 =3 0.25 At inlet: I − = ϑ − ω = ϑo − ωo = −1.435o (also at 4) At exit: I + = ϑ + ω = ϑe + ωe = 59.413o 7, 9, 10) 2ϑ4 = −1.435 + 59.413 = 57.978 At 4, then, 2ω4 = 59.413 + 1.435 = 60.848 ϑ4 = 28.989o ω4 = 30.424o (also at 4, I1+ = 1.44o Select: I2+ = 20.76o I3+ = 40.09o () () Point M I+ = ϑ + ω 1 1.1 +1.44 2 3 4 (and 4w) 5 6 7 (and 7w) 8 9 (and 9w) 10 (and 11) 1.439 1.726 2.013 1.726 2.013 2.315 2.315 2.641 .3 20.76 40.09 59.41 20.76 40.09 59.41 40.09 59.41 59.41 o I4+ = 59.41o () ϑ -1.44 1.44 -1.44 -1.44 -1.44 -20.76 -20.76 -20.76 -40.09 -40.09 -59.41 11.10 20.77 30.43 20.76 30.43 40.09 40.09 49.75 59.41 I- = ϑ − ω o ω o () () o ϑ+μ 0 65.38 65.38 9.66 19.33 28.99 0 9.67 19.33 0 9.66 0 44.02 35.41 29.79 35.41 29.79 25.59 25.59 22.25 19.47 53.68 54.74 58.78 35.41 39.46 44.92 25.59 31.91 19.47 o μ () o ϑ−μ Note the very shallow angles of the m+ lines (from 4, 7, 9, 10) which will put point 10 far to the right, and point 11 (at slope (m+) of 19.47o even farther. Locating the points geometrically − x a ) tan α + ( x c − x b ) tan β = y a − y b xc = ya − yb + xa tan α + xb tan β tan α + tan β o -65.38o -34.36 -16.08 -0.80 -35.41 -20.12 6.26 -25.59 -12.59 -19.47 Notes Points 1-4: same ϑ − ω Point 5: Here ϑ = 0 (a boundary condition) and ϑ + ω = 20.76 Point 6-7: same ϑ − ω Point 8: Here ϑ = 0 and ϑ + ω = 40.09 Point 9: same ϑ − ω Point 10: Here ϑ = 0 and ϑ − ω = 59.41 (xc () yc − yb = ( xc − xb ) tan β = tan β yc = ( ya − yb ) + ( xa − xb ) tan α tan α + tan β yb tan α + ya tan β + ( xa − xb ) tan α tan β tan α + tan β put xa, ya into 1, 2 tan α into 3 Run 2 → xc (st. in 7) xb, yb into 4.5 Run 3 → yc (st. in 7) tan β into 6 Finding the (x,y)’s is very laborious. Accuracy can be increased by averaging together the angles at the ends of each segment, which can be done because those α + αc β + βc , β = b . angles come from the first pass. For instance, α = a 2 2 Boundary with Prescribed Pressure P= Pc r −1 2⎞ ⎛ ⎜1 + 2 M ⎟ ⎝ ⎠ r , r −1 and M = M ( w ) , so P = P ( w ) (given Po ) or w = w (P ) So, if P is fixed on a boundary (contact surface), we can assign w there (just as we assign θ on a solid boundary) From known point a, θ − w = θa − wa So θ = w (P ) + θa − wa and this determines the slope of the boundary, and that of the “reflected” m+ θ + w = 2w (P ) + θa − wa More General Contact Surface Condition If the outside fluid is also supersonic, we must solve on both sides of the contact surface, making sure P and θ are common at each boundary point P1 ( w1 ) = P2 ( w2 ) θ1 = θ2 From m1+ know I1+ = θ + w1+ (P ) From m2− know I2− = θ − w2− (P ) So, I1+ − I2− = w1+ (P ) + w2− (P ) ⇒ Solve for P ⇒ w1 , w2 and I1+ + 2 I2− =θ+ w1+ (P ) + w2− (P ) known now 2 Modifications for Axisymmetric Conditions ∂ sin μ sin θ ⎫ ⎪ r ⎪ ∂m ⎬ ∂ sin μ sin θ ⎪ ( θ − w) = − ⎪⎭ r ∂m− + ( θ + w) = (1) Calculate ( x, r )c from the angle ( μ − θ )a , ( μ + θ )b ( x − xa ) tan α + ( x − xb ) tan β = ra − rb ⇒x− x a tan α + xb tan β + ra − rb tan α + tan β and r − rb = ( x − xb ) tan β (2) Δm+ = x − xa cos α Δm− = x − xb cos β (3) Advance invariants θ + w , θ − w based on μ, θ at a, b: sin μa sin θa ⎫ Δm+ ⎪ ra ⎪ ⎬ sin μb sin θb − Δm− ⎪ ⎪ rb ⎭ θ + w = ( θ + w )a + θ − w = ( θ − w )b (r at c, computed in (1)) (4) M = M ( w ) , μ = μ (M) (5) Iterate from (1) with α, β averaged between (a, b) and c: α→ ( μ − θ )a + ( μ − θ )c β→ ( μ + θ )b + ( μ + θ )c 2 2 ra + rb 2 ⎫ ⎪ ⎪ ⎪ μa + μb ⎪⎪ μ→ ⎬ better 2 ⎪ ⎪ ⎪ θ + θb ⎪ θ→ a ⎪⎭ 2 r→ (6) Continue to new point. So, computation of ( θ,M) is now coupled to that of ( x, r ) , whereas in 2-D ( θ,M) can be found first. But the actual amount of computation is not much more (only the iteration stops). On the axis, r → 0 , but θ = 0 , so 1. x = xa + ra tan α 2. Δm+ = ra sin α 3. θc = 0; wc = ( θ + w )a + sin μa sin θa ra 4. Mc = M ( wc ) , μc = μ (Mc ) 5. α → ( μ − θ )a + μ c 2 ra , θa , , go to (1) (once) 6. Continue Extension to Cases with Weak Shocks Since the entropy jump in a shock increases only as the cube of ΔPshock ρu20 , the isentropic assumption can be approximately extended when characteristics of our family show some mild convergence (in principle, that is always indicative of shock formation, because they carry conflicting information). When is the convergence too strong? Since characteristics are discretized, for weak convergence the ones of our family will converge, but not cross, and as long as they don’t, it should be OK. Of course, with finer resolution they will cross, but the loss of accuracy in ignoring that is of the same order as that increased by the coarse discretization in the first place. This allows us to calculate off-design nozzle performance, like an overexpanded plug nozzle. 2-D Spike Nozzle with Pa < Padesign Pressure forces from hot-gas bathed surfaces are the same as at design. Net thrust is increased because the Pa contribution PaA e is reduced: ( ) F = Fdes. + Pades. − Pa A e and so Fvac = Fdes. + Pades. A e (just as for a bell) “Practical” Nozzle Designs Ideal nozzle are too long, last portion has small wall angle, so small thrust contribution. With small ?, maybe negative contribution. So, options: (a) Constrain length, ask for contour that gives highest thrust/given L. Methods of calculus of variations (Raw nozzle, Ref : “ Exhaust Nozzle Contour from Optimum Thrust”, Jet Propulsion 28 (June 1958): 377-382. Exit flow non-parallel, non-uniform, computationally high. (b) Ad-hoc method (widely used) is to truncate contract an ideal nozzle. (1) Design wall contour for desired Pe Pc (2) If longer than desired L truncate it to some intermediate area ratio. (3) Contract this truncated nozzle to desired length x L x'a = xb' c = xb' desired xd L ideal (4) Translate profile to right metal kink at P smoothes out. Gives “adequate” performance but less than ? nozzle. Not really justified. (Hoffman, J. D. J. of Propulsion 3 (March-April 1987): 150156. XRS-2200 Engine Data Thrust, lbf At Sea Level In Vacuum 206,500 268,000 Specific Impulse, sec. At Sea Level In Vacuum 339 439 Propellants O2, H2 Mixture Ratio (O/F) 5.5 Chamber Pressure, psia 857 Cycle Gas Generator Area Ratio 58 Throttling, Percent Thrust 40-119 Differential Throttling +/- 15% Dimensions, Inches Forward End 133 high x 88 wide Aft End 46 high x 88 wide Forward to Aft 79