Unable to get a continuous curvature for my nozzle contour

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)

Answer 1

[Divyajyoti Nayak]

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!