1D FLOW SIMULATION OF SUPERSONIC NOZZLE USING MACORMACK METHOD

ABSTRACT:

This project main objective is the numerical solution of a 1D Supersonic Nozzle Flow using the Macormack Method.

Although analytically identical, the numerical solution of conservation and nonconservation forms of the governing equations is usually different. That is the reason why both ways are solved and compared in this project. For that, the number of iterations and computation time will be compared, as long as the final profile for the normalized mass flow rate. During this project,both a grid dependence test and a CFL-based adaptative time loop control was programmed to assure the fidelity of the numerical results.

The conservation and nonconservation equations are analytically identical. Both come from applying the integral conservation method to the control volume represented by the nozzle, and therefore the analytical solution is identical. However, in numerical analysis, the two sets of equations cant be used indistinguisible, and a deeper analysis must be carried out depending on the problem.

As a general rule when the control volume is fixed in space with the fluid moving through the governing equations will be in conservation form. However, when the control volume is moving with the fluid in a sense that same fluid particles are always remain inside the control volume, the calculation will generally performed in nonconservation form.

 THE MAIN PROGRAM CODE:

clear all
close all
clc

% inputs variables
% Value of 'n' can be varied
n = 31;
L = 3;
x = linspace(0,L,n);

% Courant number
% The smaller-than-necessary value of 'c', helps in increasing the accuracy
c = 0.5;

% To have a node exactly at the throat ,the no.of nodes must be an ODD number
throat = ((n-1)/2)+1;

% Time steps
nt = 1400;

tic
% Function for non conservative form
 [v,rho,t,m,mach,pr,v_t,rho_t,t_t,m_t,mach_t,pr_t] = non_conservative_flow(n,L,c,nt,throat);

% Calculating the simulation time
Non_conservative_simulation_time = toc

tic
% Function for conservative form
 [v2,rho2,t2,m2,mach2,pr2,v_t2,rho_t2,t_t2,m_t2,mach_t2,pr_t2] = conservative_flow(n,L,c,nt,throat);
 
% Calculating the simulation time 
Conservative_simulation_time = toc

 
%% Varaiation of flow parameters at the throat
figure(3)
a1 = subplot(4,1,1);
plot(linspace(1,nt,nt),rho_t); 
hold on
plot(linspace(1,nt,nt),rho_t2,'color','r');
hold on
plot(linspace(nt-400,nt+200,600),linspace(0.634,0.634,600),'color','k');
grid on;
ylabel("rho'",'Fontsize',12,'fontweight','bold');
h = legend(a1,{'Non-Conservative Form', 'Conservative Form','Analytical solution'});
set(h, 'Location', 'northeastoutside')
axis([0 1600]);
title('Timewise variation of flow parameters at throat','Fontsize',12,'color','r','fontweight','bold');

a2 = subplot(4,1,2);
plot(linspace(1,nt,nt),t_t); 
hold on
plot(linspace(1,nt,nt),t_t2,'color','r'); 
hold on
plot(linspace(nt-400,nt+200,600),linspace(0.833,0.833,600),'color','k');
grid on;
ylabel("T'",'Fontsize',12,'fontweight','bold');
h = legend(a2,{'Non-Conservative Form', 'Conservative Form','Analytical solution'});
set(h, 'Location', 'northeastoutside')
axis([0 1600]);


a3 = subplot(4,1,3);
plot(linspace(1,nt,nt),pr_t); 
hold on
plot(linspace(1,nt,nt),pr_t2,'color','r'); 
hold on
plot(linspace(nt-400,nt+200,600),linspace(0.528,0.528,600),'color','k');
grid on;
ylabel("Pr'",'Fontsize',12,'fontweight','bold');
h = legend(a3,{'Non-Conservative Form', 'Conservative Form','Analytical solution'});
set(h, 'Location', 'northeastoutside')
axis([0 1600]);


a4 = subplot(4,1,4);
plot(linspace(1,nt,nt),mach_t); 
hold on
plot(linspace(1,nt,nt),mach_t2,'color','r'); 
hold on
plot(linspace(nt-400,nt+200,600),linspace(1.0,1.0,600),'color','k');
grid on;
ylabel("Ma'",'Fontsize',12,'fontweight','bold');
xlabel('No of time steps','fontweight','bold');
h = legend(a4,{'Non-Conservative Form', 'Conservative Form','Analytical solution'});
set(h, 'Location', 'northeastoutside')
axis([0 1600]);




% Variation in flow parameters along x axis
figure(4)
p1 = subplot(4,1,1);
plot(x,mach); 
hold on
plot(x,mach2,'color','r','*', 'markersize', 3);
grid on;
ylabel("Ma'",'Fontsize',12,'fontweight','bold');
h = legend(p1,{'Non-Conservative Form', 'Conservative Form'});
set(h, 'Location', 'northeastoutside')
title('Variation of Flow parameters along x axis','Fontsize',12,'color','r','fontweight','bold');


p2 = subplot(4,1,2);
plot(x,pr);
hold on
plot(x,pr2,'color','r','*', 'markersize', 3); 
grid on;
ylabel("Pr'",'Fontsize',12,'fontweight','bold');
h = legend(p2,{'Non-Conservative Form', 'Conservative Form'});
set(h, 'Location', 'northeastoutside')


p3 = subplot(4,1,3);
plot(x,rho); 
hold on
plot(x,rho2,'color','r','*', 'markersize', 3);
grid on;
ylabel("rho'",'Fontsize',12,'fontweight','bold');
h = legend(p3,{'Non-Conservative Form', 'Conservative Form'});
set(h, 'Location', 'northeastoutside')


p4 = subplot(4,1,4);
plot(x,t); 
hold on
plot(x,t2,'color','r','*', 'markersize', 3);
grid on;
ylabel("T'",'Fontsize',12,'fontweight','bold');
xlabel('Non-dimensional length','fontweight','bold');
h = legend(p4,{'Non-Conservative Form', 'Conservative Form'});
set(h, 'Location', 'northeastoutside')

% mass flow rate comparison between two forms of equation
figure(5)
plot(x,m,'r', 'linewidth', 1.5)
hold on
plot(x,m2,'b', 'linewidth', 1.5)
hold on
plot(x,linspace(0.579,0.579,n),'color','k','--', 'markersize', 3)
ylabel('Non-dimensional mass flow rate','fontweight','bold','fontsize',12);
xlabel('Non-dimensional distance through nozzle(x)','Fontsize',12,'fontweight','bold');
legend({'Non-Conservative Form', 'Conservative Form','Exact analytical solution'});
title('Non-dimensional mass flow rate comparison','fontweight','bold','fontsize',12,'color','r'); 
axis([0 3 0.575 0.600]);
grid on

THE NON CONSEVATIVE FORM FUNCTION CODE:

function [v,rho,t,m,mach,pr,v_t,rho_t,t_t,m_t,mach_t,pr_t] = non_conservative_flow(n,L,c,nt,throat)

% inputs variables
n = 31;
L = 3;
x = linspace(0,L,n);
dx = x(2)-x(1);

% Initial conditions of the nondimensional numbers
% x,rho,v,t,a are nondimensional numbers
rho = 1-0.3146*x;
t = 1-0.2314*x;
v = (0.1+1.09*x).*t.^0.5;
a = 1+2.2*(x-1.5).^2;
gamma = 1.4;
c = 0.5;
nt = 1400;
throat = ((n-1)/2)+1;
% calculating the value of time step based on courant number criteria
for i = 1:n
  del_t(i) = c*(dx/(t(i)^0.5+v(i)));  
end
  dt = min(del_t);


% Time loop
for k = 1:nt

rho_old = rho;
v_old = v;
t_old = t;

  % 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 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 


% average derivative
drhodt = 0.5*(drhodt_p+drhodt_c);
dvdt = 0.5*(dvdt_p+dvdt_c);
dtdt = 0.5*(dtdt_p+dtdt_c);

% final solution update
for p = 2:n-1
    rho(p) = rho_old(p) + drhodt(p)*dt;
    v(p) = v_old(p) + dvdt(p)*dt;
    t(p) = t_old(p) + dtdt(p)*dt;
    
end
   
    
% Apply boundary conditions

  % inlet
  v(1)=2*v(2)-v(3);
  
  % 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);
  
  % calculating solutions 
   m = rho.*v.*a;
   mach = v./((t).^0.5);
   pr = rho.*t;
  
  % Throat values 
   rho_t(k) = rho(throat);
   t_t(k) = t(throat);
   pr_t(k) = pr(throat);
   mach_t(k) = mach(throat);
   v_t(k) = v(throat);
   m_t(k) = m(throat);
   
   
    % plotting mass flow rate at different time steps
     figure(1);
   if k  == 50
      plot(x, m, 'b', 'linewidth', 1.5)
    hold on
     elseif k == 100
      plot(x, m, 'r', 'linewidth', 1.5)
    hold on
     elseif k == 150
      plot(x, m, 'c', 'linewidth', 1.5)
    hold on
     elseif k == 200
      plot(x, m, 'g', 'linewidth', 1.5)
    hold on
     elseif k == 550
      plot(x, m, 'k', 'linewidth', 1.5) 
    hold on 
     elseif k == 700
      plot(x, m, 'y', 'linewidth', 1.5)
    hold on
      elseif k == 1400
      plot(x, m,'m', 'linewidth', 1.5)
    
      grid on
      ylabel('Non-dimensional mass flow rate','fontweight','bold','fontsize',12);
      xlabel('Non-dimensional distance through nozzle(x)','Fontsize',12,'fontweight','bold');
      title('Non-dimensional mass flow rate(Non-conservative form) along x, at different time steps','Fontsize',12,'fontweight','bold','color','r');
      legend({'50 time steps','100 time steps','150 time steps','200 time steps','550 time steps','700 time steps','1400 time steps'});
   end
  
  end
  
end

THE CONSERVATIVE FORM FUNCTION CODE:

function [v,rho,t,m,mach,pr,v_t,rho_t,t_t,m_t,mach_t,pr_t] = conservative_flow(n,L,c,nt,throat)


% input variables
n = 31;
L = 3;
x = linspace(0,L,n);
dx = x(2)-x(1);


% Initial conditions of the nondimensional numbers
% x,rho,v,t,a are nondimensional numbers
a = 1+2.2*(x - 1.5).^2;

for i = 1:n
    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

% Assuming steady state flow for the initial condition
  v = 0.59./(rho.*a);
  gamma = 1.4;
  c = 0.5;
% Solution vector parameters
  pr = rho.*t;
  u1 = rho.*a;
  u2 = rho.*a.*v;
  u3 = u1.*((t./(gamma-1))+((gamma/2).*(v.^2)));

% calculating the value of time step based on courant number criteria
% Time step is not allowed to change for every iteration, to avoid "Time warp"
     dt = min (c.*(dx./((t.^0.5)+v))); 
     nt = 1400; 
 % Time loop
for k = 1:nt
   
    % Copy of solution vectors   
    u1_old = u1;
    u2_old = u2;
    u3_old = u3;
    
    
    % "Pure" form of the flux terms
    f1 = u2;
    term = (((gamma-1)/gamma)*(u3-((gamma/2)*((u2.^2)./u1))));
    f2 = ((u2.^2)./u1)+ term;
    f3 = ((gamma*u2.*u3)./u1)-(gamma*(gamma-1)*(u2.^3))./(2*u1.^2);


    % Predictor Method
    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;
        dlogadx(i) = (log(a(i+1))-log(a(i)))/dx;
        
        % j_term can also be used instead of j2. Infact j_term is the exact "Pure" form
        % j_term(i) = term(i)*dlogadx(i); 
        
        dadx(i) = (a(i+1)-a(i))/dx;
        j2(i) = (1/gamma)*(rho(i)*t(i))*dadx(i); 
        
        % continuity equation
        du1dt_p(i) = -df1dx(i);
        
        % momentum equation
        du2dt_p(i) = -df2dx(i) + j2(i);
        
        % energy equation
        du3dt_p(i) = -df3dx(i);
        
        %Solution 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
    
    % Decoding primitive terms to initialize the corrector step
    rho = u1./a;
    v = u2./u1;
    t = (gamma-1)*((u3./u1)-((gamma*v.^2)/2));
    
    % "Pure" form of the flux terms for the corrector step
    f1 = u2;
    term1 = (((gamma-1)/gamma)*(u3-(gamma*u2.^2)./(2*u1)));
    f2 = ((u2.^2)./u1)+ term1;
    f3 = ((gamma*u2.*u3)./u1)-(gamma*(gamma-1)*(u2.^3))./(2*u1.^2);
    
    % Corrector Method
    for i = 2:n-1

        df1dx(i) = (f1(i)-f1(i-1))/dx;
        df2dx(i) = (f2(i)-f2(i-1))/dx;
        df3dx(i) = (f3(i)-f3(i-1))/dx;
        dlogadx(i) = (log(a(i))-log(a(i-1)))/dx;
        
        % j_term1 can also be used instead of j3. Infact j_term1 is the exact "Pure" form 
        % j_term1(i) = term1(i)*dlogadx(i); 
        
        dadx1(i) = (a(i)-a(i-1))/dx;
        j3(i) = (1/gamma)*(rho(i)*t(i))*dadx1(i); 
 
        % Continuity equation
        du1dt_c(i) = -df1dx(i);
        
        % Momentum equation
        du2dt_c(i) = -df2dx(i) + j3(i);
        
        %Energy equation
        du3dt_c(i) = -df3dx(i);
        
    end
    
    % Average derivative
    
    du1dt = 0.5*(du1dt_p + du1dt_c);
    du2dt = 0.5*(du2dt_p + du2dt_c);
    du3dt = 0.5*(du3dt_p + du3dt_c);
    
    % Final solution update
 for j = 2:n-1
    u1(j) = u1_old(j) + du1dt(j)*dt;
    u2(j) = u2_old(j) + du2dt(j)*dt;
    u3(j) = u3_old(j) + du3dt(j)*dt;
 end
    
% Apply the boundary conditions
    
  % inlet
  
    % u1 value remains same, because rho is not a floating variable and also 'a' doesnt change with time
    u1(1) = rho(1)*a(1);
    
    % velocity is the floating primitive term, which can be calculated by extrapolating u2
    u2(1) = 2*u2(2) - u2(3);
    
    % Floating variable v, is calculated for every time step
    v(1) = u2(1)/u1(1);
    
    % updating u3
    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);
    
   % Decoding primitives terms 
   rho = u1./a; 
   v = u2./u1;
   t = (gamma-1)*((u3./u1)-((gamma*v.^2)/2));
    
   
  % calculating solutions 
   m = rho.*v.*a;
   mach = v./((t).^0.5);
   pr = rho.*t;
   throat = ((n-1)/2)+1; 
   % Throat values 
   rho_t(k) = rho(throat);
   t_t(k) = t(throat); 
   pr_t(k) = pr(throat);
   mach_t(k) = mach(throat);
   v_t(k) = v(throat);
   m_t(k) = m(throat);
   
   
   % plotting mass flow rate at different time steps
    figure(2);
   if k  == 50
      plot(x, m, 'b', 'linewidth', 1.5)
    hold on
     elseif k  == 100
      plot(x, m, 'r', 'linewidth', 1.5)
    hold on
     elseif k  == 150
      plot(x, m, 'c', 'linewidth', 1.5)
    hold on
     elseif k  == 200
      plot(x, m, 'g', 'linewidth', 1.5)
    hold on
   elseif k == 550
      plot(x, m, 'y', 'linewidth', 1.5) 
    hold on  
     elseif k  == 700
      plot(x, m, 'k', 'linewidth', 1.5)
    hold on
      elseif k  == 1400
      plot(x, m,'m', 'linewidth', 1.5)
      
      grid on
      ylabel('Non-dimensional mass flow rate','fontweight','bold','fontsize',12);
      xlabel('Non-dimensional distance through nozzle(x)','Fontsize',12,'fontweight','bold');
      title('Non-dimensional mass flow rate(Conservative form) along x, at different time steps','Fontsize',12,'fontweight','bold','color','r');
      legend({'50 time steps','100 time steps','150 time steps','200 time steps','550 time steps','700 time steps','1400 time steps'});
   end
 end

 end

RESULTS:

The comparison of flow characteristics of the Quasi 1D supersonic nozzle;

Conservative form vs Non Conservative form;

1. Variations of flow parameters with respect to TIME,

2. Transient Non Dimensional MASS FLOW Rates (Non Conservative Form)

3. Variation of FLOW PARAMETERS along X axis;

4. Transinet Non Dimensional MASS FLOW Rate (Conservative Form)

5. Non Dimensional Mass Flow Rate Conservative Form vs Non Conservative Form

6. Simulation Time Comparision

The computational time for both the forms are above mentioned.

OBSERVATIONS:

1. It is clear that both the conservative and the non-conservative form return identical profiles.

2. Variation of mass flow rate along the length of nozzle, we can conclude that the plot the conservative form gives a mass flow rate variation that is close to a constant value and more accurate results at the steady state. The non-conservative form tends to show some oscillations at the inflow and outflow boundaries and is less satisfactory compared to the conservation form.

3. Although the non-conservation form and the conservation form have different initial conditions, it can be observed that variations in the mass flow rates are far less in the conservative form and closer to the exact analytical value.

4. We can clearly observe that the conservative form converges quicker and has lot less variations initially compared to non -conservative form. After convergence, it can be observed that non-conservative form returns more accurate results compared to conservative form for the given grid size.

5. Although convervative form requires more effort to code compared to non-conservative form. Once developed ,computational time for conservative form should be much lesser than the non-conservative form.