%% Solving for the Exit Flow of a Supersonic Nozzle
% 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.

% Copyright 2009-2012 The MathWorks, Inc.

%% 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 = astexpandschematic('uppernozzle');

%%
% 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.  
%
% In the previous figure, the blue lines in the 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.  
%
% $$ \theta_n = Flow\;angle\;in\;the\;n^{th}\;region $$
%

%%
% * The Prandtl-Meyer angle is the angle at which the flow changes direction from one region to another.
%
% $$ \nu_n = Prandtl-Meyer\;angle\;in\;the\;n^{th}\;region $$
%

%%
% * Mach angle is the angle between the local flow direction and the weak pressure waves that emanate from a given point.
%
% $$ \mu_n = Mach\;angle\;in\;the\;n^{th}\;region $$
%
%%
% 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 locate 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).
%
% $$ \theta_4 = \nu_4 - \nu_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, the
% calculate the change in flow angle in crossing type I characteristics
% from region 1 to region 4:
%
% $$ \Delta \theta_I = \frac{\theta_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 nor 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 is
% equal to the change in flow angle:
%
% $$ \Delta \nu_I = \Delta \theta_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.
%
% $$ \Delta \nu_I = \nu_2 - \nu_1 $$
%
%
% $$ \nu_2 = \nu_1 + \Delta \nu_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
%
% $$ \Delta \theta_{II} = \theta_5 - \theta_2 $$
%

deltaThetaII = theta_5 - theta_2;

%%
% Calculate the change in Prandtl-Meyer angle across type II
% characteristics:
%
% $$ \Delta \nu_{II} = - \Delta \theta_{II} $$
%

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,
%
% $$ \Delta \theta_I = \Delta \nu_I = \theta_6 - \theta_5 = \nu_6 - \nu_5 $$
%

%%
% Rearranged this as:
%
% $$ \nu_6 - \theta_6 = \nu_5 - \theta_5 \;\;\;\;\;(1)$$
%

%%
% A type II characteristic is crossed in going from region 3 to region 6.
% Therefore,
%
% $$ \Delta \nu_{II} = -\Delta \theta_{II} $$
%
%
% $$ \nu_6 - \nu_3 = \theta_3 - \theta_6 $$
%

%%
% Rearranging this as:
%
% $$ \nu_6 + \theta_6 = \nu_3 + \theta_3 \;\;\;\;\;(2)$$
%

%%
% Add equations (1) and (2) together, then solve for the Prandtl-Meyer
% angle in region 6.  This yields the following expression.
%
% $$ \nu_6 = \frac{(\nu_5 - \theta_5) + (\nu_3 + \theta_3)}{2}$$
%

%%
% In MATLAB(R), 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 just like
% the way 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:
%
% $$ \left(\!\frac{dy}{dx}\!\right)_{\! I} = tan(\theta - \mu)$$
%
%
% $$ \left(\!\frac{dy}{dx}\!\right)_{\!\! II} = tan(\theta + \mu)$$
%

%%
% 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(1,10);
machAngles  = zeros(1,10);
typeOne     = zeros(1,10);
typeTwo     = zeros(1,10);

for i = 1:10
    [machNumbers(i), ~, machAngles(i)] = flowprandtlmeyer(k, prandtlMeyerAngles(i), 'nu');
    typeOne(i) = flowAngles(i) - machAngles(i);
    typeTwo(i) = flowAngles(i) + machAngles(i);
end


clear table;
table(1,:) = 'Region   theta        nu        Mach        mu       type I    type II';
table(2,:) = '         (Deg)      (Deg)                  (Deg)     (Deg)      (Deg) ';
for m=1:length(machNumbers)
table(m+3,:) = sprintf('%3.0d   %8.2f      %5.2f   %8.3f    %8.2f   %8.2f  %8.2f ', ...
                m, flowAngles(m), prandtlMeyerAngles(m), machNumbers(m), ...
                machAngles(m), typeOne(m), typeTwo(m));
end

disp(table)




%%
% 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 in order 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 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,
% separates 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
% referenced 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} = Slope\;of\;the\;boundary\;between\;region\;i\;and\;region\;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 = astexpandschematic('nozzlepoints');

%%
% 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.

%%
close(upperNozzle,points)

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

%%
displayEndOfDemoMessage(mfilename)
