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? The plot which I want is attached below(plot.png) Thank you!

plot.png reference: https://doi.org/10.1063/5.0150327

I have attached the function files.

PMF = Calculates the prandtl meyer angle, mach angle and Mach Number

Arc = Generates the circular arc for the expansion region of nozzle(Region where the lines emanate)

%% 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])

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

[Cris LaPierre]

Can you provide the missing variable(s)? Please save them to a mat file and attach it to your post using the paperclip icon. We at least need PMF.

[Bamelari Jovani]

I have attached the required function files.

[Cris LaPierre]

Thank you. When I ran the code, I discovered that Convergent_new_3rd is also missing. Once you attach everything, click the green play button in the toolstrip to execute the code. If nothing is missing, your code should run without error.

[Bamelari Jovani]

The convergent_new_3rd function is not needed here in this case. You can ignore it. I only want to know how I can incorporate the Mach Number Contours on the same plot. Any help in this matter is appreciated

Answer 1

[Voss]

Adding the line

contourf(x,y,M)

to the end of your code produces the following result:

%% 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])

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

contourf(x,y,M)

Is that what's intended?

[Bamelari Jovani]

I want it to look like Plot.PNG which I have attached above

[Voss]

OK, well that plot has no numbers on either axis, so what is it?