Simulation of a 1D Super-sonic nozzle flow simulation using Macormack Method
Simulation of 1D Super-sonic nozzle flow simulation using Macormack method:
The equations which are used to simulate the flow can be found in computational fluid dynamics by John D Anderson.
Main Code :
% MATLAB program to simulate 1D supersonic nozzle flow using macormack method
clear all
close all
clc
% inputs
n = 31;
x = linspace(0,3,n);
dx = x(2)-x(1);
gamma = 1.4;
% area
a = 1+2.2*(x-1.5).^2;
% courant number
C = 0.5;
nt = 1400;% number of time steps
% conservative solver
tic
[rho_con,rho_con_throat,T_con,T_con_throat,v_con,p_con,p_throat,mass_flow_con,M_con,M_con_throat] = conserve(x,gamma,nt,n,C,dx,a)
time_c = toc;
fprintf(\'computationla time %f\',time_c);
% time wise variation of density, temperature, mach number,pressure for conservative form
time_var_con = figure(3);
subplot(4,1,1)
plot(rho_con_throat,\'b\');
title(\' time-wise variation of properties at the nozzle throat(conservative form)\');
ylabel(\'rho/rho_{0}\');
axis([0 1400 0.5 1]);
legend(\'density\');
grid minor;
subplot(4,1,2)
plot(T_con_throat,\'r\');
ylabel(\'T/T_{0}\');
axis([0 1400 0.6 1.2]);
legend(\'temperature\');
grid minor;
subplot(4,1,3)
plot(p_throat);
ylabel(\'p/p_{0}\');
axis([0 1400 0.4 1]);
legend(\'pressure\');
grid minor;
subplot(4,1,4)
plot(M_con_throat);
ylabel(\'M\');
legend(\'mach number\');
axis([0 1400 0.8 1.4])
grid minor;
xlabel(\'number of iterations\');
% properties profiles in 1D Supersonic nozzle flow
qualitative = figure(4);
hold on
subplot(4,1,1)
plot(x,rho_con);
title(\'properties profile in the 1D Supersonic Nozzle flow(conservative form)\');
ylabel(\'rho/rho_{0}\');
legend(\'density\');
grid minor;
subplot(4,1,2)
plot(x,p_con);
ylabel(\'p/p_{0}\');
legend(\'pressure\');
grid minor;
subplot(4,1,3)
plot(x,T_con);
ylabel(\'T/T_{0}\');
legend(\'temperature\');
grid minor;
subplot(4,1,4)
plot(x,M_con);
ylabel(\'M\');
legend(\'mach number\');
grid minor
% non conservative solver
tic
[rho,T,v,M_non,p_non,massflow_non,M_non_throat,p_non_throat,rho_throat,T_non_throat] = non_conservative(x,n,gamma,nt,C,dx,a);
time_non = toc;
fprintf(\'computationla time %f/n\',time_non);
% time wise variation of density, temperature, mach number,pressure for non conservative form
time_var_con = figure(5);
subplot(4,1,1)
plot(rho_throat,\'b\');
title(\'time-wise properties at the nozzle throat (non conservative form)\');
ylabel(\'rho/rho_{0}\');
axis([0 1400 0.5 1]);
legend(\'density\');
grid minor;
subplot(4,1,2)
plot(T_non_throat,\'r\');
ylabel(\'T/T_{0}\');
axis([0 1400 0.6 1.2]);
legend(\'temperature\');
grid minor;
subplot(4,1,3)
plot(p_non_throat);
ylabel(\'p/p_{0}\');
axis([0 1400 0.4 1]);
legend(\'pressure\');
grid minor;
subplot(4,1,4)
plot(M_non_throat);
ylabel(\'M\');
legend(\'mach number\');
axis([0 1400 0.8 1.4])
grid minor;
xlabel(\'number of iterations\');
% properties profiles
figure(6)
subplot(4,1,1)
plot(x,M_non,\'-r\');
title(\'properties profiles of 1D Supersonic nozzle flow( non conservative form)\');
ylabel(\'M\');
legend(\'mach number\');
grid minor;
subplot(4,1,2)
plot(x,p_non);
ylabel(\'p/p_{0}\');
legend(\'pressure\');
grid minor;
subplot(4,1,3)
plot(x,T);
ylabel(\'T/T_{0}\');
legend(\'temperature\');
grid minor;
subplot(4,1,4)
plot(x,rho);
ylabel(\'rho/rho_{0}\');
legend(\'density\');
grid minor;
% comparision of mass flow rates by conservative and non conservative forms
mass_comparision = figure(7);
hold on
plot(x,massflow_non);
plot(x,mass_flow_con);
legend(\'Conservative\',\'Non conservative\');;
grid minor;
xlabel(\'x/L\');
ylabel(\'rhova/rho_{0}v_{0}a_{0}\');
title(\'mass flow comparision between conservative and non conservative forms\');Conservative code:
function [rho_con,rho_con_throat,T_con,T_con_throat,v_con,p_con,p_throat,mass_flow_con,M_con,M_con_throat] = conserve(x,gamma,nt,n,C,dx,a)
% initial profiles
for i = 1:length(x)
if x(i)>=0 && x(i)<=0.5
rho_con(i) = 1;
T_con(i) = 1;
elseif x(i)>=0.5 && x(i)<=1.5
rho_con(i) = 1-0.366*(x(i)-0.5);
T_con(i) = 1-0.167*(x(i)-0.5);
elseif x(i)>=1.5 && x(i)<=3.5
rho_con(i) = 0.634-0.3879*(x(i)-1.5);
T_con(i) = 0.833-0.3507*(x(i)-1.5);
endif
endfor
v_con = 0.59./(rho_con.*a);
p_con = rho_con.*T_con;
% calculating initial solution vector
U1 = rho_con.*a;
U2 = rho_con.*a.*v_con;
U3 = rho_con.*a.*((T_con./(gamma-1))+(gamma/2)*(v_con.^2));
for k = 1:nt
% copy of old solution vectors
U1_old = U1;
U2_old = U2;
U3_old = U3;
% calculating flux vectors
F1 = U2;
F2 = (U2.^2./U1)+((gamma-1)/gamma)*(U3-(gamma.*U2.^2./(2.*U1)));
F3 = (gamma.*U2.*U3./U1)-(gamma*(gamma-1)/2)*(U2.^3./U1.^2);
% calculating time step
a_con = sqrt(T_con);
timestep = C*dx./(a_con+v_con);
dt = min(timestep);
% predictor step
for j = 2:n-1
J2(j) = (1/gamma).*rho_con(j).*T_con(j).*((a(j+1)-a(j))/dx);
dU1dt_p(j) = -(F1(j+1)-F1(j))/dx;
dU2dt_p(j) = -((F2(j+1)-F2(j))/dx)+J2(j);
dU3dt_p(j) = -(F3(j+1)-F3(j))/dx;
% updating solution vector
U1(j) = U1(j)+dU1dt_p(j)*dt;
U2(j) = U2(j)+dU2dt_p(j)*dt;
U3(j) = U3(j)+dU3dt_p(j)*dt;
endfor
% updating variables and flux vectors
rho_con = U1./a;
v_con = U2./U1;
T_con = (gamma-1)*(U3./U1-(gamma/2)*(U2./U1).^2);
p_con = rho_con.*T_con;
F1 = U2;
F2 = (U2.^2./U1)+((gamma-1)/gamma)*(U3-(gamma*U2.^2./(2.*U1)));
F3 = (gamma*U2.*U3./U1)-(gamma*(gamma-1)/2)*(U2.^3./U1.^2);
% corrector method
for j = 2:n-1
J2_c(j) = (1/gamma)*rho_con(j).*T_con(j).*((a(j)-a(j-1))/dx);
dU1_dt_c(j) = -(F1(j)-F1(j-1))/dx;
dU2_dt_c(j) = -((F2(j)-F2(j-1))/dx)+J2_c(j);
dU3_dt_c(j) = -(F3(j)-F3(j-1))/dx;
endfor
% calculating average time derivative
U1_avg = 0.5*(dU1dt_p + dU1_dt_c);
U2_avg = 0.5*(dU2dt_p + dU2_dt_c);
U3_avg = 0.5*(dU3dt_p + dU3_dt_c);
% final solution update
for i = 2:n-1
U1(i) = U1_old(i)+U1_avg(i)*dt;
U2(i) = U2_old(i)+U2_avg(i)*dt;
U3(i) = U3_old(i)+U3_avg(i)*dt;
endfor
% boundary conditions
% inlet
rho_con(1) = 1;
T_con(1) = 1;
U1(1) = rho_con(1).*a(1);
U2(1) = 2*U2(2)-U2(3);
v_con(1) = U2(1)./U1(1);
U3(1) = U1(1)*(T_con(1)/(gamma-1))+(gamma/2)*v_con(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);
% updating variables
rho_con = U1./a;
v_con = U2./U1;
T_con = (gamma-1)*(U3./U1-(gamma/2)*(U2./U1).^2);
p_con = rho_con.*T_con;
% mass flow rate
mass_flow_con = rho_con.*v_con.*a;
M_con = v_con./sqrt(T_con); % mach number
p_con = rho_con.*T_con;
% properties at throat
rho_con_throat(k) = rho_con(((n-1)/2)+1);
p_throat(k) = p_con(((n-1)/2)+1);
T_con_throat(k) = T_con(((n-1)/2)+1);
M_con_throat(k) = M_con(((n-1)/2)+1);
% comparision of non dimensional mass flow rate for different time steps
mass = figure(1);
hold on
if k == 100
plot(x,mass_flow_con,\'b\');
elseif k == 200
plot(x,mass_flow_con,\'g\');
elseif k == 700
plot(x,mass_flow_con,\'v\');
end
endfor
title(\'transient non dimensional mass flow (conservative form)\');
xlabel(\'x/L\');
ylabel(\'conservative mass flow rate\');
legend(\'100 time step\',\'200 time step\',\'700 time step\');
endfunction
Non Conservative Code:
function [rho,T,v,M_non,p_non,massflow_non,M_non_throat,p_non_throat,rho_throat,T_non_throat] = non_conservative(x,n,gamma,nt,C,dx,a)
% calculate input profles
rho = 1-0.3146*x;
T = 1-0.2314*x;
v = (0.1+1.09*x).*T.^(1/2);
% time steps
a_non = sqrt(T);
timestep_non = C*dx./(a_non+v);
dt = min(timestep_non);
for k2 = 1:nt
rho_old = rho;
v_old = v;
T_old = T;
% predictor method
for j = 2:n-1
% drho_dt_p = -rho(j)*(v(j+1)-v(j))/dx-rho(j)*v(j)*(log(A(j+1))-log(A(j)))/dx-v(j)*(rho(j+1)-rho(j))/dx;
dvdx = (v(j+1)-v(j))/dx;
drhodx = (rho(j+1)-rho(j))/dx;
dlogadx = (log(a(j+1))-log(a(j)))/dx;
dTdx = (T(j+1)-T(j))/dx;
% continuty equation
drho_dt_p(j) = -rho(j)*dvdx -rho(j)*v(j)*dlogadx -v(j)*drhodx;
% momentum equation
dv_dt_p(j) = - v(j)*dvdx -(1/gamma)*(dTdx+(T(j)/rho(j))*drhodx);
% energy equation
dT_dt_p(j) = -v(j)*dTdx -(gamma-1)*T(j)*(dvdx+v(j)*dlogadx);
% solution update
rho(j) = rho(j)+drho_dt_p(j)*dt;
v(j) = v(j)+dv_dt_p(j)*dt;
T(j) = T(j)+dT_dt_p(j)*dt;
end
% correcto step
for j = 2:n-1
dvdx = (v(j)-v(j-1))/dx;
drhodx = (rho(j)-rho(j-1))/dx;
dlogadx = (log(a(j))-log(a(j-1)))/dx;
dTdx = (T(j)-T(j-1))/dx;
% continuty equation
drho_dt_c(j) = -rho(j)*dvdx -rho(j)*v(j)*dlogadx -v(j)*drhodx;
% momentum equation
dv_dt_c(j) = - v(j)*dvdx -(1/gamma)*(dTdx+(T(j)/rho(j))*drhodx);
% energy equation
dT_dt_c(j) = -v(j)*dTdx -(gamma-1)*T(j)*(dvdx+v(j)*dlogadx);
end
% calculating average time derivative
drho_dt_avg = 0.5*(drho_dt_p+drho_dt_c);
dv_dt_avg = 0.5*(dv_dt_p+dv_dt_c);
dT_dt_avg = 0.5*(dT_dt_p+dT_dt_c);
% next time step
for i = 2:n-1
rho(i) = rho_old(i)+drho_dt_avg(i)*dt;
v(i) = v_old(i)+dv_dt_avg(i)*dt;
T(i) = T_old(i)+dT_dt_avg(i)*dt;
end
% apply boundary conditions
% inlet
v(1) = 2*v(2)-v(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);
% mass flow rate
massflow_non = rho.*v.*a;
p_non = rho.*T;
M_non = v./sqrt(T);
% throat values
p_non_throat(k2) = p_non(((n-1)/2)+1);
M_non_throat(k2) = M_non(((n-1)/2)+1);
T_non_throat(k2) = T(((n-1)/2)+1);
rho_throat(k2) = rho(((n-1)/2)+1);
% comparision of non dimensional mass flow rate for different time steps
mass2 = figure(2);
hold on
if k2 == 50
plot(x,massflow_non,\'r\');
elseif k2 == 100
plot(x,massflow_non,\'b\');
elseif k2 == 150
plot(x,massflow_non,\'g\');
elseif k2 == 200
plot(x,massflow_non,\'c\');
elseif k2 == 700
plot(x,massflow_non,\'y\');
endif
endfor
title(\'transient non dimensional mass flow (non conservative form)\');
xlabel(\'x/L\');
ylabel(\'non conservative mass flow rate\');
legend(\'50 timestep\',\'100 time step\',\'150 timestep\',\'200 time step\',\'700 time step\');
endfunctionRESULTS:
(i) Time wise properties at nozzle throat
(a) Conservative form
(b) non conservative form
Non dimensionalized parameters like density,pressure,temperature and mach number are plotted at throat point against the time step to check at which time step steady state reaches for both conservative and non conservative forms. Initially for non conservative it shows greater variations but soon it stable and simillarly for conservative also it shows small variations and even it also reaches steady state soon. from the above results we can observe that for both non conservative and conservative steady state reaches after 700 iterations approx.
(ii) Properties profiles of 1D Supersonic Nozzle flow with length
(a) Conservative form
(b) Non conservative form
we can observe from the plots that for both conservative and non conservative the density, pressure and temperature are dropping at the exit of the nozzle and the mach number is increasing as the acceleration is increasing. These results are almost identical to the profiles in the textbook (John D Anderson).
(iii) Non dimensional mass flow rate with various time steps
(a) Conservative form:
(b) Non conservative form:
It can be seen from both figures that the mass flow rate becomes steady after 700 time steps for both conservative and non conservative .
(iv) comparision of non dimensional mass flow rates of both forms
from the plot we can observe that non conservative form has more fluctuations when compared to conservative form thus non conservative is less satisfactory when compared to conservative
Grid Dependence test:
At throat point the values non dimensionalized parameters for both conservative and non conservative are given below
for n = 31 Pressure Mach number Temperature density
conservative 0.545 0.981 0.838 0.650
non-conservative 0.534 0.999 0.836 0.638
exact 0.528 1.000 0.833 0.634
for n = 61 Pressure Mach number Temperature density
conservative 0.531 0.9993 0.833 0.638
non-conservative 0.532 0.9997 0.835 0.637
exact 0.528 1.000 0.833 0.634
as we can see there is not much difference between the parameters when the grid value changes from n = 31 to 61 and also the percentage of error is also small between the numerical and exact solution.
Conservative runs faster when compared to non conservative
computational time for conservative is 15.788 seconds whereas the computational time for non conservative is 18.810 seconds