Supersonic Nozzle Design using Method of Characteristics

136 views (last 30 days)
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

Answers (1)

Shishir Reddy
Shishir Reddy on 6 Sep 2024
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.

Categories

Find more on Programming in Help Center and File Exchange

Products


Release

R2024a

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!