1D supersonic nozzle flow simulation using macormack method

Matlab code for solving 1D supersonic nozzle flow using macormack method:

clear all 
close all
clc

% grid space
n=31;
x=linspace(0,3,n);
dx=x(2)-x(1);
gamma = 1.4;


%function for non conserved form
tic
[rho,T,v,A,M,P,m,rho_th,T_th,P_th,M_th,rho_res,T_res,v_res,k_n] = supersonic_non_conserved(n,x,dx,gamma);
time_for_non_conservative_solver = toc
%number of cycles for convergence
c_1 = k_n

hold on

%function for conservative form
tic
[rho_c,T_c,v_c,A_c,M_c,P_c,m_c,rho_th_c,T_th_c,P_th_c,M_th_c,rho_res_c,T_res_c,v_res_c,k_c] = supersonic_conserved(n,x,dx,gamma);
time_for_conservative_solver = toc
%number of cycles for convergence
c_2 = k_c


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%--------------------plots for non conservative form----------------------%
%steady state plot
figure(3)
title(\'steady state plot\')
subplot(4,1,1)
plot(x,M,\'linewidth\',1)
xlabel(\'x/Xo\')
ylabel(\'mach number\')
subplot(4,1,2)
plot(x,P,\'k\',\'linewidth\',1)
xlabel(\'x/Xo\')
ylabel(\'pressure ratio\')
subplot(4,1,3)
plot(x,rho,\'r\',\'linewidth\',1)
xlabel(\'x/Xo\')
ylabel(\'density\')
subplot(4,1,4)
plot(x,T,\'g\',\'linewidth\',1)
xlabel(\'x/Xo\')
ylabel(\'temperature\')



%varibles at throat
figure(4)
title(\'variables at throat section\')
subplot(4,1,1)
plot(rho_th,\'b\',\'linewidth\',1)
ylabel(\'rho\')
subplot(4,1,2)
plot(T_th,\'k\',\'linewidth\',1)
ylabel(\'temp\')
subplot(4,1,3)
plot(P_th,\'c\',\'linewidth\',1)
ylabel(\'P\')
subplot(4,1,4)
plot(M_th,\'r\',\'linewidth\',1)
xlabel(\'number of time steps\')
ylabel(\'mach\')

%residual values
figure(5)
plot(rho_res,\'b\',\'linewidth\',1)
hold on
plot(v_res,\'r\',\'linewidth\',1)
hold on
plot(T_res,\'--k\',\'linewidth\',1)
xlabel(\'number of time steps\')
ylabel(\'residual\')
legend(\'density\',\'velocity\',\'temperature\')
%-------------------------------------------------------------------------%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%---------------plots for conservative form-------------------------------%
%steady state plot
figure(6)
title(\'steady state plot\')
subplot(4,1,1)
plot(x,M_c,\'linewidth\',1)
xlabel(\'x/Xo\')
ylabel(\'mach number\')
subplot(4,1,2)
plot(x,P_c,\'k\',\'linewidth\',1)
xlabel(\'x/Xo\')
ylabel(\'pressure ratio\')
subplot(4,1,3)
plot(x,rho_c,\'r\',\'linewidth\',1)
xlabel(\'x/Xo\')
ylabel(\'density\')
subplot(4,1,4)
plot(x,T_c,\'g\',\'linewidth\',1)
xlabel(\'x/Xo\')
ylabel(\'temperature\')


%varibles at throat
figure(7)
title(\'variables at throat section\')
subplot(4,1,1)
plot(rho_th_c,\'b\',\'linewidth\',1)
ylabel(\'rho\')
subplot(4,1,2)
plot(T_th_c,\'k\',\'linewidth\',1)
ylabel(\'temp\')
subplot(4,1,3)
plot(P_th_c,\'c\',\'linewidth\',1)
ylabel(\'P\')
subplot(4,1,4)
plot(M_th_c,\'r\',\'linewidth\',1)
xlabel(\'number of time steps\')
ylabel(\'mach\')

%residual values
figure(8)
plot(rho_res_c,\'b\',\'linewidth\',1)
hold on
plot(v_res_c,\'r\',\'linewidth\',1)
hold on
plot(T_res_c,\'--k\',\'linewidth\',1)
xlabel(\'number of time steps\')
ylabel(\'residual\')
legend(\'density\',\'velocity\',\'temperature\')
%-------------------------------------------------------------------------%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%comparing  mass flow rate 
%non conservative
figure(9)
plot(x,m,\'r\',\'linewidth\',1.5)
hold on
%conservative form
plot(x,m_c,\'k\',\'linewidth\',1.5)
xlabel(\'axial length of nozzle\')
ylabel(\'mass flow\')
legend(\'non conservative\',\'conservative\')

Matlab function for solving non-conservative form:

function [rho,T,v,A,M,P,m,rho_th,T_th,P_th,M_th,rho_res,T_res,v_res,k_n] = supersonic_non_conserved(n,x,dx,gamma)

% grid space
n=31;
x=linspace(0,3,n);
dx=x(2)-x(1);

% initial values
rho = 1-0.3146*x; % density
T = 1-0.2314*x; % temperature
v = (0.1+1.09*x).*((T).^0.5); % velocity
A = 1+(2.2*((x-1.5).^2)); % area
gamma = 1.4;
throat_A = find(A==1);

%defining cfl and time derivative
c=0.5;  % CFL number
tst = (c*dx)./(sqrt(T)+v);
dt_nc = min(tst);


%time loop
nt_nc = 1400;
error_1 = 9e9;
tol = 1e-9;

k_n =0;

while error_1>tol
    
   k_n=k_n+1;
   


  % copy of initial values
  rho_old = rho;
  v_old = v;
  T_old = T;
  M_old = v./sqrt(T);
  m_old = rho_old.*v_old.*A;
  P_old = rho.*T;
  

  
  % prdictor step
  for j=2:n-1
    dvdx = (v(j+1)-v(j))/dx;
    drhodx = (rho(j+1)-rho(j))/dx;
    dlnAdx = (log(A(j+1))-log(A(j)))/dx;
    dTdx = (T(j+1)-T(j))/dx;
    
    %continuity eqn
    drhodt_p(j) = -rho(j)*dvdx -rho(j)*v(j)*dlnAdx -v(j)*drhodx;
    
    %momentum eqn
    dvdt_p(j) = -v(j)*dvdx -(1/gamma)*(dTdx+(T(j)/rho(j))*drhodx);
    
    %energy eqn
    dTdt_p(j) = -v(j)*dTdx -(gamma-1)*T(j)*(dvdx+v(j)*dlnAdx);
    
    % solution update
    rho(j) = rho(j)+(drhodt_p(j))*dt_nc;
    v(j) = v(j)+(dvdt_p(j))*dt_nc;
    T(j) = T(j)+(dTdt_p(j))*dt_nc;
    
  end
  
  %corrector step
  for j =2:n-1
    dvdx = (v(j)-v(j-1))/dx;
    drhodx = (rho(j)-rho(j-1))/dx;
    dlnAdx = (log(A(j))-log(A(j-1)))/dx;
    dTdx = (T(j)-T(j-1))/dx;
    
    %continuity eqn
    drhodt_c(j) =  -rho(j)*dvdx -rho(j)*v(j)*dlnAdx -v(j)*drhodx;
    
    %momentum eqn
    dvdt_c(j) = -v(j)*dvdx -(1/gamma)*(dTdx+(T(j)/rho(j))*drhodx);
    
    %energy eqn
    dTdt_c(j) = -v(j)*dTdx -(gamma-1)*T(j)*(dvdx+v(j)*dlnAdx);
  end
  
  %average values
  drhodt_avg = 0.5*(drhodt_p+drhodt_c);
  dvdt_avg = 0.5*(dvdt_p+dvdt_c);
  dTdt_avg = 0.5*(dTdt_p+dTdt_c);
  
  %final corrected values
  for i= 2:n-1
    rho(i) = rho_old(i)+(drhodt_avg(i))*dt_nc;
    v(i) = v_old(i)+(dvdt_avg(i))*dt_nc;
    T(i) = T_old(i)+(dTdt_avg(i))*dt_nc;
  end
  
  % boundary condition
  
  %inlet
  v(1) = 2*v(2) -v(3);
  rho(1)=1;
  T(1)=1;
  
  %outlet
  v(n) = 2*v(n-1)-v(n-2);
  rho(n) = 2*rho(n-1)-rho(n-2);
  T(n) = 2*T(n-1)-T(n-2);
 
  %mach number 
  M = v./sqrt(T);
  %non dimensional pressure
  P = rho.*T;
  %mass flow 
  m = rho.*v.*A;
  
  % variables at throat
  rho_th(k_n) = rho(throat_A);
  T_th(k_n) = T(throat_A);
  P_th(k_n) = P(throat_A);
  M_th(k_n) = M(throat_A);
  
  %residual values
  rho_res(k_n) = drhodt_avg(throat_A);
  v_res(k_n) = dvdt_avg(throat_A);
  T_res(k_n) = dTdt_avg(throat_A);
  
  %error of mass flow
  error_1 = max(abs(m-m_old));
  
  %mass flow at different time step
  figure(1)
  hold on
  grid on
  if k_n==1
      plot(x,m,\'b\',\'linewidth\',1)
  else if k_n==50
          plot(x,m,\'c\',\'linewidth\',1)
      else if k_n==100
              plot(x,m,\'k\',\'linewidth\',1)
          else if k_n==150
                  plot(x,m,\'g\',\'linewidth\',1)
              else if k_n==200
                      plot(x,m,\'y\',\'linewidth\',1)
                  else if k_n==700
                          plot(x,m,\'r\',\'linewidth\',1)
                          hold on 
                          xlabel(\'x/Xo\')
                          ylabel(\'mass flow rate\')
                          legend(\'dt=0\',\'dt=50\',\'dt=100\',\'dt=150\',\'dt=200\',\'dt=700\')
                          title(\'mass flow through nozzle at different time step\')
                      end
                  end
              end
          end
      end
  end

  
end



end

Matlab function for solving conservative form:

function  [rho_c,T_c,v_c,A_c,M_c,P_c,m_c,rho_th_c,T_th_c,P_th_c,M_th_c,rho_res_c,T_res_c,v_res_c,k_c] = supersonic_conserved(n,x,dx,gamma)

% grid space
n=31;
x=linspace(0,3,n);
dx=x(2)-x(1);
gamma = 1.4;

%initial condition 

% density and temperature
for i=1:n
    
    
    if (x(i)>=0 & x(i)<=0.5 )
        rho_c(i) = 1;
        T_c(i) = 1;
    else if( x(i)>=0.5 & x(i)<=1.5)
            rho_c(i) = 1-((0.366)*(x(i)-0.5));
            T_c(i) = 1-((0.167)*(x(i)-0.5));
        else if (x(i)>=1.5 & x(i)<=3.5)
                rho_c(i) = 0.634-((0.3879)*(x(i)-1.5));
                T_c(i) = 0.833-((0.3507)*(x(i)-1.5));
            end
        end
    end
    
  
    
end

% area
A_c = 1+(2.2*((x-1.5).^2));
%nozzle throat section
throat_c = find(A_c==1);
% initial velocity
v_c = (0.59)./(rho_c.*A_c);

% initial values of u1 u2 u3
u1 = rho_c.*A_c;
u2 = rho_c.*A_c.*v_c;
u3 = rho_c.*((T_c/(gamma-1))+((gamma/2)*((v_c).^2))).*A_c;

% time step calculation
c = 0.5; %CFL number
%tsp = (c*dx)./(sqrt(T_c)+v_c);
%dt = min(tsp);

% time loop
nt =1400;
error = 9e9;
tol = 1e-9;

k_c = 0;

while error>tol
    k_c=k_c+1;
    m_c_old = rho_c.*v_c.*A_c;
    
%for k=1:nt
    
    u1_old = u1 ;
    u2_old = u2 ;
    u3_old = u3 ;
    
    tsp = (c*dx)./(sqrt(T_c)+v_c);
    dt = min(tsp);
    
    f1 = u2;
    term1 = (((gamma-1)/gamma).*(u3-((gamma/2).*(((u2).^2)./u1))));
    f2 = ((u2.^2)./u1)+ term1;
    f3 = (gamma.*(u2.*u3./u1))-((gamma.*(gamma-1)/2).*(((u2).^3)./((u1).^2)));
    
    
    %predictor step 
    for i=2:n-1
        
        df1dx(i) = ((f1(i+1))-(f1(i)))/dx;
        df2dx(i) = ((f2(i+1))-(f2(i)))/dx;
        df3dx(i) = ((f3(i+1))-(f3(i)))/dx;
        dAdx(i) = ((A_c(i+1))-(A_c(i)))/dx; 
        j2(i) = ((1/gamma)*rho_c(i)*T_c(i)*dAdx(i));
        
        %continuity eqn 
        du1dt_p(i) = (-df1dx(i));
        %momentum eqn
        du2dt_p(i) = (-df2dx(i))+j2(i);
        %energy eqn
        du3dt_p(i) = (-df3dx(i));
        
        %solutioin update 
        u1(i) = u1(i)+(du1dt_p(i)*dt);
        u2(i) = u2(i)+(du2dt_p(i)*dt);
        u3(i) = u3(i)+(du3dt_p(i)*dt);
    end
    
    %updating premitive variables
    rho_c = u1./A_c;
    v_c = u2./u1;
    T_c = (gamma-1)*((u3./u1)-((gamma/2)*((v_c).^2)));
    %P_c = rho_c.*T_c;
    
    %predicted flux values
    f1 = u2;
    term1 = (((gamma-1)/gamma).*(u3-((gamma/2).*(((u2).^2)./u1))));
    f2 = ((u2.^2)./u1)+ term1;
    f3 = (gamma.*(u2.*u3./u1))-((gamma.*(gamma-1)/2).*(((u2).^3)./((u1).^2)));
    
    %corrector step
    for j=2:n-1
        
        df1dx(j) = (f1(j)-f1(j-1))/dx;
        df2dx(j) = (f2(j)-f2(j-1))/dx;
        df3dx(j) = (f3(j)-f3(j-1))/dx;
        dAdx(j) = (A_c(j)-A_c(j-1))/dx; 
        j2(j) = (1/gamma)*rho_c(j)*T_c(j)*dAdx(j);
                     
        %continuity eqn 
        du1dt_c(j) = -df1dx(j);
        %momentum eqn
        du2dt_c(j) = -df2dx(j)+j2(j);
        %energy eqn
        du3dt_c(j) = -df3dx(j);
    end
    
    %average time derivatives
    du1dt_avg = 0.5*(du1dt_p+du1dt_c);
    du2dt_avg = 0.5*(du2dt_p+du2dt_c);
    du3dt_avg = 0.5*(du3dt_p+du3dt_c);
    
    %final corrected values
    for j=2:n-1
        u1(j) = u1_old(j)+(du1dt_avg(j))*dt;
        u2(j) = u2_old(j)+(du2dt_avg(j))*dt;
        u3(j) = u3_old(j)+(du3dt_avg(j))*dt;
    end
    
    %boundary conditions
    
    %inlet
    rho_c(1) = 1;
    T_c(1) = 1;
    u1(1) = rho_c(1)*A_c(1);
    u2(1) = 2*u2(2)-u2(3);
    v_c(1) = u2(1)/u1(1);
    u3(1) = u1(1)*((T_c(1)/(gamma-1))+((gamma/2)*(v_c(1)^2)));
    P_c(1) = rho_c(1)*T_c(1);
    %inlet fluxes
    f1(1) = u2(1);
    f2(1) = (((u2(1)).^2)./u1(1))+((gamma-1/gamma)*(u3(1)-((gamma*(u2(1)).^2)./(2*u1(1)))));
    f3(1) = (gamma*(u2(1).*u3(1)./u1(1)))-((gamma*(gamma-1)/2)*(((u2(1)).^3)./((u1(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);
    %fluxes at outlet
    f1(n) = u2(n);
    f2(n) = (((u2(n))^2)/u1(n))+((gamma-1/gamma)*(u3(n)-((gamma*(u2(n))^2)/(2*u1(n)))));
    f3(n) = (gamma*(u2(n)*u3(n)/u1(n)))-((gamma*(gamma-1)/2)*(((u2(n))^3)/((u1(n))^2)));
    %premitive variables
    rho_c(n) = u1(n)/A_c(n);
    v_c(n) = u2(n)/u1(n);
    T_c(n) = (gamma-1)*((u3(n)/u1(n))-((gamma/2)*((v_c(n))^2)));
    P_c(n) = rho_c(n)*T_c(n);
    
    %mach number 
    M_c = v_c./sqrt(T_c);
    %non dimensional pressure
    P_c = rho_c.*T_c;
    %mass flow 
    m_c = rho_c.*v_c.*A_c;
    
    % variables at throat
    rho_th_c(k_c) = rho_c(throat_c);
    T_th_c(k_c) = T_c(throat_c);
    P_th_c(k_c) = P_c(throat_c);
    M_th_c(k_c) = M_c(throat_c);
    
    %residual at throat
    rho_res_c(k_c) = du1dt_avg(throat_c);
    v_res_c(k_c) = du2dt_avg(throat_c);
    T_res_c(k_c) = du3dt_avg(throat_c);    
 
    
    %error of mass flow
    error = max(abs(m_c-m_c_old));
    
    %mass flow at different time step
    figure(2)
    hold on
    grid on
    if k_c==1
        plot(x,m_c,\'b\',\'linewidth\',1)
        hold on
    else if k_c==100
            plot(x,m_c,\'k\',\'linewidth\',1)
            hold on
        else if k_c==200
                plot(x,m_c,\'y\',\'linewidth\',1)
                hold on
            else if k_c==700
                    plot(x,m_c,\'r\',\'linewidth\',1)
                    hold on 
                    xlabel(\'x/Xo\')
                    ylabel(\'mass flow rate\')
                    legend(\'dt=0\',\'dt=100\',\'dt=200\',\'dt=700\')
                    title(\'mass flow through nozzle at different time step\')
                end
            end
        end
    end

     
     
end



end

In the above program, the function for solving 1D supersonic nozzle flow in non-conservative form and conservative form is written which calculates the premitive variables at different time step and the variation of premitive variables at throat section of the nozzle is observed through the plots and also mass flow at different time step through out the nozzle is observed in the plots.

Number of cycles required for convergence and Time required for solving:

Time take by non-conservative solver = 10.2613s

number of time steps for convergence by non-conservative solver = 1574 time steps

Time taken by conservative solver = 5.5192s

number of time steps for convergence by conservative solver = 2002 time steps

Non-conservative form plots :-

1) steady state profile of non-conservative form

we can observe from the above plot that the variation of density, pressure, temperature and mach number through out the non-dimensional distance of the supersonic nozzle is similar to the analytical solution from the J.D.Anderson book. The mach number increases at the exit of the nozzle, the pressure drops at the exit of the nozzle, the density is high at inlet of the nozzle and drops at exit of nozzle, the temperature also drops at exit of the nozzle.

2) mass flow variation at different time step

The mass flow variation at different time as a function of non-dimensional distance of the nozzle is plotted above. We can see that the mass flow through nozzle at dt=0 is changing largely is not close to the exact solution. The same large variation of mass flow through the nozzle is observed for dt=50, dt=100 and dt=150. But we can observe that the mass flow variation is small and at dt=700 the mass flow does not vary through the nozzle and is constant.

3) Variation with respect to time step at nozzle throat

The variation of density, pressure, temperature and mach number at the throat section of the nozzle is plotted with respect to time step. We can observe that the variation of variables is large at the small time steps and the variation decreases at large time steps and finaly reaches steady state condition.

4) Residual variation with respect to time step

The avereage time derivative is plotted as a function of time steps which is known as residual values. We can see that values of density,temperature and velocity is initialy high and it oscilates upto time step=600 and then reaches steady state and is constant.

Conservative form plots:-

1) steady state profile

The steady state profile of conservative solver and non-conservative solver is similar and it is almost same as the profile given in J.D.Anderson text book.

2) Mass flow variation at different time step

The mass flow variation is large at dt=0 and dt=100 and the variation reduces at dt=200 and finaly the mass flow variation is constant at dt=700.

3) Variation with respect to time steps at nozzle throat

we can see that the variation of density, temperature, pressure and mach number at throat section of the nozzle is fluctuating at initial time steps and after 600 time steps the fluctuation reduces and the values are almost constant.

4) Residual variation with respect to time

The average derivative of density, temperature and velocity are plotted against time steps. We can observe from the plot that the values of variables are high at starting time steps and gradualy the becomes constant at larger time step. When compared to non-conservative form the residuals from conservative form are higher and becomes zero as time step increases.

Mass flow through the nozzle in conservative solver and non-conservative solver:-

The mass flow across the non-dimensional distnace of the nozzle is plotted in above plot. As from the plot we can say that the conservative form solver conserves the mass flow and the non-conservative form does not conserve as the name implies. In non-conservative form solver the mass flow variation across the nozzle is large and in conservative form solver the mass flow variation across the nozzle is small and it is almost close to exact solution.

Conclution:

From the output obtained we can conclude that the conservative form solver takes less time to solve the code with time = 10.26 and the non-conservative form solver takes alittle more time to solve the code withtime = 5.59 and thus the conservative is the fastest solver.

The conservative form solver being the faster solver has more calculation effort than non-conservative form. The residuals from conservative form solver is higher when compared to the non-conservative form. Thus, both conservative and non-conservative solver has its own advantages and disdvantages.