clc; clear; close all;
1. Problem Setup
gamma = 1.4; % Specific heat ratio for air
P0 = 101325; % Stagnation (total) pressure, Pa
T0 = 300; % Stagnation (total) temperature, K
R = 287; % Gas constant for air, J/kg-K
x_vals = 0:0.05:1; % x from 0 to 1 with step of 0.05
% Given nozzle contour:
A = @(x) 11.9774*x.^2 - 15.4747*x + 6; % area distribution
% Evaluate nozzle area over x
A_of_x = A(x_vals);
% Find the throat area A*
[A_throat, idx_throat] = min(A_of_x);
x_throat = x_vals(idx_throat);
% List of back-pressure ratios to study:
Pinf_over_P0_list = [1.0 0.98 0.97 0.96085 0.95 0.9 0.8 0.7 0.6 0.5 ...
0.435 0.4 0.3 0.2 0.1 0.04];
2. Plot Nozzle Contour
figure(1);
plot(x_vals, A_of_x, 'b-o','LineWidth',2);
xlabel('x'); ylabel('A(x)'); grid on;
title('Nozzle Contour: A(x) vs x');
3. Prepare arrays to store results
nCases = length(Pinf_over_P0_list);
M_inlet = zeros(1,nCases);
M_throat = zeros(1,nCases);
M_exit = zeros(1,nCases);
Pexit_over_P0 = zeros(1,nCases);
mdot_array = zeros(1,nCases);
4. Loop Over Each Back-Pressure Ratio
for i = 1:nCases
Pinf_over_P0 = Pinf_over_P0_list(i);
% ---- (A) Check if nozzle is "choked" ----
% A simplified way: compare with approximate critical ratio
P_crit_ratio = ( (gamma+1)/2 )^(-gamma/(gamma-1));
% For gamma=1.4, P_crit_ratio ~ 0.528
isChoked = (Pinf_over_P0 < P_crit_ratio);
% ---- (B) Build M(x) along the nozzle using isentropic area relation ----
M_profile = zeros(size(x_vals));
for j = 1:length(x_vals)
% Local area ratio:
A_local = A_of_x(j);
% A* known from the nozzle geometry (the minimum)
ARatio = A_local / A_throat;
% Solve the area-Mach equation for M:
% A/A* = 1/M * [ (2/(gamma+1)) (1 + (gamma-1)/2 * M^2 ) ]^[(gamma+1)/(2(gamma-1))]
% There are two solutions: subsonic M < 1 or supersonic M > 1
% We'll pick subsonic for x < x_throat, supersonic for x > x_throat
if isChoked
% If the nozzle is choked, throat is M=1, diverging section supersonic
if j <= idx_throat
% subsonic branch
M_sub = fsolve(@(M) areaMachEq(M,ARatio,gamma), 0.2, ...
optimset('Display','off'));
M_profile(j) = M_sub;
else
% supersonic branch
M_sup = fsolve(@(M) areaMachEq(M,ARatio,gamma), 2.0, ...
optimset('Display','off'));
M_profile(j) = M_sup;
end
else
% Not choked: entire nozzle subsonic
M_sub = fsolve(@(M) areaMachEq(M,ARatio,gamma), 0.2, ...
optimset('Display','off'));
M_profile(j) = M_sub;
end
end
% Store key Mach numbers
M_inlet(i) = M_profile(1);
M_throat(i) = M_profile(idx_throat);
M_exit(i) = M_profile(end);
% ---- (C) Pressure distribution along x (isentropic) ----
% P(x)/P0 = [1 + (gamma-1)/2 M^2]^(-gamma/(gamma-1))
P_over_P0_profile = (1 + (gamma-1)/2 * M_profile.^2).^(-gamma/(gamma-1));
% BUT physically, if there's back pressure P_inf, the actual exit P might get "clamped."
% In purely isentropic approach, we ignore shock. We assume flow adjusts at the exit.
% Record exit pressure ratio (isentropic approach):
Pexit_over_P0(i) = P_over_P0_profile(end);
% ---- (D) Mass Flow Rate (if choked, depends on throat) ----
if isChoked
% Choked formula at the throat:
Tt = T0 / (1 + (gamma-1)/2 * M_throat(i)^2);
at = sqrt(gamma * R * Tt);
rho0 = P0/(R*T0);
% Velocity at throat:
Vthroat = M_throat(i)*at;
A_th = A_throat; % throat area
mdot = rho0 * Vthroat * A_th; % approximate isentropic
else
% Not choked, approximate subsonic entire nozzle
% Just take final cross-section as example (not strictly correct)
T_exit = T0 / (1 + (gamma-1)/2 * M_exit(i)^2);
a_exit = sqrt(gamma * R * T_exit);
rho0 = P0/(R*T0);
Vexit = M_exit(i)*a_exit;
A_ex = A_of_x(end);
mdot = rho0 * Vexit * A_ex;
end
mdot_array(i) = mdot;
end
5. Plot M vs x for a Sample Back Pressure
Just pick a single case to illustrate
sampleIndex = find(abs(Pinf_over_P0_list - 0.7) < 1e-6, 1); % e.g., 0.7
if ~isempty(sampleIndex)
Pinf_over_P0 = Pinf_over_P0_list(sampleIndex);
% Recompute the M_profile for that case, for plotting
[M_profile, P_profile] = computeIsentropicProfile( ...
x_vals, A_of_x, idx_throat, ...
gamma, Pinf_over_P0);
figure(2);
plot(x_vals, M_profile, 'r-o','LineWidth',2);
xlabel('x'); ylabel('M(x)');
title(['Mach Number vs x for P_{\infty}/P_{0} = ',num2str(Pinf_over_P0)]);
grid on;
figure(3);
plot(x_vals, P_profile, 'b-o','LineWidth',2);
xlabel('x'); ylabel('P(x)/P_{0}');
title(['Pressure Ratio vs x for P_{\infty}/P_{0} = ',num2str(Pinf_over_P0)]);
grid on;
end
6. Plot the Summaries vs P_{\infty}/P_{0}
figure(4);
plot(Pinf_over_P0_list, mdot_array, 'ko-','LineWidth',2);
xlabel('P_{\infty}/P_{0}');
ylabel('Mass Flow Rate (kg/s)');
title('m_{dot} vs P_{\infty}/P_{0}'); grid on;
figure(5);
plot(Pinf_over_P0_list, M_inlet, 'r-o','LineWidth',2); hold on;
plot(Pinf_over_P0_list, M_throat, 'g-o','LineWidth',2);
plot(Pinf_over_P0_list, M_exit, 'b-o','LineWidth',2);
legend('M_{inlet}','M_{throat}','M_{exit}');
xlabel('P_{\infty}/P_{0}'); ylabel('Mach Number');
title('M_{inlet}, M_{throat}, M_{exit} vs P_{\infty}/P_{0}'); grid on;
figure(6);
plot(Pinf_over_P0_list, Pexit_over_P0, 'm-o','LineWidth',2);
xlabel('P_{\infty}/P_{0}');
ylabel('P_{exit}/P_{0}');
title('Exit Pressure Ratio vs P_{\infty}/P_{0}');
grid on;
--- Supporting Functions ---
function F = areaMachEq(M, ARatio, gamma)
% Isentropic area--Mach equation residual: we want F=0
% ARatio = (1/M)*[ (2/(gamma+1))*(1+((gamma-1)/2)*M^2 ) ]^((gamma+1)/(2*(gamma-1)))
term = (2/(gamma+1))*(1+(gamma-1)/2*M^2);
F = ARatio - (1/M)*term^((gamma+1)/(2*(gamma-1)));
end
function [M_profile, P_profile] = computeIsentropicProfile( ...
x_vals, A_of_x, idx_throat, gamma, Pinf_over_P0)
% Example subfunction that returns M(x) and P(x)/P0 for one chosen back pressure ratio.
% (Similar approach to the main loop.)
A_throatLocal = min(A_of_x);
M_profile = zeros(size(x_vals));
% Quick choke check:
P_crit_ratio = ((gamma+1)/2)^(-gamma/(gamma-1));
isChoked = (Pinf_over_P0 < P_crit_ratio);
for j = 1:length(x_vals)
ARatio = A_of_x(j)/A_throatLocal;
if isChoked
% subsonic for x < throat, supersonic for x > throat
if j <= idx_throat
M_sub = fsolve(@(M) areaMachEq(M,ARatio,gamma), 0.2, ...
optimset('Display','off'));
M_profile(j) = M_sub;
else
M_sup = fsolve(@(M) areaMachEq(M,ARatio,gamma), 2.0, ...
optimset('Display','off'));
M_profile(j) = M_sup;
end
else
% fully subsonic
M_sub = fsolve(@(M) areaMachEq(M,ARatio,gamma), 0.2, ...
optimset('Display','off'));
M_profile(j) = M_sub;
end
end
% Isentropic P/P0
P_profile = (1 + (gamma-1)/2 .* M_profile.^2).^(-gamma/(gamma-1));
end
