close all clear clc % Gas properties and throat conditions gamma = 5/3; % cp/cv P0 = 3.0189E7; % [Pa] rho0 = 105.54; % [kg/m^3] % Compute throat velocity c0 = sqrt(gamma*P0/rho0); % [m/s] speed of sound u0 = c0; % [m/s] sonic % Load nozzle data dd = load('expanding_nozzle_geometry.dat'); x_vect = dd(:,1); y_vect = dd(:,2); A_vect = pi*y_vect.^2; A_crit = A_vect(1); % Create Mach-area relation M_test = linspace(1, 200, 10000); % max Mach: 200 A_test = A_crit./M_test.*( 2*(1+(gamma-1)/2*M_test.^2) / (gamma+1) ).^((gamma+1)/(2*(gamma-1))); % Find Mach number in the loaded nozzle M_vect = interp1(A_test, M_test, A_vect); % Find velocity u_vect = u0 * M_vect.*sqrt((gamma+1)/2./(1+(gamma-1)/2*M_vect.^2)); % Find pressure P_vect = P0 * ((gamma+1)/2./(1+(gamma-1)/2*M_vect.^2)).^(gamma/(gamma-1)); % Find density rho_vect = rho0 * (P_vect/P0).^(1/gamma); % PLOT figure plot(x_vect, M_vect, 'b', 'linewidth', 2) xlabel('x [m]', 'fontsize', 14) ylabel('Mach number', 'fontsize', 14) figure semilogy(x_vect, u_vect, 'k', 'linewidth', 2) xlabel('x [m]', 'fontsize', 14) ylabel('velocity', 'fontsize', 14) figure semilogy(x_vect, P_vect, 'k', 'linewidth', 2) xlabel('x [m]', 'fontsize', 14) ylabel('pressure', 'fontsize', 14) figure semilogy(x_vect, rho_vect, 'k', 'linewidth', 2) xlabel('x [m]', 'fontsize', 14) ylabel('density', 'fontsize', 14)