Supersonic Nozzle Design using Method of Characteristics

Greetings everyone, I am working on the design of supersonic nozzle using the method of characteristics from this paper Supersonic Nozzle Design | J. Appl. Mech. | ASME Digital Collection. The code for the minimum length expansion has been generated and plotted(Plot.PNG). However, I want to implement the "Less than Maximum Expansion". As you can see from the matlab code, the maximum expansion angle has been taken to be one half the prandtl meyer angle. Under this condition, the flow has fully expanded to reach the exit mach number and does not need to go further expansions or reflections from the wall(Plot.PNG). The waves gets cancelled.

However, if we choose the expansion angle to be smaller than the maximum expansion angle, the flow would not reach the exit mach number since it has not expanded fully. Thus it needs to undergo reflections from the wall until it reaches the desired exit mach number and the reflected waves get cancelled after certain reflections.

The cancellation(no reflection) criteria is that the angle of the wall point and the point just behind should be the same

The problem, I am unable to understand and implement the above "Less than expansion method" from the provided paper. How exactly would the waves reflect and get cancelled? I want to obtain the Desired.PNG plot. Any help in this matter is highly appreciated. Screenshots from the paper has also been attached!

Thank you!

%% Initialize datapoint matrices

clear

G=1.4;

Me = 2.0;

n = 4; %speed index/no of characteristics. As per paper it is 53. For visualization i took 4

display = 1;

Km = zeros(n,n); % K- vlaues (Constant along right running characteristic lines)

Kp = zeros(n,n); % K+ vlaues (Constant along left running characteristic lines)

Theta = zeros(n,n); % Flow angles relative to the horizontal

Mu = zeros(n,n); % Mach angles

M = zeros(n,n); % Mach Numbers

x = zeros(n,n); % x-coordinates

y = zeros(n,n); % y-coordinates

%% Generate the contraction region using a 3rd order polynomial. Do not run this

%{

R = 287;

T0 = 300;

% Define the mach number in contraction region

uu = 25.0; %velocity in contraction region

a = sqrt(G*R*T0 - 0.5*(G-1)*uu^2); % local speed of sound

Mconv = uu/a; % Mach number

% Calculate the area ratio of contraction region

[~,~,~,~,area] = flowisentropic(G,Mconv); % area ratio

% Inlet convergent section geometry

y0 =1;

H_in = area * y0;

L_e = 1.5; %(1.0/3.0)*xwall(end);

%

[xconv,yconv] = Convergent_3rd(G,y0,H_in,L_e,n);

%}

%% Generate the convergent portion of a nozzle

% The inlet height/area and exit height/area(divergent) are same

% Therefore we have same A/A* = 1.6875

% Corresponding to Mach = 0.3722

%% Find NuMax (maximum expansion angle)

[~, NuMax, ~] = PMF(G,Me,0,0);

ThetaMax = NuMax/2;

dT = ThetaMax/n;

%% from pucketts paper provided: this section

%theta0 = (1.0/1.687)^(2/9) * (NuMax/2);

%no_ref = (ThetaMax - theta0)/dT; % no of reflections

%% Define some flow parameters of originating characteristic lines

ThetaArc(:,1) = (0:dT:ThetaMax); %ThetaArc(:,1) = (0:dT:theta0) as per puckett

NuArc = ThetaArc;

KmArc = ThetaArc + NuArc;

[~, ~, MuArc(:,1)] = PMF(G,0,NuArc(:,1),0);

%% Coordinates of wall along curve from throat

y0 = 1; % Define throat half-height

ThroatCurveRadius = 1.5*y0; % Radius of curvature just downstream of the throat

% for larger factors, ywall deviates from A/A* preferred value is 1.1

%L_e = 1.1 * y0 * sind(ThetaMax);

[xarc, yarc] = Arc(ThroatCurveRadius,ThetaArc); % Finds x- and y-coordinates for given theta-values

yarc(:,1) = yarc(:,1) + y0; % Defines offset due to arc being above horizontal

%% Fill in missing datapoint info along first C+ line

% First centerline datapoint done manually

Km(:,1) = KmArc(2:length(KmArc),1);

Theta(:,1) = ThetaArc(2:length(KmArc),1);

Nu(:,1) = Theta(:,1);

Kp(:,1) = Theta(:,1)-Nu(:,1);

M(1,1) = 1.0001;

Nu(1,1) = 0;

Mu(1,1) = 90;

y(1,1) = 0;

x(1,1) = xarc(2,1) + (y(1,1) - yarc(2,1))/tand((ThetaArc(2,1) - MuArc(2,1) - MuArc(2,1))/2);

% Finds the information at interior nodes along first C+ line

for i=2:n

[M(i,1), Nu(i,1), Mu(i,1)] = PMF(G,0,Nu(i,1),0);

s1 = tand((ThetaArc(i+1,1) - MuArc(i+1,1) + Theta(i,1) - Mu(i,1))/2);

s2 = tand((Theta(i-1,1) + Mu(i-1,1) + Theta(i,1) + Mu(i,1))/2);

x(i,1) = ((y(i-1,1)-x(i-1,1)*s2)-(yarc(i+1,1)-xarc(i+1,1)*s1))/(s1-s2);

y(i,1) = y(i-1,1) + (x(i,1)-x(i-1,1))*s2;

end

%% Find flow properties at remaining interior nodes

for j=2:n;

for i=1:n+1-j;

Km(i,j) = Km(i+1,j-1);

if i==1;

Theta(i,j) = 0;

Kp(i,j) = -Km(i,j);

Nu(i,j) = Km(i,j);

[M(i,j), Nu(i,j), Mu(i,j)] = PMF(G,0,Nu(i,j),0);

s1 = tand((Theta(i+1,j-1)-Mu(i+1,j-1)+Theta(i,j)-Mu(i,j))/2);

x(i,j) = x(i+1,j-1) - y(i+1,j-1)/s1;

y(i,j) = 0;

else

Kp(i,j) = Kp(i-1,j);

Theta(i,j) = (Km(i,j)+Kp(i,j))/2;

Nu(i,j) = (Km(i,j)-Kp(i,j))/2;

[M(i,j), Nu(i,j), Mu(i,j)] = PMF(G,0,Nu(i,j),0);

s1 = tand((Theta(i+1,j-1)-Mu(i+1,j-1)+Theta(i,j)-Mu(i,j))/2);

s2 = tand((Theta(i-1,j)+Mu(i-1,j)+Theta(i,j)+Mu(i,j))/2);

x(i,j) = ((y(i-1,j)-x(i-1,j)*s2)-(y(i+1,j-1)-x(i+1,j-1)*s1))/(s1-s2);

y(i,j) = y(i-1,j) + (x(i,j)-x(i-1,j))*s2;

end

end

end

%% Find wall node information

xwall = zeros(2*n,1);

ywall = xwall;

ThetaWall = ywall;

xwall(1:n,1) = xarc(2:length(xarc),1);

ywall(1:n,1) = yarc(2:length(xarc),1);

ThetaWall(1:n,1) = ThetaArc(2:length(xarc),1);

for i=1:n-1

ThetaWall(n+i,1) = ThetaWall(n-i,1); % criteria for stopping the reflection from the wall

end

%% Location of wall points

for i=1:n

s1 = tand((ThetaWall(n+i-1,1) + ThetaWall(n+i,1))/2);

s2 = tand(Theta(n+1-i,i)+Mu(n+1-i,i));

xwall(n+i,1) = ((y(n+1-i,i)-x(n+1-i,i)*s2)-(ywall(n+i-1,1)-xwall(n+i-1,1)*s1))/(s1-s2);

ywall(n+i,1) = ywall(n+i-1,1) + (xwall(n+i,1)-xwall(n+i-1,1))*s1;

end

%% Provide wall geometry to user

assignin('caller','xwall',xwall)

assignin('caller','ywall',ywall)

assignin('caller','Coords',[xwall ywall])

%% Generate the convergent portion of nozzle % Dont run this

H_in = ywall(end);

L_e = (xwall(end)*(1.0/3.0));

[xconv,yconv] = Convergent_new_3rd(y0,H_in,L_e,n);

%%

% Draw contour and characteristic web

if display == 1

plot(xwall,ywall,'-')

axis equal

axis([0 ceil(xwall(length(xwall),1)) 0 ceil(ywall(length(ywall),1))])

hold on

plot(xarc,yarc,'k-')

for i=1:n-1

plot(x(1:n+1-i,i),y(1:n+1-i,i))

end

for i=1:n

plot([xarc(i,1) x(i,1)],[yarc(i,1) y(i,1)])

plot([x(n+1-i,i) xwall(i+n,1)],[y(n+1-i,i) ywall(i+n,1)])

end

for c=1:n

for r=2:n+1-c

plot([x(c,r) x(c+1,r-1)],[y(c,r) y(c+1,r-1)])

end

end

%hold on

%contourf(x, y, M, 20, 'LineColor', 'none'); % Draw Mach number contours

%colorbar;

xlabel('Length [x/y0]')

ylabel('Height [y/y0]')

end

Answer 1

[Shishir Reddy]

Hi Bamelari,

To implement the "Less than Maximum Expansion" method in the supersonic nozzle design, your approach needs to be adjusted to account for the reflections and cancellations of expansion waves when the expansion angle is less than the maximum.

In a nozzle designed for less than maximum expansion, the initial expansion angle is set lower than what would be required for direct expansion to the desired Mach number. This means the flow needs additional adjustments, achieved through reflections.

Each reflection of the expansion waves further adjusts the flow properties, incrementally increasing the Mach number. The reflections act as additional expansion processes that help the flow reach the desired conditions.

Implementation steps:

1. Adjust Initial Expansion Angle: Modify the initial expansion angle to be less than the maximum. You can set it to a fraction of ThetaMax.

theta0 = (1.0/1.687)^(2/9) * (NuMax/2);
dT = theta0 / n; 

2. Calculating Number of Reflections:

no_ref = ceil((ThetaMax - theta0) / dT);

3. Validation: Consider adding code for validation (e.g. checking if the exit Mach number is achieved)

if abs(exitMach - desiredExitMach) < 0.01
    fprintf('Validation successful: The exit Mach number is within acceptable limits.\n');
else
    fprintf('Validation failed: The exit Mach number is not within acceptable limits.\n');

By adjusting the code according to these steps, the "Less than Maximum Expansion" method can be implemented for the supersonic nozzle design.

I hope this helps.