How to plot Mach Number Contours
Greetings everyone, below is the code which generates the diverging portion of a supersonic nozzle using the Method Of Characteristics. How can I plot the Mach Number contours within the same plot along with a colorbar? Thank you!
%% Initialize datapoint matrices
clear
G=1.4;
Me = 2;
n = 55;
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;
%% Define some flow parameters of originating characteristic lines
dT = ThetaMax/n;
ThetaArc(:,1) = (0:dT:ThetaMax);
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.0;
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);
end
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 % do not 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
xlabel(‘Length [x/y0]’)
ylabel(‘Height [y/y0]’)
endGreetings everyone, below is the code which generates the diverging portion of a supersonic nozzle using the Method Of Characteristics. How can I plot the Mach Number contours within the same plot along with a colorbar? Thank you!
%% Initialize datapoint matrices
clear
G=1.4;
Me = 2;
n = 55;
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;
%% Define some flow parameters of originating characteristic lines
dT = ThetaMax/n;
ThetaArc(:,1) = (0:dT:ThetaMax);
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.0;
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);
end
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 % do not 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
xlabel(‘Length [x/y0]’)
ylabel(‘Height [y/y0]’)
end Greetings everyone, below is the code which generates the diverging portion of a supersonic nozzle using the Method Of Characteristics. How can I plot the Mach Number contours within the same plot along with a colorbar? Thank you!
%% Initialize datapoint matrices
clear
G=1.4;
Me = 2;
n = 55;
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;
%% Define some flow parameters of originating characteristic lines
dT = ThetaMax/n;
ThetaArc(:,1) = (0:dT:ThetaMax);
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.0;
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);
end
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 % do not 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
xlabel(‘Length [x/y0]’)
ylabel(‘Height [y/y0]’)
end #contour plot, nozzle, mach number MATLAB Answers — New Questions
