AE/ME 339 Computational Fluid Dynamics (CFD) K. M. Isaac

November 30, 2004

topic19_nozzle_flow

Computational Fluid Dynamics (AE/ME 339)

1

K. M. Isaac MAEEM Dept., UMR

Quasi One-Dimensional Nozzle Flow

Insert Nozzle Flow Schematic Schlieren pictures of the flow

November 30, 2004

topic19_nozzle_flow

2

1

Computational Fluid Dynamics (AE/ME 339)

K. M. Isaac MAEEM Dept., UMR

Consider the steady isentropic flow through a convergent-divergent nozzle shown in the Figure 19.1 The nozzle is attached to a large reservoir in which the conditions remains constant, denoted by pressure p0 and temperature T0. p0 and T0 are known as the stagnation pressure and stagnation temperature, respectively. If the nozzle exhausting into a region in which the pressure is much smaller (exact pressure ratios will be established as we proceed) compared to the reservoir pressure. The flow can be analyzed as a quasi one-dimensional flow (flow properties are a function of the axial location only; there is no variation in the radial or azimuthal direction). November 30, 2004

topic19_nozzle_flow

Computational Fluid Dynamics (AE/ME 339)

3

K. M. Isaac MAEEM Dept., UMR

Under the right conditions, the flow accelerates gradually along the nozzle axis and attains a Mach number greater than one at the exit. The flow variables are given by close form solutions which can be used to compare any numerical results.

November 30, 2004

topic19_nozzle_flow

4

2

Computational Fluid Dynamics (AE/ME 339)

K. M. Isaac MAEEM Dept., UMR

Isentropic one-dimensional flow We assume that the flow variables can be represented as a function of the axial coordinate alone. For steady flow we can write ρ Au = constant where A is the nozzle cross-sectional area Taking the logrithm and differentiating yields dρ

dA du + =0 ρ A u The momentum equation can be written as d

+

u 2 dp + =0 2 ρ

November 30, 2004

topic19_nozzle_flow

Computational Fluid Dynamics (AE/ME 339)

5

K. M. Isaac MAEEM Dept., UMR

Rewriting the momentum equation gives d ρ dp dρ = udu + a 2 =0 udu + ρ dρ ρ dρ substitute for from continuity equation

ρ

⎛ du dA ⎞ udu − a 2 ⎜ + ⎟=0 A⎠ ⎝ u dA u 2 du du du = 2 − = ( M 2 − 1) A a u u u November 30, 2004

topic19_nozzle_flow

6

3

Computational Fluid Dynamics (AE/ME 339)

K. M. Isaac MAEEM Dept., UMR

The above equation can be rewritten as 1 dA ⎛ M 2 − 1 ⎞ du =⎜ ⎟ A dx ⎝ u ⎠ dx

x

This equation can be used to discuss nozzle shape Case 1: subsonic region (M < 1, (M 2 - 1) < 0)

November 30, 2004

topic19_nozzle_flow

Computational Fluid Dynamics (AE/ME 339)

7

K. M. Isaac MAEEM Dept., UMR

dA < 0), then the dx du flow will accelerate because > 0. dx Similar arguments indicate that when the Mach number is > 1, the flow will accelerate in a diverging passage. If the passage is converging (

The opposite will be true for decelerating flow.

November 30, 2004

topic19_nozzle_flow

8

4

Computational Fluid Dynamics (AE/ME 339)

K. M. Isaac MAEEM Dept., UMR

Adiabatic relation for an ideal gas T0 ⎡ γ − 1 2 ⎤ = 1+ M ⎥ T ⎢⎣ 2 ⎦ Isentropic relations for an ideal gas γ

p 0 ⎡ γ − 1 2 ⎤ γ −1 = 1+ M ⎥ p ⎢⎣ 2 ⎦ 1

ρ0 ⎡ γ − 1 2 ⎤ γ −1 = 1+ M ⎥ 2 ρ ⎢⎣ ⎦ November 30, 2004

topic19_nozzle_flow

Computational Fluid Dynamics (AE/ME 339)

9

K. M. Isaac MAEEM Dept., UMR

Equations of Nozzle Flow. Mass flow rate: m = ρ Au Can be transformed using isentropic relation as: 1/ γ

⎛ p⎞ ρ m = ρ 0 uA = ρ0 uA ⎜ ⎟ ρ0 ⎝ p0 ⎠

Substitution for the throat conditions gives m =

A p0 ⎛ 2 ⎞ γ⎜ ⎟ RT0 ⎝ γ +1⎠ *

November 30, 2004

γ −1 γ

topic19_nozzle_flow

10

5

Computational Fluid Dynamics (AE/ME 339)

K. M. Isaac MAEEM Dept., UMR

The area ratio can be expressed in terms of the Mach number as follows γ +1

1 ⎡ 2 ⎛ γ − 1 2 ⎞ ⎤ γ −1 ⎛ A⎞ M ⎟⎥ ⎜ *⎟ = 2⎢ ⎜1 + 2 M ⎣γ +1⎝ ⎝A ⎠ ⎠⎦ Note that the area ratio between any given location in the nozzle and the throat is only a function of the Mach number. Once the 2

area ratio is known, the corresponding Mach number can be determined. Note also that the flow must be isentropic (i. e., no shocks present in the nozzle passage) for these results to hold. November 30, 2004

topic19_nozzle_flow

Computational Fluid Dynamics (AE/ME 339)

11

K. M. Isaac MAEEM Dept., UMR

Insert Figure 7.2 November 30, 2004

topic19_nozzle_flow

12

6

Computational Fluid Dynamics (AE/ME 339)

K. M. Isaac MAEEM Dept., UMR

CFD Solution of Nozzle Flow

Insert Figure 7.3 November 30, 2004

topic19_nozzle_flow

Computational Fluid Dynamics (AE/ME 339)

Continuity: Momentum: Energy:

November 30, 2004

13

K. M. Isaac MAEEM Dept., UMR

∂ ∂ ( ρ A) + ( ρ Au ) = 0 ∂t ∂x ∂ ∂ ∂p ( ρ uA) + ( ρ u 2 A) = − A ∂t ∂x ∂x ∂ ∂ ∂ ( ρ eA) + ( ρ euA) = − p ( Au ) ∂t ∂x ∂x

topic19_nozzle_flow

14

7

Computational Fluid Dynamics (AE/ME 339)

K. M. Isaac MAEEM Dept., UMR

Equations in non-conservation form for ideal gas ∂ ∂u ∂ρ =0 Continuity: ( ρ A) + ρ A + uA ∂t ∂x ∂x ∂u ∂u ∂ρ ⎞ ⎛ ∂T + ρu = −R ⎜ ρ +T Momentum: ρ ⎟ ∂t ∂x ∂x ⎠ ⎝ ∂x ∂T ∂T ∂ (ln A) ⎤ ⎡ ∂u ρ cv + ρ uc v = − ρ RT ⎢ + u Energy: ∂t ∂x ∂x ⎥⎦ ⎣ ∂x

November 30, 2004

topic19_nozzle_flow

Computational Fluid Dynamics (AE/ME 339)

15

K. M. Isaac MAEEM Dept., UMR

The finite difference form We divide the nozzle axis (x-axis) into a number of segments, I, each of length ∆x. Therefore, beginning with first node (1), we will have a total of (I+1) nodes. Nodes 1 and (I+1) are the boundary nodes. MacCormack's method can now be used for the numerical solution. Applying MacCormack's method to the non-conservation form yields the following:

November 30, 2004

topic19_nozzle_flow

16

8

Computational Fluid Dynamics (AE/ME 339)

K. M. Isaac MAEEM Dept., UMR

⎛ un − un ⎞ ⎛ ρn − ρn ⎞ n − ln A ⎞ ⎛ ln A ⎛ ∂ρ ⎞ n n n n ⎜ ⎟ + + 1 1 i i i i −ρ u ⎜ ⎟ − u ⎜ i +1 i ⎟ ⎜ ⎟ = − ρi ⎜ ⎟ ⎟ i i ⎜ ⎟ i ⎜⎜ ∆x ∆x ∆x ⎝ ∂t ⎠i ⎜ ⎟ ⎟ ⎝ ⎠ ⎝ ⎠ ⎝ ⎠ n n⎞ ⎛ un − un ⎞ n⎛ n −T ⎞ ⎛T ⎛ T ⎞ ⎜ ρi + 1 − ρi ⎟ ⎛ ∂u ⎞ n ⎜ ⎟ + + 1 1 i i i i − R⎜ ⎟ − R⎜ ⎟ ⎜ ⎟ = −ui ⎜ ⎟ ⎟ ⎜ ⎟ ∆x ∆x ∆x ⎝ ∂t ⎠i ⎝ ρ ⎠i ⎜⎜ ⎜ ⎟ ⎟ ⎝ ⎠ ⎝ ⎠ ⎝ ⎠ ⎛ T n −T n ⎞ ⎡ un − un n ln An − ln An ⎤ ⎛ ∂T ⎞ n n ⎢ i +1 n ⎜ ⎟ i i i i +1 i ⎥ + 1 − (γ − 1)Ti +u ⎜ ⎟ = −ui ⎜ ⎟ ⎢ ⎥ i ∂ x x x t ∆ ∆ ∆ ⎝ ⎠i ⎜ ⎟ ⎝ ⎠ ⎣⎢ ⎦⎥

November 30, 2004

topic19_nozzle_flow

Computational Fluid Dynamics (AE/ME 339)

17

K. M. Isaac MAEEM Dept., UMR

Predictor Step n +1 i

⎛ ∂ρ ⎞ = ρ +⎜ ⎟ ∆t ∂ t ⎝ ⎠i

n +1 i

⎛ ∂u ⎞ = u + ⎜ ⎟ ∆t ⎝ ∂t ⎠i

n

ρ

n i

n

u

n i

⎛ ∂T ⎞ = Ti + ⎜ ⎟ ∆t ⎝ ∂t ⎠i n

Ti

November 30, 2004

n +1

n

topic19_nozzle_flow

18

9

Computational Fluid Dynamics (AE/ME 339)

K. M. Isaac MAEEM Dept., UMR

n n ⎞ ⎛ un −un ⎞ ⎛ ln A − ln A ⎞ n ⎛⎜ ρ − ρ n n n ⎜ ⎟ − − −1 ⎟ 1 1 i i i i i i = −ρ −ρ u ⎜ ⎟−u ⎟ ⎜ ⎟ i ⎜⎜ i i i ⎜ ⎟ ∆x ∆x ∆x ⎟ ⎜ ⎟ ⎝ ⎠ ⎝ ⎠ ⎝ ⎠ n n ⎞ ⎛un −un ⎞ n⎛ n +1 ⎛ T −T ⎞ ⎛ T ⎞ ⎜ ρi − ρi − 1 ⎟ ⎛ ∂u ⎞ n ⎜ ⎟ 1 1 i i i i − − ⎟ − R⎜ ⎟ − R⎜ ⎜ ⎟ = −ui ⎜ ⎟ ⎟ ⎜ ⎟ ∆x ∆x ∆x ⎝ ∂t ⎠i ⎝ ρ ⎠i ⎜⎜ ⎜ ⎟ ⎟ ⎝ ⎠ ⎝ ⎠ ⎝ ⎠ ⎛ T n −T n ⎞ ⎡u n − u n n +1 ln An − ln An ⎤ ⎛ ∂T ⎞ n ⎢ i n n ⎜ ⎟ i i i i i −1 ⎥ − − 1 1 − (γ − 1)Ti +u ⎜ ⎟ = −ui ⎜ ⎟ ⎢ ⎥ i x x x t ∆ ∆ ∆ ∂ ⎝ ⎠i ⎜ ⎟ ⎝ ⎠ ⎣⎢ ⎦⎥ n +1

⎛ ∂ρ ⎞ ⎜ ⎟ ⎝ ∂t ⎠i

November 30, 2004

topic19_nozzle_flow

Computational Fluid Dynamics (AE/ME 339)

19

K. M. Isaac MAEEM Dept., UMR

Find average of the derivatives n n+1 1 ⎡⎛ ∂ρ ⎞ ⎛ ∂ρ ⎞ ⎤ ⎛ ∂ρ ⎞ ⎜ ⎟ = ⎢⎜ ⎟ +⎜ ⎟ ⎥ ⎝ ∂t ⎠av 2 ⎢⎣⎝ ∂t ⎠i ⎝ ∂t ⎠i ⎦⎥ n n +1 1 ⎡⎛ ∂u ⎞ ⎛ ∂u ⎞ ⎤ ⎛ ∂u ⎞ ⎜ ⎟ = ⎢⎜ ⎟ + ⎜ ⎟ ⎥ ⎝ ∂t ⎠ av 2 ⎢⎣⎝ ∂t ⎠i ⎝ ∂t ⎠ i ⎥⎦ n +1 n 1 ⎡⎛ ∂T ⎞ ⎛ ∂T ⎞ ⎤ ⎛ ∂T ⎞ ⎜ ⎟ = ⎢⎜ ⎟ +⎜ ⎟ ⎥ ⎝ ∂t ⎠ av 2 ⎢⎣⎝ ∂t ⎠i ⎝ ∂t ⎠ i ⎥⎦ November 30, 2004

topic19_nozzle_flow

20

10

Computational Fluid Dynamics (AE/ME 339)

K. M. Isaac MAEEM Dept., UMR

Corrected values, time level (n+1) n +1

n +1 i

⎛ ∂ρ ⎞ = ρ +⎜ ⎟ ∆t ∂ t ⎝ ⎠ av

n +1 i

⎛ ∂u ⎞ = u + ⎜ ⎟ ∆t ⎝ ∂t ⎠ av

ρ

n i

n +1

u

n i

n +1

Ti

n +1

⎛ ∂T ⎞ = Ti + ⎜ ⎟ ∆t ⎝ ∂t ⎠av

November 30, 2004

n

topic19_nozzle_flow

Computational Fluid Dynamics (AE/ME 339)

21

K. M. Isaac MAEEM Dept., UMR

∆t , can be used to determine ∆x the time step that satisfies the stability criterion. C ≤ 1 ensures that the solution will be stable. Note the above stability criterion

Courant Number

C = (u + a )

is an empirical relation and therefore, it would be safe to use a time step for which the value of C will be well below unity.

November 30, 2004

topic19_nozzle_flow

22

11

Computational Fluid Dynamics (AE/ME 339)

K. M. Isaac MAEEM Dept., UMR

Boundary conditions The boundary conditions can be established using the theory of the method of characteristics (Tannehill, et al. 1997). At the subsonic inflow boundary we set the values of T and ρ and use linear extrapolation for u. At the nozzle exit, T, ρ , and u are obtained by using linear extrapolation. November 30, 2004

topic19_nozzle_flow

23

Program Completed University of Missouri-Rolla Copyright 2002 Curators of University of Missouri November 30, 2004

topic19_nozzle_flow

24

12