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)

(1)

(2)

ur (can also get this from ∇ × u = o )

Continuity:

16.512, Rocket Propulsion Prof. M. Martinez-Sanchez

∂ ( ρ u 2πr δn ) = o ∂s

Lecture 8 Page 1 of 20

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 ρ

16.512, Rocket Propulsion Prof. M. Martinez-Sanchez

⎛ 2 ⎛ ∂p ⎞ ⎞ ⎜⎜ c = ⎜ ⎟ ⎟⎟ ⎝ ∂ρ ⎠s ⎠ ⎝

Lecture 8 Page 2 of 20

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 μ =

1 M2 − 1

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

16.512, Rocket Propulsion Prof. M. Martinez-Sanchez

du u

(to be integrated later)

Lecture 8 Page 3 of 20

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 ϑ + μ )

16.512, Rocket Propulsion Prof. M. Martinez-Sanchez

Lecture 8 Page 4 of 20

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ω =

M2 − 1

du = u

⎛ dM 1 dT ⎞ M2 − 1 ⎜ + 2 T ⎟⎠ ⎝ M

but 16.512, Rocket Propulsion Prof. M. Martinez-Sanchez

Lecture 8 Page 5 of 20

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:

16.512, Rocket Propulsion Prof. M. Martinez-Sanchez

Lecture 8 Page 6 of 20

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, μ − ϑ + ϕ =

ϑt = −ωa

ϕ=

π + K tan−1 2

ϕ=

M2 − 1

−

π + ωa 2

π +ω−μ 2

M2 − 1 − tan−1 M 2 − 1 − tan−1 K

16.512, Rocket Propulsion Prof. M. Martinez-Sanchez

1

π ϑt 2

μ − ω + ωa + ϕ =

Sub. ϑ = ω − ωa

r dϕ = dr

1 2

M −1

⎛ ⎜⎜ K = ⎝

γ +1⎞ ⎟ γ − 1 ⎠⎟

π 2

Lecture 8 Page 7 of 20

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 −

16.512, Rocket Propulsion Prof. M. Martinez-Sanchez

1 Ma2

Lecture 8 Page 8 of 20

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 Ma

ha 1 = 34.41 ht 3.554

16.512, Rocket Propulsion Prof. M. Martinez-Sanchez

ra = 34.41 ht

ha =9.684 ht

Lecture 8 Page 9 of 20

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 16.512, Rocket Propulsion Prof. M. Martinez-Sanchez

Lecture 8 Page 10 of 20

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

16.512, Rocket Propulsion Prof. M. Martinez-Sanchez

(also at 4,

ϑ4 = 28.989o ω4 = 30.424o

Lecture 8 Page 11 of 20

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 =

16.512, Rocket Propulsion Prof. M. Martinez-Sanchez

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

()

ya − yb + xa tan α + xb tan β tan α + tan β

Lecture 8 Page 12 of 20

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

16.512, Rocket Propulsion Prof. M. Martinez-Sanchez

Lecture 8 Page 13 of 20

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

16.512, Rocket Propulsion Prof. Manuel Martinez-Sanchez

Lecture 4-5 Page 14 of 20

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) =

16.512, Rocket Propulsion Prof. Manuel Martinez-Sanchez

Lecture 4-5 Page 15 of 20

(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

16.512, Rocket Propulsion Prof. Manuel Martinez-Sanchez

Lecture 4-5 Page 16 of 20

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. 16.512, Rocket Propulsion Prof. Manuel Martinez-Sanchez

Lecture 4-5 Page 17 of 20

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.

16.512, Rocket Propulsion Prof. Manuel Martinez-Sanchez

Lecture 4-5 Page 18 of 20

Gives “adequate” performance but less than ? nozzle. Not really justified. (Hoffman, J. D. J. of Propulsion 3 (March-April 1987): 150156.

16.512, Rocket Propulsion Prof. Manuel Martinez-Sanchez

Lecture 4-5 Page 19 of 20

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

16.512, Rocket Propulsion Prof. Manuel Martinez-Sanchez

Lecture 4-5 Page 20 of 20