Unable to get a continuous curvature for my nozzle contour

3 次查看（过去 30 天）
Bamelari Jovani
Bamelari Jovani 2024-9-13
回答： Divyajyoti Nayak 2024-9-13
Greetings everyone, I am trying to fit a Bezier Curve as shown in the paper provide(fig1), but I am unable to get the desired plot(fig 2 and fig3). I have attached the codes as well as the paper and plot. Any help on this is highly appreciated!.
PS: The inflection point is the last point on xarc, i.e, xarc(end) and Nozzle_WithCircARc.m is the main file. It would be better to run the file cell by cell.
%function Nozzle_WithCircARc(G,Me,n,display)
%% Initialize datapoint matrices
clearvars;close all;clc;
G=1.4;
Me = 2.0;
n = 53; % speed index from pucketts paper
display = 0;
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 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;
%%
%ThetaMax0 = (1.0/1.687)^(2/9) * (NuMax/2);
%% Define some flow parameters of originating characteristic lines
dT = ThetaMax/n;
%no_ref = (ThetaMax - ThetaMax_ref)/dT;
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.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
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
%% Non-scaled/non-dimensionalized plot
figure (1)
plot(xwall,ywall,'.b')
title("Non dimensionalized plot")
axis equal
hold on
plot(xarc,yarc,'-')
hold on
plot(xconv,yconv,'-')
xlabel('xwall')
ylabel('ywall')
axis equal
%%
%{
plot(xwall,ThetaWall)
xlabel('xwall')
ylabel('ThetaWall')
%%
x1 = linspace(min(xwall),max(xwall),1000);
%y1 = spline(xw,yw,x1);
y1 = interp1(xwall,ywall,x1,'spline');
figure;
title('Subplot 2: Cubic spline interpolation')
plot(x1,y1,'r')
%}
%%
%{
for i = 1:length(M)
M_centerline(i) = M(1,i);
end
for i =1:length(x)
x_axial(i) = x(1,i);
end
plot(x_axial', Mcenterline','r','LineWidth',2)
%}
%% Find the scaling factor
Mexit = Me;
[~,~,~,~,area] = flowisentropic(G,Mexit);
Hexit = 11.2798; %mm
Full_Throat_height = Hexit/area;
half_yt = Full_Throat_height/2.0;
%% Scaling of factors
xarc = half_yt.*xarc;
yarc = half_yt.*yarc;
xwall = half_yt.*xwall;
ywall = half_yt.*ywall;
xconv = half_yt.*xconv;
yconv = half_yt.*yconv;
%% Scaled up coordinates in mm
%{
title("Dimensionalized Plot")
plot(xwall,ywall,'-')
hold on
plot(xarc,yarc,'-')
hold on
plot(xconv,yconv,'-')
xlabel('xwall')
ylabel('ywall')
axis equal
%}
%% Combine Coordinates
% Combine xconv and xwall
coords_new_x = [xconv; xwall];
coords_new_y = [yconv; ywall];
%%
figure (2)
hold on;
plot(xconv, yconv, '.b', 'LineWidth', 1.5) % Convergent section
plot(xwall, ywall, '.r', 'LineWidth', 1.5) % Wall section
plot(xarc,yarc,'.g','LineWidth',1.5) % arc
plot(xconv, -1.*yconv, '.b', 'LineWidth', 1.5) % Convergent section
plot(xwall, -1.*ywall, '.r', 'LineWidth', 1.5) % Wall section
plot(xarc,-1.*yarc,'.g','LineWidth',1.5) % arc
xlabel('x [mm]');
ylabel('y [mm]');
title('Scaled up CD nozzle');
legend('Convergent Section', 'Straightening Section', 'Initial Expansion (Arc)');
grid on;
axis equal
%% Export coordinates
x1 = linspace(min(coords_new_x ),max(coords_new_x),1000);
%y1 = spline(xw,yw,x1);
y1 = interp1(coords_new_x,coords_new_y,x1,'spline');
x2 = x1';
y2 = y1';
%figure;
%plot(x1,y1,'r')
%writematrix(x2,'Spline.xlsx','Sheet',1,'Range','A1');
%writematrix(y2,'Spline.xlsx','Sheet',1,'Range','B1');
%%
%{
plot(x2,y2)
%%
plot(xwall,ywall,'-')
axis equal
axis([0 ceil(xwall(length(xwall),1)) 0 ceil(ywall(length(ywall),1))])
hold on
plot(xarc,yarc,'k-')
%}
%% Coordinates of wall
%{
writematrix(coords_new_x,'coordinates.xlsx','Sheet',1,'Range','A1');
writematrix(coords_new_y,'coordinates.xlsx','Sheet',1,'Range','B1');
%}
%%
%{
plot(x2,y2,'-k')
hold on
plot(x2,-1.*(y2),'-r')
grid on
title("CD Nozzle using a Circular Arc")
axis equal
%}
%% Analyze the discontinuity
% Find second-order derivative of wall contour
dy2_dx2 = zeros(size(xwall));
% Loop over all points except the last two
for i = 1:length(xwall)-2
% Calculate the second-order derivative using the forward difference formula
dy2_dx2(i) = (ywall(i+2) - 2*ywall(i+1) + ywall(i)) / (xwall(i+1) - xwall(i))^2;
end
dy2_dx2(end-1:end) = NaN; % last two points
%% Bezier Curve
threshold = 0.05; % Define a threshold for discontinuity detection (can be adjusted)
discontinuities = find(abs(diff(dy2_dx2)) > threshold);
%%
if ~isempty(discontinuities)
for idx = 1:length(discontinuities)
point_idx = discontinuities(idx); % Index of discontinuity
% Skip a few points around the discontinuity
region_start = max(1, point_idx - 2); % Upstream region
region_end = min(length(xwall), point_idx + 2); % Downstream region
% Use Bézier curve fitting only in this region to smooth the discontinuity
P0 = [xwall(region_start), ywall(region_start)]; % Start point
P1 = [xwall(point_idx), ywall(point_idx)]; % Control point (inflection/discontinuity)
P2 = [xwall(region_end), ywall(region_end)]; % End point
% Parameter t varies between 0 and 1 to interpolate the Bézier curve
t = linspace(0, 1, region_end - region_start + 1)';
xwall(region_start:region_end) = (1-t).^2 * P0(1) + 2*(1-t).*t * P1(1) + t.^2 * P2(1);
ywall(region_start:region_end) = (1-t).^2 * P0(2) + 2*(1-t).*t * P1(2) + t.^2 * P2(2);
end
end
%%
plot(xwall,ywall,'.b')
axis equal
%% Checking for continuous second order derivative of wall contour
cont_dy2_dx2 = zeros(size(xwall));
% Loop over all points except the last two (since forward difference needs i+2)
for i = 1:length(xwall)-2
% Calculate the second-order derivative using the forward difference formula
cont_dy2_dx2(i) = (ywall(i+2) - 2*ywall(i+1) + ywall(i)) / (xwall(i+1) - xwall(i))^2;
end
cont_dy2_dx2(end-1:end) = NaN;
plot(xwall./xwall(end), dy2_dx2)

请先登录，再回答此问题。

回答（1 个）

Divyajyoti Nayak
Divyajyoti Nayak 2024-9-13
Hi @Bamelari Jovani
From what I can understand, you are not able to plot the second order derivative curves. I think this is because the forward difference formula in the code cannot be used directly in this case. The formula used in the code is an approximation when the difference in x values (delta x) are small and can be considered equal, but in this case the values of "xwall" are not equally spaced. Initially the difference is small but after the inflection point it is huge, hence the sharp rise in "dy2_dx2".
To find the second order derivative using difference formulation, it would be better to find the first derivative using difference formula and then differentiate that to find the second derivative. Here's a sample code for that:
dy_dx = zeros(size(xwall));
for i = 1:length(xwall)-1
dy_dx(i) = (ywall(i+1) - ywall(i))/(xwall(i+1)-xwall(i));
end
dy_dx(end) = NaN;
for i = 1:length(xwall)-2
dy2_dx2(i) = (dy_dx(i+1)-dy_dx(i))/(xwall(i+1)-xwall(i));
end
dy2_dx2(end-1:end) = NaN; % last two points
This gives the correct derivatives but it still doesn't match the graphs in the given pictures. That means there's something wrong in the values of "xwall" or "ywall". If you could provide the whole research paper or the section which is used to construct the wall data, that would be helpful.
Hope this helps!

标签

产品


版本

R2024a

Community Treasure Hunt

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

Start Hunting!

Translated by