Iterative Numerical Solution for Mach, P, and T distribution along the

SSME Nozzle

Source https://smallsat.wordpress.com

clear all;
clc;
close all;
Pe = 16.8727;       %kPa
P0 = 20.4*1000;     %kPa
T0 = 3500;          %K

%Gas Properties
MW = 22;
gamma = 1.22;
%X Axis
X = [0.00 4.00 8.00 12.00 16.00 20.00 25.00 30.00 50.00 70.00 90.00 100.00...
120.00 140.00 160.00 180.00 200.00 220.00 240.00 260.00 280.00 300.00 305.00];
%Z Axis
Z = [17.50 15.50 13.00 12.25 13.00 15.50 18.50 22.00 33.00 41.50 50.50 54.50...
61.00 68.00 74.50 80.50 86.00 91.00 97.00 101.00 105.00 107.50 107.85];
%Diameter
Dia = 2*Z;

% Compute Mach Number at Exit
Me = (2/(gamma - 1)*(((P0/Pe)^((gamma - 1)/gamma)) - 1))^0.5;

fprintf('Mach Number at Exit = %8.6f \n',Me);
%Compute ratio of geometric throat area to Sonic throat area
AeAt = 77.5;                %SSME Nozzle geometry;
AtAs = (1/AeAt)*(1/Me*(2/(gamma + 1)*(1 + (gamma - 1)/2*Me^2))...
    ^((gamma + 1)/(2*(gamma - 1))));
fprintf('At/A* = %8.6f \n',AtAs);
% At/A* = 1, Nozzle is Choked

%Using  Iterative Solver at each Nozzle station
%Computing A/A* for each section
D_exit = 215.7;         %Exit Diameter
A_As = 77.5*(Dia/D_exit).^2;

oldMach = 4;
eps = 0.001;

for i =1:23
    Mx(i) = Solver( oldMach,A_As(i),eps,gamma);
 end
plot(X,Mx,'b');
hold on;
plot(X,Mx,'*r');
title('Mach Number Distribution in Nozzle');
xlabel('X axis [cm]');
ylabel('Mach Number');
grid on;

%Solving for Pressure Distribution in Nozzle

Px = P0./(( 1 + ((gamma - 1)/2)*Mx.*Mx).^(gamma/(gamma - 1)));
figure(2);
plot(X,Px,'b');
hold on;
plot(X,Px,'*r');
title('Pressure Distribution in Nozzle');
xlabel('X axis [cm]');
ylabel('Pressure[kPa]');
grid on;

%Pressure Distribution in Nozzle(normalized)
figure(3);
Pnx = 1./( 1 + ((gamma - 1)/2)*Mx.^2).^(gamma/(gamma - 1));
plot(X,Pnx,'b');
hold on;
plot(X,Pnx,'*r');
title('Pressure Distribution in Nozzle (normalized)');
xlabel('X axis [cm]');
ylabel('P/P_0');
grid on;

%Temperature Distribution in Nozzle
figure(4);
Tx = T0./( 1 + ((gamma - 1)/2)*Mx.^2);
plot(X,Tx,'b');
hold on;
plot(X,Tx,'*r');
title('Temperature Distribution in Nozzle');
xlabel('X axis [cm]');
ylabel('T  [ K ]');
grid on;

%Computing Exit Temperature, velocity based on Exit mach
T_exit = T0/( 1 + ((gamma - 1)/2)*Me^2);
fprintf('T_exit = %8.2f K\n',T_exit);

V_exit = (T_exit*gamma*8314.4126/MW)^0.5*Me;
fprintf('V_exit = %8.2f m/s \n',V_exit);

Mach Number at Exit = 4.858225 
At/A* = 1.000152 
T_exit =   973.23 K
V_exit =  3254.40 m/s

function Mach = Solver( oldMach,A_As,eps,gamma)
maxit = 100000;
for i= 1: maxit
newMach = oldMach -(FOFM(oldMach,gamma)- A_As)/DFDM(oldMach,gamma);
oldMach = newMach;
fofm = FOFM(oldMach,gamma);

if( abs(fofm – A_As)/A_As <= eps)
%display(‘break’);
break;
end
end
Mach = oldMach;
end

function FOFM = FOFM(Mach,gamma)

FOFM = 1./Mach.*(2/(gamma + 1)*(1 + (gamma – 1)/2.*Mach.^2)).^…
((gamma + 1)/(2*(gamma – 1)));
end

function DFDM = DFDM(Mach,gamma)

DFDM = (2^((1 – 3*gamma)/(2 – 2*gamma)))*(Mach.^2 – 1)./…
(Mach.^2.*( 2 + Mach.^2*(gamma – 1))).*…
( 1 + (gamma – 1)/2*Mach.^2/(gamma + 1)).^((gamma + 1)/(2*(gamma -1)));
end

Published with MATLAB® 7.10