Solving for the Exit Flow of a Supersonic Nozzle

Open Live Script

This example shows how to use the method of characteristics and Prandtl-Meyer flow theory to solve a problem in supersonic flow involving expansions. Solve for the flow field downstream of the exit of a supersonic nozzle.

Problem Definition

This section describes the problem to be solved. It also provides necessary equations and known values.

Solve for the flow field downstream of a supersonic nozzle using the method of characteristics. The Mach number at the exit plane is 1.5 and the pressure at the exit plane is 200 kilopascals. The back pressure is 100 kilopascals.

Assumptions:

Flow is isentropic.
Variation in flow properties depend on the interaction of expansion waves that occur throughout the wake of the nozzle.
The geometry of the nozzle and the flow is symmetric.

Model the expansion fan as three characteristics. Due to symmetry, arbitrarily choose to work only on the top half of the flow. Following is a figure of the nozzle exit.

upperNozzle = plotExpansionSchematic('uppernozzle');

Figure contains an axes object. The axes object contains 25 objects of type line, text.

The given information in the problem is:

exitMach = 1.5; % Mach number at the exit plane [dimensionless]
exitPres = 200; % Static pressure at the exit plane [kPa]
backPres = 100; % Pressure downstream of the nozzle, outside of the expansion wake

The fluid is assumed to be air that behaves like a perfect gas with the following constant specific heat ratio.

k = 1.4; % Specific heat ratio [dimensionless]

Method of Characteristics

The method of characteristics is a theory for supersonic flow that analyzes the irrotational potential flow equation in fully nonlinear form. Isentropic flow is assumed. The definition of characteristics are the curves in the flow where the velocity is continuous but the first derivative of velocity is discontinuous.

The blue lines in the previous figure are approximate characteristics. Characteristics of type I make a negative acute angle with the flow direction. Characteristics of type II make a positive acute angle with the flow direction. A detailed derivation of the method is outside the scope of this example analysis. This example analysis uses the region-to-region procedure. It is assumed that you are familiar with this procedure.

In Prandtl-Meyer flow and the method of characteristics, calculate the important angles for all regions of the flow.

Flow angle is the direction in which the air is moving.

$θ_{n} = F l o w a n g l e i n t h e n^{t h} r e g i o n$

The Prandtl-Meyer angle is the angle at which the flow changes direction from one region to another.

$ν_{n} = P r a n d t l - M e y e r a n g l e i n t h e n^{t h} r e g i o n$

Mach angle is the angle between the local flow direction and the weak pressure waves that emanate from a given point.

$μ_{n} = M a c h a n g l e i n t h e n^{t h} r e g i o n$

Compute the Mach number in each region and solve for the angles of both types of characteristic in all of the regions. Solve for the geometric boundary of all of the regions by calculating the slopes of all of the characteristics and location of all of the intersections of characteristics.

Computing the Flow Properties Through the First Expansion Fan

Determine the Mach number outside the wake (region 4). The Mach number at this location can be found using isentropic ratios for pressure and the given values for pressure. The pressure ratio at the exit plane is easily solved for using flowisentropic.

[~, ~, exitPresRatio] = flowisentropic(k, exitMach);

The back pressure ratio is the ratio of the back pressure to the stagnation pressure. The isentropic pressure ratio at the outer wake region is:

$\frac{p_{4}}{p_{0}} = \frac{p_{4}}{p_{1}} \frac{p_{1}}{p_{0}}$

backPresRatio = backPres / exitPres * exitPresRatio;

Calculate the Mach number in region 4 using flowisentropic.

backMach = flowisentropic(k, backPresRatio, 'pres');

The 'pres' string input indicates that the function is in pressure ratio input mode. The flow angle in the back pressure region is the difference in Prandtl-Meyer angles from the exit plane region (region 1) to the back pressure region (region 4).

$θ_{4} = ν_{4} - ν_{1}$

[~, nu_1]    = flowprandtlmeyer(k, exitMach);
[~, nu_4]    = flowprandtlmeyer(k, backMach);
theta_4 = nu_4 - nu_1;

Because we are approximating the flow with three characteristics, calculate the change in flow angle in crossing type I characteristics from region 1 to region 4:

$Δ θ_{I} = \frac{θ_{4}}{3}$

deltaThetaI = theta_4 / 3;

Note that flow in region 1 is parallel to the horizontal and therefore:

theta_1 = 0;

In fact, the flow in any region that straddles the centerline is parallel to the centerline. This is because the centerline is considered to be a boundary for this symmetric flow. In addition, there are no sources or sinks at the boundary.

theta_5  = 0;
theta_8  = 0;
theta_10 = 0;

The flow angles of regions 2 and 3 follow simply.

theta_2 = theta_1 + deltaThetaI;
theta_3 = theta_2 + deltaThetaI;

Across type I characteristics, the change in Prandtl-Meyer angle equals to the change in flow angle:

$Δ ν_{I} = Δ θ_{I}$

deltaNuI = deltaThetaI;

Calculate the Prandtl-Meyer angle in region 2 by using the Prandtl-Meyer angle in region 1 and deltaNuI, the change in Prandtl-Meyer angle through type I characteristics. Calculate the Prandtl-Meyer angle in region 3 in a similar manner to that of region 2.

$Δ ν_{I} = ν_{2} - ν_{1}$

$ν_{2} = ν_{1} + Δ ν_{I}$

nu_2 = nu_1 + deltaNuI;
nu_3 = nu_2 + deltaNuI;

Calculating Flow Properties in the Interference Regions

The flow angle in region 5 is known to be zero from the centerline boundary condition. Therefore, the change in angle from region 2 to region 5 is

$Δ θ_{I I} = θ_{5} - θ_{2}$

deltaThetaII = theta_5 - theta_2;

Calculate the change in Prandtl-Meyer angle across type II characteristics:

$Δ ν_{I I} = - Δ θ_{I I}$

deltaNuII = -deltaThetaII;

Then, calculate the Prandtl-Meyer angle in region 5. You already know the region 2 Prandtl-Meyer angle and the change in Prandtl-Meyer angle across type II characteristics.

nu_5 = nu_2 + deltaNuII;

To calculate the properties in region 6, use the fact that the properties in region 3 and region 5 are known. Note also that which characteristic the flow crosses define the changes in properties. From region 5 to region 6, a type I characteristic is crossed. Therefore,

$Δ θ_{I} = Δ ν_{I} = θ_{6} - θ_{5} = ν_{6} - ν_{5}$

Rearranging this as:

$ν_{6} - θ_{6} = ν_{5} - θ_{5} (1)$

A type II characteristic is crossed in going from region 3 to region 6. Therefore,

$Δ ν_{I I} = - Δ θ_{I I}$

$ν_{6} - ν_{3} = θ_{3} - θ_{6}$

Rearranging this as:

$ν_{6} + θ_{6} = ν_{3} + θ_{3} (2)$

Add equations (1) and (2) together, then solve for the Prandtl-Meyer angle in region 6. This yields the following expression.

$ν_{6} = \frac{(ν_{5} - θ_{5}) + (ν_{3} + θ_{3})}{2}$

In MATLAB®, use:

nu_6 = ((nu_5 - theta_5) + (nu_3 + theta_3))/2;

From equation (1), the flow angle in region 6 is

theta_6 = nu_6 - (nu_5 - theta_5);

For region 7, a type one characteristic is crossed and all information is available in region 6.

nu_7    = nu_6 + deltaNuI;
theta_7 = theta_6 + deltaThetaI;

Region 8 is on the centerline; its flow angle is zero. Going from region 6 to region 8 requires crossing a type II characteristic. Therefore, calculate the Prandtl-Meyer angle in region 8 as:

nu_8 = nu_6 + deltaNuII;

Calculate the Prandtl-Meyer angle and flow angle in region 9 as you did for region 6. Region 8 is the upstream region across the type I characteristic. Region 7 is the upstream region across the type II characteristic.

nu_9    = ((nu_8 - theta_8) + (nu_7 + theta_7))/2;
theta_9 = nu_9 - (nu_8 - theta_8);

Region ten is on the centerline. The flow is parallel and so the flow angle is zero. Use the Prandtl-Meyer angle in region 9 and the crossing over a type II characteristic to calculate the Prandtl-Meyer angle in region 10.

nu_10 = nu_9 + deltaNuII;

Preparing and Tabulating the Flow Parameter Results

For upcoming calculations, combine the flow angles into one vector and the Prandtl-Meyer angles into another vector.

flowAngles         = [theta_1; theta_2; theta_3; theta_4; theta_5; theta_6; theta_7; theta_8; theta_9; theta_10];

prandtlMeyerAngles = [nu_1; nu_2; nu_3; nu_4; nu_5; nu_6; nu_7; nu_8; nu_9; nu_10];

To calculate the Mach numbers and Mach angles in each region, using the flowprandtlmeyer function with the prandtlMeyerAngles as the input. You can use the results from this function to find the angle that the type I and type II characteristics make with the horizontal inside each region. You can then use these angles to calculate slopes in the x-y plane, where the centerline is the x-axis and the exit plane of the nozzle is the y-axis. For type I characteristics and type II characteristics, respectively, the slopes are:

${(\frac{d y}{d x})}_{I} = t a n (θ - μ)$

${(\frac{d y}{d x})}_{I I} = t a n (θ + μ)$

Note, the values in the following table for type I and type II are the angles with the horizontal, not the slopes.

% Preallocation for speed
machNumbers = zeros(size(flowAngles));
machAngles  = zeros(size(flowAngles));

for i = 1:numel(flowAngles)
    [machNumbers(i), ~, machAngles(i)] = flowprandtlmeyer(k, prandtlMeyerAngles(i), 'nu');
end

typeOne = flowAngles - machAngles;
typeTwo = flowAngles + machAngles;

m = (1:numel(machNumbers)).';
disp(table(m,round(flowAngles,2),round(prandtlMeyerAngles,2),round(machNumbers,3),round(machAngles,2),round(typeOne,2),round(typeTwo,2),...
    'VariableNames',["Region","theta (Deg)","nu (Deg)","Mach","mu (Deg)","type I (Deg)","type II (Deg)"]))

    Region    theta (Deg)    nu (Deg)    Mach     mu (Deg)    type I (Deg)    type II (Deg)
    ______    ___________    ________    _____    ________    ____________    _____________

       1             0        11.91        1.5     41.81         -41.81           41.81    
       2          4.45        16.35       1.65     37.29         -32.85           41.74    
       3          8.89         20.8      1.803      33.7          -24.8           42.59    
       4         13.34        25.24      1.959     30.69         -17.35           44.03    
       5             0         20.8      1.803      33.7          -33.7            33.7    
       6          4.45        25.24      1.959     30.69         -26.25           35.14    
       7          8.89        29.69      2.122     28.11         -19.22              37    
       8             0        29.69      2.122     28.11         -28.11           28.11    
       9          4.45        34.14      2.294     25.84          -21.4           30.29    
      10             0        38.58      2.477     23.81         -23.81           23.81

Note the following:

The flow angles increase away from the centerline.
The Prandtl-Meyer angles increase as the flow moves downstream.
The Mach number also increases as the flow moves downstream.

Solving for the Flow Geometry

The flow properties are known in all regions, but to solve for the flow field, you must calculate the actual geometry of each region. The last two columns of the above table contain the angles that each type of characteristic makes with the horizontal. Because straight lines approximate the characteristics of the flow in each region, the boundary between any two regions is approximated by the average of the angles that each makes in the bordering regions. Because the waves bend through the expansion fan, begin the analysis from the point from which the characteristics originate. The characteristics originate at the lip of the nozzle and work downstream.

Assume the intersection of the centerline and the exit plane of the nozzle is the origin of our coordinate system. Also assume lengths are normalized to half of the exit height of the nozzle. The positive x-axis is taken horizontal along the centerline in the downstream direction. The positive y-axis is vertically up in the exit plane of the nozzle. The lip of the nozzle is at the point (0,1).

All three characteristics that propagate from the upper lip are type I characteristics. Analyze the steepest sloping characteristic first because no waves interfere with the steepest wave until the steepest wave intersects the centerline. In the symmetric half-nozzle model, the steepest wave reflects back into the fan and interferes with the other waves in the expansion fan.

This wave that "reflects" from the centerline is actually the steepest sloping type II characteristic that propagates from the bottom lip. However, the analysis considers the centerline to be a boundary due to symmetry. This produces the same results that you would get if you worked both halves of the nozzle.

The steepest sloping line from the lip is a type I characteristic that separates region 1 and region 2. To calculate the angle that the steepest sloping wave makes with the horizontal, use the average of angles that the type I characteristics make in each region. To calculate the slope, use trigonometry.

avgAngle12 = (typeOne(1) + typeOne(2)) / 2;
slope12    = tand(avgAngle12);

With the following information known:

Slope of the first type I wave in x-y space.
The y-intercept of the wave ( y = 1 at the lip).
The wave intersects the centerline (y = 0) without interference.

Calculate the x-location of the point using the equation of a line in slope-intercept form. Rearrange y = m*x+b for the x-location with y = 0 to produce x = -b / m. This is the x-location of the first downstream point, point 1.

y1 = 0; % On the centerline
x1 = -1 / slope12;

From point 1, the first type II characteristic propagates and interferes with the fan. The other type I characteristics that originate from the nozzle lip are disturbed by the type II wave, but not before reaching that wave. Therefore, calculate the points of intersection of the steepest type II characteristic and the flatter type I waves from the lip. The type II characteristic coming up from the centerline separates region 2 and region 5. The average of the two angles and associated slopes are given by:

avgAngle25 = (typeTwo(2) + typeTwo(5)) / 2;
slope25    = tand(avgAngle25);

The second steepest type I characteristic is the from the nozzle lip, separating region 2 and region 3. The average angle with the horizontal and the associated slope of the wave are given by:

avgAngle23 = (typeOne(2) + typeOne(3)) / 2;
slope23    = tand(avgAngle23);

Calculate the point of intersection of the region 2-3 boundary and the region 2-5 boundary. You need this point because the characteristics interfere with each other at this point. The slopes of both boundaries and a point on each line are known. Point 1 and the nozzle lip (to be reference point 0) are known. Solve for the unknown x-coordinate of the point of intersection. Use that x-location in the equation of either of the two lines to find the y-location of the point of intersection. The point-slope form equation of a line through point p with slope m is:

$y - y_{p} = m (x - x_{p})$

The advantage of this form of a line is that you need only one point and the slope to completely define the line. The x and y without subscripts can be any point on the line. However, the point of intersection of two lines must be unique. Calling this point of intersection point 2, the equation of both lines are the following.

$y_{2} - y_{0} = m_{2, 3} (x_{2} - x_{0}) (3)$

$y_{2} - y_{1} = m_{2, 5} (x_{2} - x_{1}) (4)$

where

$m_{i, j} = S l o p e o f t h e b o u n d a r y b e t w e e n r e g i o n i a n d r e g i o n j$

Subtract and rearrange:

$x_{2} = \frac{x_{1} m_{2, 5} - x_{0} m_{2, 3} + y_{0} - y_{1}}{m_{2, 5} - m_{2, 3}}$

Knowing some of the values exactly due to the axes intercepts simplifies this expression to:

x2 = (x1 * slope25 + 1) / (slope25 -slope23);

Below the y-location of point 2 is found by plugging the x-location of point 2 into equation (4) above, but plugging into equation (3) works just as well.

y2 = (x2 - x1) * slope25;

Use the slope-intercept formula and the procedure above to calculate all points. To calculate the third point in the flow, first calculate the intersection of the region 3-4 boundary and the 3-6 boundary. The angles of the boundary lines are calculated using the average of the angles with the horizontal. You can then use trigonometry to find the slope, which is now computed in one step.

slope34 = tand( (typeOne(3) + typeOne(4)) / 2 );
slope36 = tand( (typeTwo(3) + typeTwo(6)) / 2 );

Because the boundary between region 3 and region 4 is a type I characteristic and the boundary between region 3 and region 6 is a type II characteristic, be careful to take the angles for the appropriate type. Use the point-slope form of these boundary lines subtracted from each other to calculate the x-location of point 3.

x3 = (y2 - 1 - x2 * slope36) / (slope34 - slope36);
y3 = (x3 - x2) * slope36 + y2;

The first type II characteristic now propagates beyond the fan and does not interfere with any other characteristics. The angle that defines the direction in which the steepest type II propagates is an angle which is the average of type II waves in regions 4 and 7.

slope47 = tand( (typeTwo(4) + typeTwo(7)) / 2);

Solve also the second steepest type I characteristic to continue downstream. Start from the known location of point 2 to calculate the x-intercept of the second type I characteristic. Again, the solution uses the point-slope form of the line. However, there is no interference until the centerline boundary is reached. We only need to consider one line to find the x-location of "point 4". The slope of the boundary between region 5 and region 6 is a type I wave:

slope56 = tand( (typeOne(5) + typeOne(6)) / 2 );

The rearranged point-slope form (knowing y = 0 at the centerline) is used to find point 4.

x4 = ( slope56 * x2 - y2 ) / slope56;
y4 = 0;

Calculate point 5 in the same manner as point 2. The region 6-7 boundary is type I and the region 6-8 boundary is type II.

slope67 = tand( ( typeOne(6) + typeOne(7)) / 2);
slope68 = tand( (typeTwo(6) + typeTwo(8)) / 2);

The known point on the region 6-7 boundary is point 3. The known point on the region 6-8 boundary is point 4. Use this information in the slope-intercept form, subtracting the equations, and rearrange to yield the location of the next point.

x5 = (-x4 * slope68 + x3 * slope67 + y4 -y3) / (slope67 - slope68);
y5 = (x5 - x4) * slope68 + y4;

The second type II characteristic propagates beyond the fan at an angle averaged between the region 7 and region 9 angles for type II waves.

slope79 = tand( (typeTwo(7) + typeTwo(9)) / 2);

The last point of interest is the x-intercept of the flattest type I wave. Calculate this point by knowing the location of point 5 and finding the slope of the type I wave between region 8 and region 9.

slope89 = tand( (typeOne(8) + typeOne(9)) / 2);
y6      = 0;
x6      = (slope89 * x5 - y5) / slope89;

The final type II wave propagates away at an angle averaged between region 9 and region 10.

slope910 = tand( (typeTwo(9) + typeTwo(10)) /2);

With all the points calculated and the slope of the freely propagating lines known, connect the dots:

points = plotExpansionSchematic('nozzlepoints');

Figure contains an axes object. The axes object contains 39 objects of type line, text. One or more of the lines displays its values using only markers These objects represent Point 0 (0,1), Point 1 (1.3113,0), Point 2 (1.5216,0.16263), Point 3 (1.7322,0.33233), Point 4 (1.8036,0), Point 5 (2.0959,0.17997), Point 6 (2.4862,0).

The method of characteristics is a potent method for solving supersonic gas dynamics problems. Note, that this method represents an approximation for the characteristic lines. The approximation approaches the exact case for an infinite number of characteristic lines.

Reference

[1] James, J. E. A., "Gas Dynamics, Second Edition", Allyn and Bacon, Inc, Boston, 1984.