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');
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+ 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.
