How to plot Mach Number Contours

Question

Bamelari Jovani 2024-8-22，13:16

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2147224-how-to-plot-mach-number-contours

评论： Voss 2024-8-22，16:56

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

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

xlabel('Length [x/y0]')

ylabel('Height [y/y0]')

end

4 个评论
显示 2更早的评论隐藏 2更早的评论

Cris LaPierre 2024-8-22，13:57

编辑：Cris LaPierre 2024-8-22，13:59

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 2024-8-22，14:14

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 2024-8-22，16:21

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2147224-how-to-plot-mach-number-contours#answer_1503614

在 MATLAB Online 中打开

Arc.m
PMF.m

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

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

xlabel('Length [x/y0]')

ylabel('Height [y/y0]')

end

contourf(x,y,M)

Is that what's intended?

2 个评论
显示无隐藏无

Bamelari Jovani 2024-8-22，16:50

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

Voss 2024-8-22，16:56

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

请先登录，再进行评论。

How to plot Mach Number Contours

4 个评论
显示 2更早的评论隐藏 2更早的评论

回答（1 个）

2 个评论
显示无隐藏无

另请参阅

标签

产品

版本

Community Treasure Hunt

How to plot Mach Number Contours

4 个评论 显示 2更早的评论隐藏 2更早的评论

回答（1 个）

2 个评论 显示 无隐藏 无

另请参阅

标签

产品

版本

Community Treasure Hunt

WeChat

4 个评论
显示 2更早的评论隐藏 2更早的评论

2 个评论
显示无隐藏无