Simulation of a 1D Super-sonic nozzle flow simulation using Macormack Method

Objective: Simulation of 1D supersonic nozzle flow simulation using the Macormack Method using Non-conservative and conservative forms of equations.

Non-conservative Flow:

In this non-conservative flow, a small control volume is taken inside the pipe flow but with a small change, the control volume here is moving with the instead of being fixed as in case of conservative flow.

McCormack method is used for solving the governing equations. This concept is used in programming to give the exact solution:

(A) Governing Equations:

(i) Continuity Equation:

`{delrho\'}/{delt\'} = -rho\'*{delv\'}/{delx\'} - rho\'*v\'*{del(ln A\')}/{delx\'} - v\'*{delrho\'}/{delx\'}`

(ii) Momentum Equation:

`{delv\'}/{delt\'} = - v\'*{delv\'}/{delx\'} -1/gamma*({delT\'}/{delx\'} + {T\'}/{rho\'}*{delrho\'}/{delx\'})`

(iii) Energy Equation:

`{delT\'}/{delt\'} = - v\'*{delT\'}/{delx\'} -(gamma-1)*T\'*({delv\'}/{delx\'} + v\'*{del(lnA\')}/{delx\'})`

(B) Initial conditions of Profile at time, t =0:

Density, `rho = 1-0.3146*x;`

Temperature, `T = 1- 0.2314*x`

velocity, `v = (0.1 + 1.09*x)*t^0.5`

Area,`A = 1 + 2.2*(x-1.5)^2`

(C) Boundary conditions

Inlet:

`v(1) = 2v(2) - v(3)`

Outlet:

`v(n) = 2v(n-1)-v(n-2)`

`rho(n) = 2rho(n-1)-rho(n-2)`

`T(n) = 2T(n-1)-T(n-2)`

conservative Flow:

In this conservative flow, when a fluid is flowing through a pipe a small control volume is taken in some area inside that pipe and this control volume is fixed, so that the equation can be solved find out the properties of the flow through that point.

A) Governing Equations:

(i) Continuity Equation:

`{delU_1}/{delt\'} = -{delF_1}/{delx\'}`

(ii) Momentum Equation:

`{delU_2}/{delt\'} = -{delF_2}/{delx\'} + J_2`

(iii) Energy Equation:

`{delU_3}/{delt\'} = -{delF_3}/{delx\'}`

Solution vector:U

Flux vector:F

Source term:`J_2`

`U_1= rho\'A\'`

`U_2=rho\'A\'v\'`

`U_3=rho\'({e\'}/{gamma-1} + gamma/2v\'^2)A\'`

`F_1= rho\'A\'v\'`

`F_2= rho\'A\'v\' +(1/gamma)p\'A\'`

`F_3=rho\'({e\'}/{gamma-1} + gamma/2v\'^2)v\'A\' +p\'A\'v\'`

`J_2 = 1/gamma*p\'*{dA\'}/{dx\'}`

(B) Initial conditions of Profile at time, t=0:

`x(i)>=0 && x(i)<=0.5`

`rho(i) = 1`

`T(i) = 1`

`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);`

`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)`

Area,`A = 1 + 2.2*(x-1.5)^2`

Assuming steady state flow

`v = 0.59/(rho*A)`

(C) Boundary conditions

Inlet:

`U_1(1)=rho(1)*A(1)`

`U_2(1)= 2U_2(2)-U_2(3)`

`v(1)= {U_2(1)}/{U_1(1)}`

`T(1) = 1`

`U_3(1) = U_1(1)*(({T(1)}/(gamma-1))+(gamma/2)*(v(1)^2))`

Outlet:

`U_1(n) = 2U_1(n-1) - U_1(n-2)`

`U_2(n) = 2U_2(n-1) - U_2(n-2)`

`U_3(n) = 2U_3(n-1) - U_3(n-2)`

The flow field are Non-Dimensionalized, where the flow variable are referenced to their reservoir values.

Non-Dimensional Temperature:`T\'=T/T_o` where `T_o`= reservoir Temperature

Non-Dimensional Density:`rho\' = rho/rho_o`where `rho_o`= reservoir Density

Non-Dimensional length:`x\'=x/L` where L = length of nozzle

Non-Dimensional velocity:`v\'= v/a_o` where `a_o` = speed of sound

Non-Dimensional Time:`t\'=t/(L/a_o)`

Non-Dimensional Area`A\'=A/(A^*)` where `A^*`= Throat area

Calling a Non-conservative function:

function [v,rho,t,m,p,mach,v_t,rho_t,t_t,m_t,p_t,mach_t,k,time_execution_1] = non_conservative(n,x,dx,c,gamma,a,nt,throat)

tic
% Inputs
        % Length of the domain
               L = 3
        % no of the grid points
               n = 31
        % x array along the length of the domain 
              x =linspace(0,L,n)
        % Grid size
              dx = x(2) - x(1)
        % other parameters   
              gamma = 1.4
% calculate Initial profile, which are non dimensional 
     rho = 1-0.3146*x;
     t = 1- 0.2314*x;   % t =temperature
     v = (0.1 + 1.09*x).*t.^0.5;
     a = 1 + 2.2*(x-1.5).^2; % Area, a

% Time steps 
     nt = 1600
% calculating the value of time step using courant criteria
for i = 1:n
    del_t(i)= c*(dx/(t(i)^0.5+v(i)));
end
    dt = min(del_t(i));
  
% outer time loop
for  k = 1:nt
    rho_old = rho;
    v_old = v;
    t_old = t;
    a_old = a;
        
% 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
    
    % 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)
    
    % Final solution update
    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
    
        % 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
        % Mass flow rate
             m = rho.*a.*v;
        % pressure
             p = rho.*t;
        % Mach number
             mach = v./((t).^0.5);
             
   % Throat values
           rho_t(k) = rho(throat);
           v_t(k)=v(throat);
           t_t(k)= t(throat);
           p_t(k)=p(throat);
           m_t(k)=m(throat);
           mach_t(k)= mach(throat);
             
   % plotting mass flow rate at different time steps
       figure(1);
       if k == 50
          plot(x,m,\'b\',\'linewidth\',2)
       hold on
       else
           if k == 100
              plot(x,m,\'r\',\'linewidth\',2)
           else
               if k == 200
                  plot(x,m,\'g\',\'linewidth\',2)
               else
                   if k == 400
                       plot(x,m,\'y\',\'linewidth\',2)
                   else
                       if k == 800
                          plot(x,m,\'c\',\'linewidth\',2)
                       else
                           if k == 1600
                               plot(x,m,\'m\',\'linewidth\',2)
                           end
                       end
                   end
               end
           end
       end 
       grid on
       
       xlabel(\'Non-Dimensional distance through nozzle(x)\',\'fontweight\',\'bold\',\'fontsize\',12)
       ylabel(\'Non-Dimensional mass flow rate through nozzle(m)\',\'fontweight\',\'bold\',\'fontsize\',12)
       title(\'Non-Dimensional mass flow rate(Non-conservative form) along x at different time steps\',\'fontweight\',\'bold\',\'fontsize\',12,\'color\',\'r\')
       legend({\'50 time steps\',\'100 time steps\',\'200 time steps\',\'400 time steps\',\'800 time steps\',\'1600 time steps\'})
       
       time_execution_1=toc
end

Calling a conservative function:

function [v,rho,t,m,p,mach,v_t,rho_t,t_t,m_t,p_t,mach_t,k,time_execution_2] = conservative(n,x,dx,c,gamma,a,nt,throat)

tic
% Inputs
        % Length of the domain
               L = 3
        % no of the grid points
               n = 31
        % x array along the length of the domain 
              x =linspace(0,L,n)
        % Grid size
              dx = x(2) - x(1)
        % other parameters
              gamma = 1.4
% calculate Initial profile, which are non dimensional 
        for i = 1:n
            if x(i)>=0 && x(i)<=0.5 
               rho(i) = 1;
               t(i) = 1;
            else
                if 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);
                else        
                    if 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
            end
        end 
                
 % Assuming steady state flow for the initial profile
           v = 0.59./(rho.*a);
 
 % Time steps 
           nt = 1600
 
 % calculating the value of time step using courant criteria
for i = 1:n
    del_t(i)= c*(dx/(t(i)^0.5+v(i)));
end
    dt = min(del_t(i));
      
  % solution vector parameters         
           p = rho.*t;
           U1 = rho.*a;
           U2 = rho.*a.*v;
           U3 = rho.*((t./(gamma-1))+((gamma/2).*(v.^2))).*a
          
 % outer Time loop 
      for k = 1:nt;
          U1_old = U1;
          U2_old = U2;
          U3_old = U3;
 % Defining flux variables
          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
                    

       dF1dx = (F1(j+1)-F1(j))/dx;       
       dF2dx = (F2(j+1)-F2(j))/dx;
       dF3dx = (F3(j+1)-F3(j))/dx;
       dlogadx = (log(a(j+1))-log(a(j)))/dx
       dadx = (a(j+1)-a(j))/dx  
       J2 = (1/gamma).*rho(j).*t(j).*dadx
       
       % continuity equation
       dU1dt_p(j) = -dF1dx
       % Momentum equation
       dU2dt_p(j) = -dF2dx + J2
       % Energy equation
       dU3dt_p(j) = -dF3dx
       
       % 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
       % Decoding primitive terms to initialize the corrector step
       rho = U1./a;
       v = U2./U1;
       t = (gamma-1)*((U3./U1)-((gamma/2)*((U2./U1).^2)));
       p = rho.*t;
       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
      dF1dx = (F1(j)-F1(j-1))/dx;       
       dF2dx = (F2(j)-F2(j-1))/dx;
       dF3dx = (F3(j)-F3(j-1))/dx;
       dlogadx = (log(a(j))-log(a(j-1)))/dx
       dadx = (a(j)-a(j-1))/dx  
       J2 = (1/gamma).*rho(j).*t(j).*dadx
       
       % continuity equation
       dU1dt_c(j) = -dF1dx
       % Momentum equation
       dU2dt_c(j) = -dF2dx + J2
       % Energy equation
       dU3dt_c(j) = -dF3dx
 end
 
 % Average derivatives
      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 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
 
 % Apply the Boundary conditions
 % Inlet
      U1(1)=rho(1)*a(1);
      U2(1)= 2*U2(2)-U2(3);
      v(1)= U2(1)/U1(1);
      t(1) = 1;
      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);
      
      rho = U1./a;
      v = U2./U1;
      t = (gamma-1)*((U3./U1) - ((gamma/2)*((U2./U1).^2)))
      
 % caculating solutions
       % Mass flow rate    
              m = rho.*v.*a;
       % pressure
             p = rho.*t;
       % Mach number
             mach = v./((t).^0.5);
 % Throat values
           rho_t(k) = rho(throat);
           v_t(k)=v(throat);
           t_t(k)= t(throat);
           p_t(k)=p(throat);
           m_t(k)=m(throat);
           mach_t(k)= mach(throat);
% plotting mass flow rate at different time steps
       figure(2);
       if k == 50
          plot(x,m,\'b\',\'linewidth\',2)
       hold on
       else
           if k == 100
              plot(x,m,\'r\',\'linewidth\',2)
           else
               if k == 200
                  plot(x,m,\'g\',\'linewidth\',2)
               else
                   if k == 400
                       plot(x,m,\'y\',\'linewidth\',2)
                   else
                       if k == 800
                          plot(x,m,\'c\',\'linewidth\',2)
                       else
                           if k == 1600
                               plot(x,m,\'m\',\'linewidth\',2)
                           end
                       end
                   end
               end
           end
       end 
       grid on
       
       xlabel(\'Non-Dimensional distance through nozzle(x)\',\'fontweight\',\'bold\',\'fontsize\',12)
       ylabel(\'Non-Dimensional mass flow rate through nozzle(m)\',\'fontweight\',\'bold\',\'fontsize\',12)
       title(\'Non-Dimensional mass flow rate(conservative form) along x at different time steps\',\'fontweight\',\'bold\',\'fontsize\',12,\'color\',\'r\')
       legend({\'50 time steps\',\'100 time steps\',\'200 time steps\',\'400 time steps\',\'800 time steps\',\'1600 time steps\'})
       
       time_execution_2=toc
end

Main Matlab Program:

clear all
close all
clc

% Inputs
        % Length of the domain
                 L = 3;
        % Number of the grid points
                 n = 31;
        % x array along the length of the domain      
                 x = linspace(0,L,n);
        % Grid spacing
                 dx = x(2) - x(1);
        % Area,a
                 a = 1+2.2*(x-1.5).^2;
                 gamma = 1.4
        % Exact node at the throat, numbers of nodes is an odd number
                 throat = ((n-1)/2)+1;
         
        % courant number
                 c = 0.5
                 
          % No. of time steps 
          nt = 1600;
          
        

% non conservative form
[v1,rho1,t1,m1,p1,mach1,v_t1,rho_t1,t_t1,m_t1,p_t1,mach_t1,k1,time_execution_1]=non_conservative(n,x,dx,c,gamma,a,nt,throat);
% conservative form
[v2,rho2,t2,m2,p2,mach2,v_t2,rho_t2,t_t2,m_t2,p_t2,mach_t2,k2,time_execution_2]=conservative(n,x,dx,c,gamma,a,nt,throat);

% plots comparing nonconservative and conservative
figure(3)
s1=subplot(4,1,1);
plot(linspace(1,nt,nt), rho_t1);
hold on
plot(linspace(1,nt,nt), rho_t2);
hold on
plot(linspace(nt-600,nt+200,800),linspace(0.634,0.634,800),\'color\',\'g\');
grid on
ylabel(\"rho\'\",\'fontsize\',12,\'fontweight\',\'bold\')
h=legend(s1,{\'Non-conservative Form\', \'conservative Form\',\'Analytical solution\'})
set(h, \'Location\',\'northeastoutside\')
title(\'Timewise variation of flow parameters at throat\',\'fontsize\',12,\'fontweight\',\'bold\')

s2=subplot(4,1,2)
plot(linspace(1,nt,nt), t_t1);
hold on
plot(linspace(1,nt,nt), t_t2);
hold on
plot(linspace(nt-600,nt+200,800),linspace(0.833,0.833,800),\'color\',\'g\');
grid on
ylabel(\"T\'\",\'fontsize\',12,\'fontweight\',\'bold\')
h=legend(s2,{\'Non-conservative Form\', \'conservative Form\',\'Analytical solution\'})
set(h, \'Location\',\'northeastoutside\')


s3=subplot(4,1,3)
plot(linspace(1,nt,nt), p_t1);
hold on
plot(linspace(1,nt,nt), p_t2);
hold on
plot(linspace(nt-600,nt+200,800),linspace(0.528,0.528,800),\'color\',\'g\');
grid on
ylabel(\"pr\'\",\'fontsize\',12,\'fontweight\',\'bold\')
h=legend(s3,{\'Non-conservative Form\', \'conservative Form\',\'Analytical solution\'})
set(h, \'Location\',\'northeastoutside\')
legend(\'Non-conservative Form\', \'conservative Form\',\'Analytical solution\')


s4=subplot(4,1,4)
plot(linspace(1,nt,nt), mach_t1);
hold on
plot(linspace(1,nt,nt), mach_t2);
hold on
plot(linspace(nt-600,nt+200,800),linspace(1,1,800),\'color\',\'g\');
grid on
xlabel(\"Number of time steps\",\'fontsize\',12,\'fontweight\',\'bold\')
ylabel(\"Ma\'\",\'fontsize\',12,\'fontweight\',\'bold\')
h=legend(s4,{\'Non-conservative Form\', \'conservative Form\',\'Analytical solution\'})
set(h, \'Location\',\'northeastoutside\')


% variation in flow parameters along x axis
figure(4)
q1=subplot(4,1,1)
plot(x,rho1,\'r\')
hold on
plot(x,rho2,\'b\')
grid on
ylabel(\"rho\'\",\'fontsize\',12,\'fontweight\',\'bold\')
h = legend(q1,{\'Non conservative Form\',\'conservative Form\'});
set(h,\'location\',\'northeastoutside\')
title(\'Variation of flow parameters along x axis\',\'fontsize\',12,\'fontweight\',\'bold\')

q2=subplot(4,1,2)
plot(x,t1,\'r\')
hold on
plot(x,t2,\'b\')
grid on
ylabel(\"T\'\",\'fontsize\',12,\'fontweight\',\'bold\')
h = legend(q2,{\'Non conservative Form\',\'conservative Form\'});
set(h,\'location\',\'northeastoutside\')

q3=subplot(4,1,3)
plot(x,p1,\'r\')
hold on
plot(x,p2,\'b\')
grid on
ylabel(\"pr\'\",\'fontsize\',12,\'fontweight\',\'bold\')
h = legend(q3,{\'Non conservative Form\',\'conservative Form\'});
set(h,\'location\',\'northeastoutside\')

q4=subplot(4,1,4)
plot(x,mach1,\'r\')
hold on
plot(x,mach2,\'b\')
grid on
xlabel(\'Non-Dimensional Length\',\'fontsize\',12,\'fontweight\',\'bold\')
ylabel(\"Ma\'\",\'fontsize\',12,\'fontweight\',\'bold\')
h = legend(q4,{\'Non conservative Form\',\'conservative Form\'});
set(h,\'location\',\'northeastoutside\')

Results and conclusion:

Plot of Variation of flow parameters along x-axis:

From the plots it is clear that the both the conservative and non-conservative form produce same plot after steady state, 1600 time steps are carried out by time marching approach, at which the steady state is almost obtained and the values of primitive variables are good in approximation with exact analytical solution.

Plot of Time variation of flow parameters at throat:

• From the plot it is clear that the non-conservative form lines are closer to the exact values of rho\',T\',pr\', and Mach\' when compared to conservative form of solution
• The Residue (difference between the numerical and exact solution) exists after 1600 time steps.
• Non-conservative form leads to smaller residue, and the quality index is determined by the rate of decay of the residue. So, non-conservative form obtains a better quality index when compared to the conservative form.
• From the plots, it is clear that the steady state is obtained with a good level of approximation after 400 time steps.

Plot of Non-Dimensional mass flow rate(Non-conservative form):

• The mass flowrate obtained for the last three time steps (i.e., 400, 800, 1600) are almost overlying in the same region of space
• This plot substantiates our previous assumption from the observation of \"Time variation of flow parameters at throat\", i.e., after 400 time steps almost steady state is obtained .

Plot of Non-Dimensional mass flow rate(conservative form):

• Even for this conservative form the mass flowrate obtained for the last three time steps (i.e., 400, 800, 1600) are almost overlying in the same region of space.
• The Initial conditions for the conservative form and non conservative form are different so the comparison is inappropriate.
• The steady state is attained after infinite steps, but from the plot after 400 the variation is too small with the exact solution.

Grid dependence:

• If the number of grid points is doubled, and there is only a marginal difference in the solution, then the calculations with the initial grid points are independent.
• In this case n =31, when modifying the value of \'n\' the number of time steps should also be modified, so that the elapsed non dimensional time is same for both cases.
• Grid independence test can also be carried out by comparing the values obtained by conservative and non-conservative forms at the throat, for n=31 and 61. It has to be noted that the non-conservative form results are more accurate when compared with the conservative form results.
• After doubling the number of grid points, the solutions of the conservative form are slightly closer to the exact value, but still attained the same accuracy of non-conservative form with half as many grid points i.e., n = 31. Therefore for practical purposes, n =31 can be treated as grid-independent.

Flow parameters	Non conservative form	Conservative form	Exact value
n = 31	n = 61	n = 31	n = 61
rho\'	0.639	0.638	0.652	0.639	0.634
t\'	0.836	0.835	0.841	0.836	0.833
Pr\'	0.534	0.532	0.548	0.533	0.528
mach\'	0.999	1.0	0.978	0.998	1.00

WHICH METHOD IS FASTER?

For n = 31, c = 0.5

Non-Conservative simulation Time = 93.448

Conservative simulation Time = 171.8982

For n = 61, c = 0.5

Non-Conservative simulation Time = 170.474

Conservative simulation Time = 737.3677