Week 7 - Simulation of a 1D Super-sonic nozzle flow simulation using Macormack Method
Quasi 1D analysis of a supersonic nozzle using in both the conservation and non-conservation forms of the governing equations and sovle them using the MacCormack's technique:
Convergent divergent nozzle or the supersonic nozzle or the delaval nozzle is used for supersonic flows. This has three sections: Subsonic section, throat, supersonic section.
Fluid of low subsonic speed enter the subsonic section and the speed increases to the high subsonic speed at the end of the section. Fluid further increases the speed and reaches sonic speed at the throat. In the supersonic diverging section speed of the fluid becomes supersonic(Ma>1).
This nozzle is modelled using MacCormack's Method in both conservative and non-conservative forms of govering equations.
Main program:
close all
clear all
clc
%inputs
n=31;
x=linspace(0,3,n);
dx=x(2)-x(1);
gamma=1.4;
A = 1+ (2.2*(x-1.5).^2);
n_throat = find((A-1)<0.001);
CFL=0.5;
%Non conservative form
[rho,v,t,p,M,m,p_throat,rho_throat,v_throat,t_throat,M_throat,k,t_simulation]=project2_nonconservative(n,x,dx,A,gamma,n_throat);
time_nonconservative = t_simulation;
fprintf('Computational time:Non-Conservative = %f sec',t_simulation);
fprintf('\n Number of iterations = %f',k);
%Conservative form
[rho_ini,v_ini,T_ini,p_ini,M_ini,mdot,rho_throat_ini,v_throat_ini,T_throat_ini,M_throat_ini,p_throat_ini,t_sim]=project2_conservative(n,x,dx,A,gamma,n_throat);
time_conservative = t_sim;
fprintf('Computational time:Non-Conservative = %f sec',t_sim);
fprintf('\n Number of iterations = %f',k);Non Conservative form:
%non-conservative form of 1D supersonic nozzle flow equations
function [rho,v,t,p,M,m,p_throat,rho_throat,v_throat,t_throat,M_throat,k,t_simulation]=project2_nonconservative(n,x,dx,A,gamma,n_throat)
%calculate initial profiles
rho=1-0.3146*x;
t=1-0.2134*x; %t=temp
v=(0.1+1.09*x).*t.^0.5;
a=1+(2.2*(x-1.5).^2);
%time steps
nt=700;
dt=0.01;
t_simulation=(nt-1)*dt;
time=0:dt:t_simulation;
%outer time loop
for k=1:nt
rho_old=rho;
v_old=v;
t_old=t;
%predictor time step
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 method
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
%compute 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);
for i=2:n-1
%final solution update
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 boundary condition
%inlet
rho(1)= 2*rho(2)- rho(3);
v(1)= 2*v(2)- v(3);
t(1)= 2*t(2)- t(3);
%outlet
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);
%Defining the non dimensional values of pressure and mass flow rate,Mach
%number
p=rho.*t;
m = rho.*a.*v; %mass flow rate
M = v./sqrt(t);
%mdot=rho.*v.*A;
%Parameters at the throat
t_throat(k) = t(n_throat);
p_throat(k) = p(n_throat);
rho_throat(k) = rho(n_throat);
v_throat(k) = v(n_throat);
M_throat(k)= M(n_throat);
%Throat section residual values
rho_r(k) = drhodt(n_throat);
t_r(k)= dTdt(n_throat);
v_r(k)=dvdt(n_throat);
%Comparision of mass flow rates at different times
mass=figure(1);
hold on;
if k ==50
plot(x,m,'r')
elseif k ==100
plot(x,m);
elseif k == 150
plot(x,m,'b')
elseif k==200
plot(x,m,'g')
elseif k==700
plot(x,m)
end
end
legend('t=0','t=50','t=100','t=150','t=200');
xlabel('x/L')
ylabel('mass flow rate')
title('Change in mass flow rate with time')
%throat section
figure(2);
%density ratio
subplot(4,1,1);
plot(rho_throat,'b');
title('Transient Analysis of Properties at Throat');
legend('Density Ratio');
ylabel('rho/rho0');
axis([0 700 0.5 1]);
%Temperature
subplot(4,1,2);
plot(t_throat,'m');
legend('Temperature');
ylabel('T/T0');
axis([0 700 0.5 1]);
%Pressure
subplot(4,1,3);
plot(p_throat,'g');
legend('Pressure');
ylabel('P/P0');
axis([0 700 0.3 1.5]);
%Mach number
subplot(4,1,4);
plot(M_throat,'r');
legend('Mach number');
ylabel('M/M0');
axis([0 700 0.8 1.5]);
xlabel('Number of Iterations');
%Steady State Profiles
figure(3);
%Density Ratio
subplot(4,1,1);
plot(x,rho,'g');
title('Steaty state profiles');
legend('Density Ratio');
ylabel('rho/rho0');
%Temperature
subplot(4,1,2);
plot(x,t,'m');
legend('Temperature');
ylabel('T/T0');
%Pressure
subplot(4,1,3);
plot(x,p,'r');
legend('Pressure');
ylabel('p/p0');
%Mach number
subplot(4,1,4);
plot(x,M,'b');
legend('Mach number');
ylabel('M/M0');
xlabel('length');
Change in mass flow rate is studied with respect to length at different times. Non-dimensional variables are taken for this purpose.
Changes in various properties(temperature, pressure, mach number and density) with respect to time are studied.
From the above properties, we observe that these variables increase first and then decreases with time.
Steady state profiles of the flow are as shown below:
In this profile we can observe clearly the change in properties with respect to length. From the Mach number profile, it can be stated that velocity increases as we move forward in x direction. Till the throat section, the increase in velocity is not much observed. But as it enters the supersonic section, velocity increases much faster than the subsonic and throat section. In the similar way, pressure,density and temperature decreases.
Total simulation time is 6.99 seconds
Conservative form:
%Conservative form of 1D supersonic nozzle flow equations
function [rho_ini,v_ini,T_ini,p_ini,M_ini,mdot,rho_throat_ini,v_throat_ini,T_throat_ini,M_throat_ini,p_throat_ini,t_sim]=project2_conservative(n,x,dx,A,gamma,n_throat)
%calculate initial profiles
dt = 0.001;
nt=1400;
CFL=0.5;
t_sim=(nt-1)*dt;
%time=0:dt:t_sim;
for i=1:length(x)
if x(i)>=0 && x(i)<0.5
rho_ini(i) = 1;
T_ini(i) = 1;
elseif x(i)>= 0.5 && x(i)< 1.5
rho_ini(i) = 1- 0.366*(x(i)-0.5);
T_ini(i)= 1- 0.167*(x(i)-0.5);
elseif x(i)>= 1.5 && x(i) <=3
rho_ini(i) = 0.634 - 0.3879*(x(i)-1.5);
T_ini(i) = 0.833 - 0.3507*(x(i)-1.5);
end
end
rho=rho_ini;
T=T_ini;
v_ini = (0.59)./(rho_ini.*A);
U=v_ini;
rho_old_ini = rho_ini;
v_old_ini = v_ini;
T_old_ini = T_ini;
% Initial conditions
U1 = rho_ini.*A;
U2 = rho_ini.*A.*v_ini;
U3 = rho_ini.*A.*((T_ini./(gamma-1))+((gamma/2)*(v_ini.^2)));
%Time loop
for j = 1:nt
U1_old = U1;
U2_old = U2;
U3_old = U3;
dt_ini = CFL*(dx./(sqrt(T_ini)+v_ini));
dt = min(dt_ini);
F1 = U2;
F2 = ((U2.^2)./U1) + (((gamma-1)/gamma)*(U3 - ((gamma*0.5*(U2.^2))./U1)));
F3 = ((gamma*U2.*U3)./U1) - ((gamma*0.5*(gamma-1)*(U2.^3))./(U1.^2));
%Predictor step
for k = 2:n-1
dAdxj(k) = (1/gamma)*rho_ini(k)*T_ini(k)*((A(k+1)-A(k))/dx);
%Continuity equation
dU1dt_ini_p(k) = -(F1(k+1)-F1(k))/dx;
%Momentum equation
dU2dt_ini_p(k) = -((F2(k+1)-F2(k))/dx) + dAdxj(k);
%Energy equation
dU3dt_ini_p(k) = -(F3(k+1)-F3(k))/dx;
%Predictor step update
U1(k) = U1(k) + (dU1dt_ini_p(k)*dt);
U2(k) = U2(k) + (dU2dt_ini_p(k)*dt);
U3(k) = U3(k) + (dU3dt_ini_p(k)*dt);
end
%Predicted values
rho_ini = U1./A;
v_ini = U2./U1;
T_ini = (gamma-1)*((U3./U1)-((gamma/2)*((U2./U1).^2)));
p_ini = rho_ini.*T_ini;
F1 = U2;
F2 = ((U2.^2)./U1) + (((gamma-1)/gamma)*(U3 - ((gamma*0.5*(U2.^2))./U1)));
F3 = ((gamma*U2.*U3)./U1) - ((gamma*0.5*(gamma-1)*(U2.^3))./(U1.^2));
%corrector step
for k = 2:n-1
dAdxj(k) = (1/gamma)*rho_ini(k)*T_ini(k)*((A(k)-A(k-1))/dx);
% Continuity Equation
dU1dt_ini_c(k) = -(F1(k)-F1(k-1))/dx;
%Momentum Equation
dU2dt_ini_c(k) = -(F2(k)-F2(k-1))/dx + dAdxj(k);
%Energy Equation
dU3dt_ini_c(k) = -(F3(k)-F3(k-1))/dx;
end
%Average Time derivative
dU1dt_ini_avg = 0.5*(dU1dt_ini_p + dU1dt_ini_c);
dU2dt_ini_avg = 0.5*(dU2dt_ini_p + dU2dt_ini_c);
dU3dt_ini_avg = 0.5*(dU3dt_ini_p + dU3dt_ini_c);
%Final Solution Update
for i = 2:n-1
U1(i) = U1_old(i) + (dU1dt_ini_avg(i)*dt);
U2(i) = U2_old(i) + (dU2dt_ini_avg(i)*dt);
U3(i) = U3_old(i) + (dU3dt_ini_avg(i)*dt);
end
%Boundary Conditions
%Inlet
rho_ini = U1./A;
v_ini = U2./U1;
T_ini = (gamma-1)*(U3./U1-(gamma/2)*(v_ini.^2));
U1(1) = rho_ini(1)*A(1);
U2(1) = 2*U2(2)-U2(3);
U3(1) = U1(1)*(T_ini(1)/(gamma-1)+(gamma/2)*(v_ini(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);
%Calculation of properties
p_ini = rho_ini.*T_ini;
M_ini = v_ini./sqrt(T_ini); %Mach number
mdot = rho_ini.*A.*v_ini;
%Properties at Throat
rho_throat_ini(j) = rho_ini(n_throat);
T_throat_ini(j) = T_ini(n_throat);
M_throat_ini(j) = M_ini(n_throat);
v_throat_ini(j) = v_ini(n_throat);
p_throat_ini(j) = p_ini(n_throat);
%Plotting change in mass-flow
figure(4);
hold on;
if j==1
plot(x,mdot,'*');
elseif j==50
plot(x,mdot,'--');
elseif j==100
plot(x,mdot,'b');
elseif j==150
plot(x,mdot,'r');
elseif j==200
plot(x,mdot,'g');
elseif j==500
plot(x,mdot,'y');
elseif j==700
plot(x,mdot);
end
end
title('Mass flow rate, mdot:Conservative Approach')
xlabel('x/L');
ylabel('m/m0');
legend('t=0','t=50','t=100','t=150','t=200','t=500','t=700')
axis([0 3 0 2]);
%grid minor;
%Change in properties at throat
figure(5);
% %Density ratio
subplot(4,1,1);
plot(rho_throat_ini,'b');
title('Change in properties at throat with time');
legend('Density ratio');
ylabel('rho/rho0');
axis([0 700 0.5 1]);
% %Temperature
subplot(4,1,2);
plot(T_throat_ini,'g');
legend('Temperature');
ylabel('T/T0');
axis([0 700 0.3 1]);
%Pressure
subplot(4,1,3);
plot(p_throat_ini,'m');
legend('Pressure');
ylabel('P/P0');
axis([0 700 0.2 1]);
%Mach number
subplot(4,1,4);
plot(M_throat_ini,'m');
legend('Mach number');
ylabel('M/M0');
axis([0 700 0.8 1.5]);
xlabel('Number of iterations');
%Steady state profiles of nozzle
figure(6);
%Density Ratio
subplot(4,1,1);
plot(x,rho_ini);
title('Steady state profiles of nozzle');
legend('Density ratio');
ylabel('rho/rho0');
%Temperature
subplot(4,1,2);
plot(x,T_ini);
legend('Temperature');
ylabel('T/T0');
%Pressure
subplot(4,1,3);
plot(x,p_ini);
legend('Pressure');
ylabel('p/p0');
%Mach number
subplot(4,1,4);
plot(x,M_ini);
legend('Mach Number');
ylabel('M/M0');
xlabel('Length');
end Simulation time=0.699sec
Variation of mass flow rate with distance:
Plotting properties at throat:
From the above plot, it can be concluded that the solution is converged.
Steady state profiles of nozzle:
Density, Temperature and pressure ratios decrease with length and mach number increases. This shows that the velocity increases and supersonic flow come out of the nozzle.
Grid independence test:
Grid independence test means, the solution must not show much variation with the change in grid points. This is calculated for 700 iterations for both conservative and non-conservative form.
Conservative:
| Grid points | Pressure ratio | Density ratio | Temperature ratio | Velocity ratio |
| n=31 | 0.5470 | 0.6509 | 0.8403 | 0.8989 |
| n=61 | 0.5322 | 0.6381 | 0.8340 | 0.9111 |
| n=91 | 0.5289 | 0.6346 | 0.8335 | 0.9180 |
Difference in the variables is 1e-1 for grid points between 61 and 91. Based on the accuracy required and computational time, number of grid points are selected.
Non-conservative:
| Grid points | Pressure ratio | Density ratio | Temperature ratio | Velocity ratio |
| n=31 | 0.0660 | 0.1436 | 0.4533 | 0.6716 |
| n=61 | 0.3279 | 0.4493 | 0.7285 | 0.8420 |
| n=91 | 0.4333 | 0.5470 | 0.7921 | 0.8847 |
| n=111 | 0.4668 | 0.5766 | 0.8094 | 0.8960 |
Solution variation is less between grid points 91 and 111. Any feasible number can be choosen.