Week 7 - Simulation of a 1D Supersonic nozzle flow simulation using Macormack Method
Introduction to 1D Supersonic Nozzle Flow Simulation:
Supersonic nozzle flow simulation is a computational study aimed at understanding the complex behaviour of high speed compressible gas flow through a converging -diverging nozzle. This type of simulation is crucial in aerospace engineering,populsion systems, and fluid dynamics research.
In a 1D ( One-dimensional) supersonic flow simulation, the focus is on the flow behaviour along the axial direction of the nozzle. This simplification allows researchers to model and analyze the flow characteristics without the computational complexity associated with three dimensional simulations. The simulation typically involves solving a set of governing equations, such as Euler or Navier Stokes equations, which describe the conservation of mass, momentum and energy
Key aspects of 1D supersonic nozzle flow simulation include:
1. Converging-Diverging Nozzle Geometry:
- The nozzle features a carefully designed shape that includes a converging section, a throat, and a diverging section.
- The converging section accelerates the flow, leading to supersonic velocities at the throat.
- The diverging section allows for expansion, optimizing the flow for efficient supersonic propulsion.
2. Compressible Flow:
- The simulation considers the compressibility of the gas, which becomes significant at high velocities.
- Shock waves, expansion waves, and other compressible flow phenomena are expected and need to be accurately captured.
3. Shock Waves and Throat Conditions:
- Supersonic flows often involve shock waves, discontinuities that occur when the flow transitions from subsonic to supersonic or vice versa.
- The throat, where the cross-sectional area is smallest, is a critical point where the flow reaches its maximum velocity and conditions are most challenging.
4. Numerical Methods:
- Computational fluid dynamics (CFD) techniques, such as finite volume or finite difference methods, are commonly employed for 1D supersonic nozzle flow simulations.
- Time-marching algorithms and predictor-corrector methods are often used to solve the governing equations over discrete time steps.
5. Applications:
- Understanding supersonic nozzle flow is vital in the design and optimization of propulsion systems, including rocket engines and jet nozzles.
- These simulations contribute to the development of efficient and high-performance aerospace vehicles.
Matlab code :
clear all
close all
clc
% Inputs
% number of grids
n = 31;
% length
L = 3;
x = linspace(0,L,n);
dx = x(2)-x(1);
% To find the analytical solution
analytical_solution = 0.579*ones(1,n);
% comparing conversative and non conservative form
[x, M_nc, rho_nc, v_nc, p_nc, T_nc, mfr_nc, rho_throat_nc,v_throat_nc, p_throat_nc, T_throat_nc, M_throat_nc, mfr_throat_nc] = non_conser();
[x, M_c, rho_c, v_c, p_c, T_c, mfr_c, rho_throat_c,v_throat_c, p_throat_c, T_throat_c, M_throat_c, mfr_throat_c] = conservative();
% Plotting non -conseravtive form
figure(3)
subplot (5,1,1)
plot(x,v_nc,'color','r','linewidth',0.5)
xlabel('Length of the nozzle','FontSize',10)
ylabel('velcoity','FontSize',10)
title ('Time wise variation of primitive variables for Non-Conservative form','FontSize',12)
grid on
subplot(5,1,2)
plot(x,rho_nc,'color','g','linewidth',0.5)
xlabel('Length of the nozzle','FontSize',10)
ylabel('Density','FontSize',10)
grid on
subplot(5,1,3)
plot(x,p_nc,'color','b','linewidth',0.5)
xlabel('Length of the nozzle','FontSize',10)
ylabel('Pressure','FontSize',10)
grid on
subplot(5,1,4)
plot(x,T_nc,'color','c','linewidth',0.5)
xlabel('Length of the nozzle','FontSize',10)
ylabel('Temperature','FontSize',10)
grid on
subplot(5,1,5)
plot(x,M_nc,'color','m','lineWidth',0.5)
xlabel('Length of the nozzle','FontSize',10)
ylabel('Mach Number','FontSize',10)
grid on
figure(4)
subplot(5,1,1)
plot(1:1400,p_throat_nc,'color','k','linewidth',1)
xlabel('Number of Time Steps','FontSize',10)
ylabel('Pressure','FontSize',10)
title('Variation of primitive at throat for time steps 1400 in Non-conservartion form','FontSize',12)
grid on
subplot(5,1,2)
plot(1:1400,rho_throat_nc,'color','k','linewidth',1)
xlabel('Number of Time Steps','FontSize',10)
ylabel('Density','FontSize',10)
grid on
subplot(5,1,3)
plot(1:1400,v_throat_nc,'color','r','linewidth',1)
xlabel('Number of Time Steps','FontSize',10)
ylabel('velocity','FontSize',10)
grid on
subplot(5,1,4)
plot(1:1400,M_throat_nc,'color','k','linewidth',1)
xlabel('Number of Time Steps','FontSize',10)
ylabel('Mach number','FontSize',10)
grid on
subplot(5,1,5)
plot(1:1400,T_throat_nc,'color','k','linewidth',1)
xlabel('Number of Time Steps','FontSize',10)
ylabel('Temperature','FontSize',10)
grid on
% Plotting conseravtive form
figure(5)
subplot (5,1,1)
plot(x,v_c,'color','r','linewidth',1.5)
xlabel('Length of the nozzle')
ylabel('velcoity')
title ('Time wise variation of primitive variables for Conservative form','FontSize',12)
grid on
subplot(5,1,2)
plot(x,rho_c,'color','r','linewidth',1.5)
xlabel('Length of the nozzle')
ylabel('Density')
grid on
subplot(5,1,3)
plot(x,p_c,'color','b','linewidth',1.5)
xlabel('Length of the nozzle')
ylabel('Pressure')
grid on
subplot(5,1,4)
plot(x,T_c,'color','g','linewidth',1.5)
xlabel('Length of the nozzle')
ylabel('Temperature')
grid on
subplot(5,1,5)
plot(x,M_c,'color','g','linewidth',1.5)
xlabel('Length of the nozzle')
ylabel('Mach Number')
grid on
figure(6)
subplot(5,1,1)
plot(1:1400,p_throat_c)
xlabel('number of time steps')
ylabel('pressure')
grid on
subplot(5,1,2)
plot(1:1400,rho_throat_c)
xlabel('number of time steps')
ylabel('density')
grid on
subplot(5,1,3)
plot(1:1400,v_throat_c)
xlabel('number of time steps')
ylabel('velocity')
grid on
subplot(5,1,4)
plot(1:1400,M_throat_c)
xlabel('number of time steps')
ylabel('mach number')
grid on
subplot(5,1,5)
plot(1:1400,T_throat_c)
xlabel('numbe of time steps')
ylabel('temperature')
grid on
figure (7)
plot(x,mfr_c,'color','k')
hold on
plot(x,mfr_nc,'color','r')
hold on
plot(x,analytical_solution,'color','g')
xlabel('length of the nozzle')
ylabel('mass flow rate')
grid on
Non Conservative function code :
function [x, M, rho, v, p, T, mfr, rho_throat,v_throat, p_throat, T_throat, M_throat, mfr_throat] = non_conser()
% Inputs
% length of the domain
L = 3;
% number of grid points
n = 31;
x = linspace(0,L,n);
dx = x(2)-x(1);
% Ratio of specfic heats (Cp/Cv) for air standard conditions
gamma = 1.4;
% calculate the initial profile
% density
rho = 1-0.3146*x;
% temperature
T = 1-0.2314*x;
% velocity
v = (0.1+1.09*x).*T.^0.5;
% nozzle shape
a = 1+2.2*(x-1.5).^2;
% Pressure
p = rho.*T;
% mass flow rate
mfr = rho.*v.*a;
% mach number
M = v./(T.^0.5);
% courrant number
C = 0.5;
% Throat values
throat = find (a==1);
% Time step
% number of time steps
nt = 1400;
% To calculate the time step
% dt = min(C*dx)./(a+v);
% alternate form
dt = min((C*dx)./((sqrt(T)+v)));
% Outer time loop
for k = 1:nt
rho_old = rho;
v_old = v;
T_old = T;
M_old = M;
p_old = p;
% Predictor method
for j = 2:n-1
dvdx = (v(j+1)-v(j))/(dx);
drhodx = (rho(j+1)-rho(j))/(dx);
dTdx = (T(j+1)-T(j))/(dx);
dlogadx = (log(a(j+1))-log(a(j)))/(dx);
% continuity equation
drhodt_p(j) = -rho(j)*dvdx-rho(j)*v(j)*dlogadx-v(j)*drhodx;
% momentum equation
dvdt_p(j) = -v(j)*dvdx-(1/gamma)*(dTdx+(T(j)/rho(j))*drhodx);
% energy equation
dTdt_p(j) = -v(j)*dTdx-((gamma-1)*T(j))*(dvdx+v(j)*dlogadx);
% solution update
rho(j) = rho(j)+drhodt_p(j)*dt;
v(j) = v(j)+dvdt_p(j)*dt;
T(j) = T(j)+dTdt_p(j)*dt;
end
% corrector step
for j = 2:n-1
dvdx = (v(j)-v(j-1))/(dx);
drhodx = (rho(j)-rho(j-1))/(dx);
dTdx = (T(j)-T(j-1))/(dx);
dlogadx = (log(a(j))-log(a(j-1)))/(dx);
% continuity equation
drhodt_c(j) = -rho(j)*dvdx-rho(j)*v(j)*dlogadx-v(j)*drhodx;
% momentum equation
dvdt_c(j) = -v(j)*dvdx-(1/gamma)*(dTdx+(T(j)/rho(j))*drhodx);
% energy equation
dTdt_c(j) = -v(j)*dTdx-((gamma-1)*T(j))*(dvdx+v(j)*dlogadx);
end
% calculate the average time derivative
drhodt = 0.5*(drhodt_p+drhodt_c);
dvdt = 0.5*(dvdt_p+dvdt_c);
dTdt = 0.5*(dTdt_p+dTdt_c);
% calculate the final solution
for i = 2:n-1
rho(i) = rho_old(i)+drhodt(i)*dt;
v(i) = v_old(i)+dvdt(i)*dt;
T(i) = T_old(i)+dTdt(i)*dt;
end
% apply the boundary conditions
% for inlet
v(1) = 2*v(2)-v(3);
% for oulet
rho(n) = 2*rho(n-1)-rho(n-2);
v(n) = 2*v(n-1)-v(n-2);
T(n) = 2*T(n-1)-T(n-2);
% calculate the pressure,mach number and mass flow rate
% Pressure
p = rho.*T;
% mass flow rate
mfr = rho.*v.*a;
% mach number
M = v./(T.^0.5);
% calculate the throat value
rho_throat(k)=rho(throat);
v_throat(k) = v(throat);
p_throat(k) = p(throat);
T_throat(k) =T(throat);
M_throat(k) = M(throat);
mfr_throat(k)= mfr(throat);
% For each time step definig the mass flow rate
figure(1)
if k == 50
plot(x,mfr)
hold on
elseif k == 100
plot(x,mfr)
hold on
elseif k == 150
plot(x,mfr)
hold on
elseif k == 200
plot(x,mfr)
hold on
elseif k == 700
plot(x,mfr)
hold on
elseif k == 1400
plot(x,mfr)
hold on
xlabel('Length of Nozzle','FontSize',12,'linewidth',1.5)
ylabel('Mass Flow Rate','FontSize',12,'linewidth',1.5)
title ('Mass Flow Rate Evolution for Non Conservative Form','FontSize',14)
grid on
end
end
end
Conservative function code :
function [x, M, rho, v, p, T, mfr, rho_throat,v_throat, p_throat, T_throat, M_throat, mfr_throat] = conservative()
% Inputs
% length of the domain
L = 3;
% number of grid points
n = 31;
x = linspace(0,L,n);
dx = x(2)-x(1);
% Ratio of specfic heats (Cp/Cv) for air standard conditions
gamma = 1.4;
% Initial conditions
% rho = 1;
% T = 1;
%
% rho = 1-0.366*(x-0.5);
% T = 1-0.167*(x-0.5);
%
% rho = 0.634-0.3879*(x-1.5);
% T = 0.833 - 0.3507*(x-1.5);
for i = 1:length(x)
if x(i)>=0 && x(i)<=0.5
rho(i) = 1;
T(i) = 1;
elseif x(i)>=0.5 && x(i)<=1.5
rho(i) = 1-0.366*(x(i)-0.5);
T(i) = 1-0.167*(x(i)-0.5);
elseif x(i)>=1.5 && x(i)<=3.5
rho(i) = 0.634-0.3879*(x(i)-1.5);
T(i) = 0.833-0.3507*(x(i)-1.5);
end
end
% Initial profile
% nozzle shape
a = 1+2.2*(x-1.5).^2;
% velocity
v = 0.59./(rho.*a);
% Pressure
p = rho.*T;
% mass flow rate
mfr = rho.*v.*a;
% mach number
M = v./(T.^0.5);
% courrant number
C = 0.5;
% To find the throat
throat =find(a==1);
% calculate the initial profile
U1 = rho.*a;
U2 = rho.*a.*v;
U3 = rho.*(T/(gamma-1)+(gamma/(2)).*v.^2).*a;
% time steps
nt = 1400;
dt = min((C*dx)./((sqrt(T)+v)));
% outer time loop
for k = 1:nt
U1_old = U1;
U2_old = U2;
U3_old = U3;
rho_old = rho;
v_old = v;
T_old = T;
% calculate Flux vectors
F1 = U2;
F2 =((U2.^2)./(U1))+((gamma-1)/(gamma))*(U3-((gamma/2)*((U2.^2)./(U1))));
F3 = ((gamma*U2.*U3)./(U1)) -((gamma*(gamma-1))/(2))*((U2.^3)./(U1.^2));
% predictor step
for j = 2:n-1
% source term
J2_p(j) = (1/gamma)*(rho(j)).*(T(j)).*(a(j+1)-a(j))/(dx);
% continuity equation
dU1dt_p(j) = -(F1(j+1)-F1(j))/(dx);
% momentum equation
dU2dt_p(j) = -(F2(j+1)-F2(j))/(dx)+J2_p(j);
% energy equation
dU3dt_p(j) = -(F3(j+1)-F3(j))/(dx);
% Solution update
U1(j) = U1(j)+dU1dt_p(j)*dt;
U2(j) = U2(j)+dU2dt_p(j)*dt;
U3(j) = U3(j)+dU3dt_p(j)*dt;
end
% Updating the primitive variables
rho = U1./a;
v = U2./U1;
T = (gamma-1)*(U3./(U1)-(gamma/(2))*v.^2);
% Updating the flux vectors values
F1 = U2;
F2 =((U2.^2)./(U1))+((gamma-1)/(gamma))*(U3-((gamma/2)*((U2.^2)./(U1))));
F3 = ((gamma*U2.*U3)./(U1))-((gamma*(gamma-1))/(2))*((U2.^3)./(U1.^2));
% Corrector step
for j = 2:n-1
% source term
J2_c(j) = (1/gamma)*(rho(j)).*(T(j)).*(a(j)-a(j-1))/(dx);
% continuity equation
dU1dt_c(j) = -(F1(j)-F1(j-1))/(dx);
% momentum equation
dU2dt_c(j) = -(F2(j)-F2(j-1))/(dx)+J2_c(j);
% energy equation
dU3dt_c(j) = -(F3(j)-F3(j-1))/(dx);
end
% calculate the average time derivative
dU1dt = 0.5*(dU1dt_p+dU1dt_c);
dU2dt = 0.5*(dU2dt_p+dU2dt_c);
dU3dt = 0.5*(dU3dt_p+dU3dt_c);
% final solution
for i = 2:n-1
U1(i) = U1_old(i)+dU1dt(i)*dt;
U2(i) = U2_old(i)+dU2dt(i)*dt;
U3(i) = U3_old(i)+dU3dt(i)*dt;
end
% Applying the boundary conditions
% Inlet
U1(1) = rho(1).*a(1);
U2(1) = 2*U2(2)-U2(3);
U3(1) = U1(1)*(T(1)./(gamma-1))+(gamma/(2))*(v(1).^2);
% Outlet
U1(n) = 2*U1(n-1)-U1(n-2);
U2(n) = 2*U2(n-1)-U2(n-2);
U3(n) = 2*U3(n-1)-U3(n-2);
% calculate the primitive variables
rho = U1./a;
v = U2./U1;
T = (gamma-1)*(U3./(U1)-(gamma/(2))*v.^2);
% Pressure
p = rho.*T;
% mass flow rate
mfr = rho.*v.*a;
% mach number
M = v./(T.^0.5);
% calculate the throat value
rho_throat(k)=rho(throat);
v_throat(k) = v(throat);
p_throat(k) = p(throat);
T_throat(k) =T(throat);
M_throat(k) = M(throat);
mfr_throat(k)= mfr(throat);
% For each time step definig the mass flow rate
figure(2)
if k == 50
plot(x,mfr)
hold on
elseif k == 100
plot(x,mfr)
hold on
elseif k == 150
plot(x,mfr)
hold on
elseif k == 200
plot(x,mfr)
hold on
elseif k == 700
plot(x,mfr)
hold on
elseif k == 1400
plot(x,mfr)
hold on
xlabel('Length of Nozzle','FontSize',12,'linewidth',1.5)
ylabel('Mass Flow Rate','FontSize',12,'linewidth',1.5)
title ('Mass Flow Rate Evolution for Conservative Form ','FontSize',14)
grid on
end
end
end
Results :
Figure 1 (Mass Flow Rate Evolution for Non-Conservative Form):
- Mass fow rates are plotted against the length of the nozzle for various time steps during the time marching process using the non-conservative form.
- The x-axis replresents the length of the nozzle, and the y-axis represents the mass flow rate.
- The plot is updated at different time steps (e.g., at time steps 50, 100, 150, 200, 700, and 1400).
Observations:
- Check how the mass flow rate evolves over time, particularly focusing on the nozzle's throat region.
- Observe whether the mass flow rate reaches a steady state or if there are any transient behaviors.
Figure 2 (Mass Flow Rate Evolution for Conservative Form):
- Similar to Figure 1, this figure represents the evolution of mass flow rates for various time steps but using the conservative form.
- The x-axis represents the length of the nozzle, and the y-axis represents the mass flow rate.
- The plot is updated at different time steps during the time marching process.
Observations:
- Compare the mass flow rate evolution between the non-conservative and conservative forms.
- Analyze whether the conservative form yields different results or converges differently compared to the non-conservative form.
Non-Conservative Form Results (Figure 3 and Figure 4):
- Primitive Variables: Velocity, Density, Pressure, Temperature, and Mach Number are plotted against the length of the nozzle and for various time steps.
- The plots show the evolution of these variables over time using the non-conservative form of the governing equations.
- Throat Conditions: The primitive variables at the throat are also visualized for different time steps.
Conservative Form Results (Figure 5 and Figure 6):
- Similar to the non-conservative form, primitive variables are plotted against the length of the nozzle and for various time steps.
- Throat Conditions: The primitive variables at the throat are visualized over time using the conservative form.
Mass Flow Rate Variation (Figure 7):
- Mass flow rates are plotted against the length of the nozzle for both non-conservative and conservative forms.
- The analytical solution is also included for reference.
- Observations: Check whether the mass flow rates converge to a steady-state value and how well the numerical results match the analytical solution.
