Generating a third order polynomial for the contraction region of supersonic nozzle

Greetings everyone, I have generated a third order polynomial to make the covergent section wall of a supersonic nozzle. But the result is not what I expected:

My result is attached below:(Poly3.PNG)

What I want:(Desired.jpeg)

The code which I have made:

Also, is there a way where I can ensure continuity throughout from converging section, at the throat and at the diverging section of nozzle?

Thank you!

%%function [xconv, yconv] = Convergent_3rd(G,y0,H_in,L_e,n)

%% Constants

%clear

G = 1.4;

%% Define Flow Properties

R = 287;

T0 = 300;

%% Define the mach number in contraction region

uu = 12.5; %velocity in contraction region in m/s

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;

%% Generating the polynomial

% Polynomial coefficients (a4, a3, a2, a1, a0)

syms x a0 a1 a2 a3

Y = a0 + a1*x + a2*x^2 + a3*x^3; % third order polynomial

% Define the boundary conditions

cond1 = subs(Y, x, L_e) == y0;

dy_dx = diff(Y, x);

cond2 = subs(dy_dx, x, L_e) == 0;

cond3 = subs(dy_dx, x, 0) == 0;

%d2y_dx2 = diff(dy_dx, x);

cond4 = subs(Y,x,0) == H_in;

% Solve for the coefficients

coeffs = solve([cond1, cond2, cond3, cond4], [a0, a1, a2, a3]);

a0 = coeffs.a0;

a1 = coeffs.a1;

a2 = coeffs.a2;

a3 = coeffs.a3;

fprintf('Polynomial coefficients:\n');

fprintf('a3 = %f\n', double(a3));

fprintf('a2 = %f\n', double(a2));

fprintf('a1 = %f\n', double(a1));

fprintf('a0 = %f\n', double(a0));

% Define the polynomial function using the obtained coefficients

Y_poly = matlabFunction(coeffs.a3*x.^3 + coeffs.a2*x.^2 + coeffs.a1*x + coeffs.a0);

% Plotting the polynomial

x_values = linspace(0, L_e, 55+1); % n+1

y_values = Y_poly(x_values);

plot(x_values,y_values)

xconv = x_values';

yconv = y_values';

%%end

Answer 1

[Shishir Reddy]

Hi Bamelari

As per my understanding, your code produces a plot that starts converging from “x = 0” which is the y - axis and you would like to adjust the polynomial such that the convergent section starts at “x = -L_e” and ends at “x = 0”. For this, the boundary conditions need to be adjusted accordingly.

The first step for this would be domain shift. i.e., “x” should be replaced by “x = x + L_e” in the polynomial expression:

Y_poly = matlabFunction(coeffs.a3*x.^3 + coeffs.a2*x.^2 + coeffs.a1*x + coeffs.a0);

Then, the x – values should be defined from -1.5 to 0 so that the final plot reflects this range.

The resulting code after making these changes will be as follows –

%%function [xconv, yconv] = Convergent_3rd(G,y0,H_in,L_e,n)

%% Constants

%clear

G = 1.4;

%% Define Flow Properties

R = 287;

T0 = 300;

%% Define the mach number in contraction region

uu = 12.5; %velocity in contraction region in m/s

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;

%% Generating the polynomial

% Polynomial coefficients (a4, a3, a2, a1, a0)

syms x a0 a1 a2 a3 

Y = a0 + a1*x + a2*x^2 + a3*x^3; % third order polynomial

% Define the boundary conditions

cond1 = subs(Y, x, L_e) == y0;

dy_dx = diff(Y, x);

cond2 = subs(dy_dx, x, L_e) == 0;

cond3 = subs(dy_dx, x, 0) == 0;

%d2y_dx2 = diff(dy_dx, x);

cond4 = subs(Y,x,0) == H_in;

% Solve for the coefficients

coeffs = solve([cond1, cond2, cond3, cond4], [a0, a1, a2, a3]);

a0 = coeffs.a0;

a1 = coeffs.a1;

a2 = coeffs.a2;

a3 = coeffs.a3;

fprintf('Polynomial coefficients:\n');

Polynomial coefficients:

fprintf('a3 = %f\n', double(a3));

a3 = 8.938651

fprintf('a2 = %f\n', double(a2));

a2 = -20.111965

fprintf('a1 = %f\n', double(a1));

a1 = 0.000000

fprintf('a0 = %f\n', double(a0));

a0 = 16.083974

% Define the polynomial function using the obtained coefficients

Y_poly = matlabFunction(coeffs.a3*(x+ L_e).^3 + coeffs.a2*(x+L_e).^2 + coeffs.a1*(x+L_e) + coeffs.a0);

% Plotting the polynomial

x_values = linspace(-L_e, 0, 55+1); % n+1

y_values = Y_poly(x_values);

plot(x_values,y_values

To ensure continuity, the polynomial transition should be smooth at the throat, this means that the first derivative(slope) should match at the junction. That means, throat should be specifically defined where the convergent and divergent sections meet.

I hope this helps.

[Alan Stevens]

I'd just do it like this (as I don't have the Aerospace Toolbox I've just taken your indicated value of H_in):

y0 = 1; H_in = 12.5; L_e = 1.5;

y = @(a2,a3,x) H_in + a2*x.^2 + a3*x.^3;

% dydx = 2*a2*x + 3*a3*x^2
% Clearly, y(0) = H_in from cond4
% and cond3 means a1 = 0
% y0 = Hin + a2*L_e^2 + a3*L_e^3  from cond1
%  0 = 2*a2*L_e + 3*a3*L_e^2        from cond2

% Find a2 and a3 as follows:
M = [L_e^2 L_e^3; 2*L_e 3*L_e^2];
K = [y0 - H_in; 0];

a = M\K;  a1 = a(1); a2 = a(2);
disp(a)

  -15.3333
    6.8148

x = linspace(0,L_e);
plot(x-L_e, y(a1,a2,x)),grid
xlabel('x')
ylabel('y')