Simulation of a 1D Super-sonic nozzle flow simulation Macormack Method for Conservative and Non - conservative governing equations

Types of Governing Equations :

• Conservative equation
• Non - conservative equation

Conservative Equation (Continuity equation, Momentum equation, Energy equation)

Continuity Equation:

`(del(rho'A'))/(delt') + (del(rho'A'V'))/(delx')`

Momentum Equation:

`(del(rho'A'V'))/(delt') + (del[rho'A'V'^2 + (1/gamma)p'A'))/(delx') = (1/gamma)p'(delA')/(delx')`

Energy Equation:

`(del[rho'((e')/(gamma - 1) + gamma/2V'^2)A'])/(delt') + (del[rho'((e')/(gamma - 1) + gamma/2V'^2)V'A'+p'A'V'])/(delx') = 0`

Non - conservative Equation (Continuity equation, Momentum equation, Energy equation)

Continuity Equation:

`(delrho')/(delt') = -rho'(delV')/(delx') - rho'V'(del(lnA'))/(delx') - V'(delrho')/(delx')`

Momentum Equation:

`(delV')/(delt') = -V'(delV')/(delx') - 1/gamma ((delT')/(delx') + (T')/(rho')(delrho')/(delx'))`

Energy Equation:

`(delT')/(delt') = -V'(delT')/(delx') -(gamma-1)T'[(delV')/(delx') + V'(del(lnA'))/(delx')]`

where,

`x' = x/L`

`V' = V/(a_0)`

`t' = t/(L/a_0)`

`A' = A/(A^*)`

MacCormak's Method:

Considering continuity equation for non-conservative form

Predictor step:

`((delrho_i)/(delt))_p=((rho_(i)^(N+1)) - (rho_i^(N) ))/(trianglet) =rho_i^(N)[((V_(i+1)^N)-(V_i^N))/(trianglex)] - rho_i^NV_i^N[(ln(A_(i+1)) - ln(A_(i)))/(trianglex)]-V_i^N[((rho_(i+1)^N)-(rho_i^N))/(trianglex)]`

`rho_i^(updated) = rho_i^(1) + (delrho)/(delt_(i-p))*dt`

Corrector step:

`((delrhoi)/(delt))_c = rho_i^(u)[((V_(i)^u)-(V_(i-1)^u))/(trianglex)] - rho_i^uV_i^u[(ln(A_(i)) - ln(A_(i-1)))/(trianglex)]-V_i^u[((rho_(i)^u)-(rho_(i-1)^u))/(trianglex)]`

`rho_i^(updated) = rho_i^(1) + (delrho)/(delt_(i-c))*dt`

`((delV)/(delt))_(avg) =0.5*[ ((delV)/(delt))_(p) + ((delV)/(delt))_(c)]`

`rho^2 = rho^1 + ((delrho)/(delt))_(avg) * dt`

similarly, we need to solve for momentum as well as energy equations.

Main code

clear all
close all
clc
%Defining inputs for the function
n = 31;
x = linspace(0,3,n);
g = 1.4;
z = linspace(1,1400,1400);

%Calling the functions
[rho_c,t_c,p_c,~,v_c,M_c, rho_c_t, t_c_t, p_c_t, ~, v_c_t, M_c_t] = conservative_function(n,x,g);
[rho_nc,t_nc,p_nc,~,v_nc,M_nc,rho_nc_t, t_nc_t, p_nc_t, ~, v_nc_t, M_nc_t] = non_conservative_function (n,x,g);


%plotting steady state distribution of primitives
figure(1)
     subplot(3,1,1)
     plot(x,rho_c,'.-','color','r','LineWidth',1)
     hold on
     plot(x,rho_nc,'o','color','b','LineWidth',1)
     title({'Number of time steps = 1400';'Comparison of steady state distribution of primitives between Conservative and Non - conservative form'})
     ylabel('Density')
     xlabel('Length along the nozzle (x)')
     legend('Conservative','Non - Conservative')
     
     subplot(3,1,2)
     plot(x,t_c,'.-','color','r','LineWidth',1)
     hold on
     plot(x,t_nc,'o','color','b','LineWidth',1)
     ylabel('Temperature')
     xlabel('Length along the nozzle (x)')
     legend('Conservative','Non - Conservative')
     
     subplot(3,1,3)
     plot(x,p_c,'.-','color','r','LineWidth',1)
     hold on
     plot(x,p_nc,'o','color','b','LineWidth',1)
     ylabel('Pressure')
     xlabel('Length along the nozzle (x)')
     legend('Conservative','Non - Conservative')
     
figure(2)
     subplot(2,1,1)
     plot(x,M_c,'.-','color','r','LineWidth',1)
     hold on
     plot(x,M_nc,'o','color','b','LineWidth',1)
     ylabel('Mach number')
     xlabel('Length along the nozzle (x)')
     legend('Conservative','Non - Conservative')
     title({'Number of time steps = 1400';'Comparison of steady state distribution of primitives between Conservative and Non - conservative form'})
     
     subplot(2,1,2)
     plot(x,v_c,'.-','color','r','LineWidth',1)
     hold on
     plot(x,v_nc,'o','color','b','LineWidth',1)
     ylabel('Velocity')
     xlabel('Length along the nozzle (x)')
     legend('Conservative','Non - Conservative')
     
figure(3)
    subplot(3,1,1)
    plot(z,(rho_c_t),'.','color','r')
    hold on
    plot(z,(rho_nc_t),'.','color','b')
    xlabel('Number of time steps')
    ylabel('Density')
    legend('Conservative','Non - Conservative')
    title('Comprison of Density, Temperature and Pressure for Conservative and Non - Conservative Governing equations')

    subplot(3,1,2)
    plot(z,(t_c_t),'.','color','r')
    hold on
    plot(z,(t_nc_t),'.','color','b')
    xlabel('Number of time steps')
    ylabel('Temperature')
    legend('Conservative','Non - Conservative')

    subplot(3,1,3)
    plot(z,(p_c_t),'.','color','r')
    hold on
    plot(z,(p_nc_t),'.','color','b')
    xlabel('Number of time steps')
    ylabel('Pressure')
    legend('Conservative','Non - Conservative')

figure(4)
    
    subplot(2,1,1)
    plot(z,(M_c_t),'.','color','r')
    hold on
    plot(z,(M_nc_t),'.','color','b')
    xlabel('Number of time steps')
    ylabel('Mach number')
    legend('Conservative','Non - Conservative')
    title('Comprison of Mach Number, and Velocity for Conservative and Non - Conservative Governing equations')
    
    subplot(2,1,2)
    plot(z,(v_c_t),'.','color','r')
    hold on
    plot(z,(v_nc_t),'.','color','b')
    xlabel('Number of time steps')
    ylabel('Velocity')
    legend('Conservative','Non - Conservative')
    
%Grid independence test
n = [31 61 81];
g = 1.4;

for q = 1:length(n)
    if n(q) == 31
    x = linspace(0,3,n(q));
    throat1 = ((n(q)-1)/2)+1;
    tic
    [rho_c1,t_c1,p_c1,m_c1,~,~, ~, ~, ~, ~,~, ~] = conservative_function(n(q),x,g);
    t1 = toc
    tic
    [rho_nc1,t_nc1,p_nc1,m_nc1,~,~,~, ~, ~, ~,~, ~] = non_conservative_function (n(q),x,g);
    t2 = toc
    
    figure(7)
    plot(x,m_c1,'color','r','LineWidth',1)
    hold on
    plot(x,m_nc1,'color','b','LineWidth',1)
    xlabel('Length along the nozzle (x)')
    ylabel('Mass flow rate')
    legend('Conservative','Non - Conservative')
    title('Mass flow rate comparison for Conservative and Non - conservative')
    
    elseif n(q) == 61
    x = linspace(0,3,n(q));
    throat2 = ((n(q)-1)/2)+1;
    tic
    [rho_c2,t_c2,p_c2,m_c2,~,~, ~, ~,~,~,~,~] = conservative_function(n(q),x,g);
    t3 = toc
    tic
    [rho_nc2,t_nc2,p_nc2,m_nc2,~,~,~, ~, ~,~, ~, ~] = non_conservative_function (n(q),x,g);
    t4=toc    
elseif n(q) == 81
    x = linspace(0,3,n(q));
    throat3 = ((n(q)-1)/2)+1;
    tic
    [rho_c3,t_c3,p_c3,m_c3,~,~, ~, ~, ~, ~, ~, ~] = conservative_function(n(q),x,g);
    t5 = toc
    tic
    [rho_nc3,t_nc3,p_nc3,m_nc3,~,~,~, ~, ~,~, ~,~] = non_conservative_function (n(q),x,g);
    t6=toc
    end
end

Temperature_ratio_1c  = t_c1(throat1)/t_c1(1)
Temperature_ratio_2c  = t_c2(throat2)/t_c2(1)
Temperature_ratio_3c  = t_c3(throat3)/t_c3(1)
Temperature_ratio_1nc = t_nc1(throat1)/t_nc1(1)
Temperature_ratio_2nc = t_nc2(throat2)/t_nc2(1)
Temperature_ratio_3nc = t_nc3(throat3)/t_nc3(1)

Pressure_ratio_1c     = p_c1(throat1)/p_c1(1)
Pressure_ratio_2c     = p_c2(throat2)/p_c2(1)
Pressure_ratio_3c     = p_c3(throat3)/p_c3(1)
Pressure_ratio_1nc    = p_nc1(throat1)/p_nc1(1)
Pressure_ratio_2nc    = p_nc2(throat2)/p_nc2(1)
Pressure_ratio_3nc    = p_nc3(throat3)/p_nc3(1)

Mass_ratio_1c         = m_c1(throat1)/m_c1(1)
Mass_ratio_2c         = m_c2(throat2)/m_c2(1)
Mass_ratio_3c         = m_c3(throat3)/m_c3(1)
Mass_ratio_1nc        = m_nc1(throat1)/m_nc1(1)
Mass_ratio_2nc        = m_nc2(throat2)/m_nc2(1)
Mass_ratio_3nc        = m_nc3(throat3)/m_nc3(1)

Density_ratio_1c      = rho_c1(throat1)/rho_c1(1)
Density_ratio_2c      = rho_c2(throat2)/rho_c2(1)
Density_ratio_3c      = rho_c3(throat3)/rho_c3(1)
Density_ratio_1nc     = rho_nc1(throat1)/rho_nc1(1)
Density_ratio_2nc     = rho_nc2(throat2)/rho_nc2(1)
Density_ratio_3nc     = rho_nc3(throat3)/rho_nc3(1)

Conservative Function

function [rho_c,t_c,p_c,m_c,v_c,M_c, rho_c_t, t_c_t, p_c_t, m_c_t, v_c_t, M_c_t] = conservative_function(n,x,g)

dx = x(2) - x(1);
throat = ((n-1)/2)+1;
%calculate initial profiles for rho and t  
for i = 1:n
    if x(i) >= 0 && x(i)<= 0.5
        rho_c(i) = 1;
        t_c(i)   = 1;
        
    elseif x(i) >= 0.5 && x(i) <= 1.5 %elseif used to define conditions
        rho_c(i) = 1 - 0.366 * (x(i)- 0.5);
        t_c(i)   = 1 - 0.167 * (x(i)- 0.5);
        
    elseif 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

%area
a = 1+2.2*(x-1.5).^2; % transpose of a matrix is same as area in anderson textbook

v_c = 0.59./(rho_c.*a);

%define number of time steps
nt = 1400; %number of time steps
c = 0.5; %courant number
dt = min((c .* dx./(sqrt(t_c)+v_c)));  


%Assigning terms for easiness of the code
term1 = rho_c.*a ;
term2 = rho_c.*a.*v_c ;
term3 = rho_c .*a.* ((t_c./(g-1)) + ((g/2).*(v_c.^2)));

%outer time loop
  for k = 1:nt
      
    term1_old = term1; %keeping a copy
    term2_old = term2;
    term3_old = term3;
      
    F1 = term2;
    F2 = ((term2.^2)./term1)+(((g-1)/g)*(term3-((g/2)*((term2.^2)./term1))));
    F3 = ((g*term2.*term3)./term1)-(g*(g-1)*(term2.^3))./(2*term1.^2);
      
%predictor method using forward differencing
    for j = 2:n-1

    da_dx(j) = (a(j+1) - a(j))/dx;
    J2 (j) = (1/g) * (rho_c(j)*t_c(j)) *da_dx(j);
    
    dF1_dx(j) = (F1(j+1)-F1(j))/dx;
    dF2_dx(j) = (F2(j+1)-F2(j))/dx;
    dF3_dx(j) = (F3(j+1)-F3(j))/dx;

    
    %for continuity equation
    d_term1_dt_p(j) = -dF1_dx(j);
    
    d_term2_dt_p(j) = -dF2_dx(j) + J2(j);
    
    d_term3_dt_p(j) = -dF3_dx(j);
    
    %solution update for predictor step
    term1(j) = term1(j) + d_term1_dt_p(j)*dt;
    term2(j) = term2(j) + d_term2_dt_p(j)*dt;
    term3(j) = term3(j) + d_term3_dt_p(j)*dt;
    end

    rho_c = term1./a;
    v_c   = term2./term1;
    t_c   = (g-1)*((term3./term1)-((g*v_c.^2)/2));
    
    %Flux terms for corrector loop
      
    F1 = term2;
    F2 = ((term2.^2)./term1)+(((g-1)/g)*(term3-(g*term2.^2)./(2*term1)));
    F3 = ((g*term2.*term3)./term1)-(g*(g-1)*(term2.^3))./(2*term1.^2);
    
 %Corrector method using backward differencing  
    for j = 2:n-1

    da_dx(j) = (a(j) - a(j-1))/dx;
    J2 (j) = (1/g) * rho_c(j)*t_c(j) *da_dx(j);
    
    dF1_dx(j) = (F1(j) - F1(j-1))/dx;
    dF2_dx(j) = (F2(j) - F2(j-1))/dx;
    dF3_dx(j) = (F3(j) - F3(j-1))/dx;
    
    %for continuity equation
    d_term1_dt_c(j) = -dF1_dx(j);
    
    d_term2_dt_c(j) = -dF2_dx(j) + J2(j);
    
    d_term3_dt_c(j) = -dF3_dx(j);
   
    end
   
    %calculating average
    d_term1_dt = 0.5 * (d_term1_dt_p + d_term1_dt_c);
    d_term2_dt = 0.5 * (d_term2_dt_p + d_term2_dt_c);
    d_term3_dt = 0.5 * (d_term3_dt_p + d_term3_dt_c);
    
    %solution update for corrector step
    for j = 2:n-1
    term1(j) = term1_old(j) + d_term1_dt(j)*dt;
    term2(j) = term2_old(j) + d_term2_dt(j)*dt;
    term3(j) = term3_old(j) + d_term3_dt(j)*dt;
    end
    %Apply Boundary conditions
    %inlet
    term1(1) = rho_c(1)*a(1);
    term2(1) = 2*term2(2) - term2(3);
    term3(1) = term1(1)*((t_c(1)/(g-1))+((g/2)*v_c(1)^2));
    
    %outlet
    term1(n) = 2*term1(n-1) - term1(n-2);
    term2(n) = 2*term2(n-1) - term2(n-2);
    term3(n) = 2*term3(n-1) - term3(n-2);
    
    rho_c = term1./a;
    v_c   = term2./term1;
    t_c   = (g-1)*((term3./term1)-((g*v_c.^2)/2));
    
    m_c = rho_c.*v_c.*a;
    M_c = v_c./sqrt(t_c);
    p_c = rho_c.*t_c;
    
    %Throat values 
    rho_c_t(k) = rho_c(throat);
    t_c_t(k)   = t_c(throat); 
    p_c_t(k)   = p_c(throat);
    m_c_t(k)   = m_c(throat);
    v_c_t(k)   = v_c(throat);
    M_c_t(k)   = M_c(throat);
    %time marching plot for mass flowrate along the length of nozzle
    figure(5)
   if n == 31
    if     k == 50
        plot(x,m_c,'LineWidth',1)
        ylabel('Mass flow rate')
        xlabel('Length along the nozzle (x)')
        title({'Variation of Mass flow rate distribution inside the nozzle at different time steps during the time-marching process','Conservative'})
        hold on
    elseif k == 100 
        plot(x,m_c,'LineWidth',1)
        hold on
    elseif k == 200 
        plot(x,m_c,'LineWidth',1)
        hold on
    elseif k == 400
        plot(x,m_c,'LineWidth',1)
        hold on
    elseif k == 800
        plot(x,m_c,'LineWidth',1)
        hold on
    elseif k == 1000
        plot(x,m_c,'LineWidth',1)
        hold on
    elseif k == 1200
        plot(x,m_c,'LineWidth',1)
        legend('k = 50','k = 100','k = 200','k = 400','k = 800','k = 1000','k = 1200')
    end
    continue
    end
  end
end
     
  

Non - conservative Function

function [rho_nc,t_nc,p_nc,m_nc,v_nc,M_nc,rho_nc_t, t_nc_t, p_nc_t, m_nc_t, v_nc_t, M_nc_t] = non_conservative_function (n,x,g)
 %inputs
  dx     = x(2) - x(1);
  throat = ((n-1)/2)+1;
  %calculate initial profiles
  
  rho_nc = 1 - 0.3146*x; %density
  t_nc   = 1-0.2314*x; %temperature
  v_nc   = (0.1 + 1.09*x).*t_nc.^0.5; %velocity
  
  %area
  a = 1+2.2*(x-1.5).^2; % transpose of a matrix is same as area in anderson textbook
  %plot(a)
  
  %define number of time steps
  nt = 1400; %number of time steps
  c  = 0.5; %courant number
  dt = min((c .* dx./(sqrt(t_nc)+v_nc)));  
  
  %outer time loop
  for k = 1:nt
      
      rho_old = rho_nc; %keeping a copy
      v_old   = v_nc;
      t_old   = t_nc;
      
      %predictor method  Forward differencing
      
    for j = 2:n-1
    dvdx       = (v_nc(j+1) - v_nc(j))/dx;
    drhodx     = (rho_nc(j+1) - rho_nc(j))/dx;
    dtdx       = (t_nc(j+1)-t_nc(j))/dx;
    d_log_a_dx = (log(a(j+1)) - log(a(j)))/dx  ;    
        
    %for continuity equation
    d_rho_dt_p(j) = -rho_nc(j) * dvdx -rho_nc(j)*v_nc(j)* d_log_a_dx - v_nc(j) * drhodx;
  
    %for Momentum equation
    dvdt_p(j) = -v_nc(j)*dvdx - (1/g) * (dtdx + (t_nc(j)/rho_nc(j))* drhodx);
  
    %for energy equation
    dtdt_p(j) = -v_nc(j) * dtdx - (g -1)*t_nc(j)*(dvdx + v_nc(j)* d_log_a_dx);
    
    end
    
    for j = 2:n-1
    %solution update
    rho_nc(j) = rho_nc(j) + d_rho_dt_p(j)*dt;
    v_nc(j)   = v_nc(j) + dvdt_p(j)*dt;
    t_nc(j)   = t_nc(j) + dtdt_p(j)*dt;
    end  
    
    %corrector method  Backward differencing
      
    for j = 2:n-1
  
    %for continuity equation
    dvdx       = (v_nc(j) - v_nc(j-1))/dx;
    drhodx     = (rho_nc(j) - rho_nc(j-1))/dx;
    dtdx       = (t_nc(j)-t_nc(j-1))/dx;
    d_log_a_dx = (log(a(j)) - log(a(j-1)))/dx  ;
    
    d_rho_dt_c(j) = -rho_nc(j) * dvdx -rho_nc(j)*v_nc(j)* d_log_a_dx - v_nc(j) * drhodx;

    %for Momentum equation
    dvdt_c(j) = -v_nc(j)*dvdx - (1/g) * (dtdx + (t_nc(j)/rho_nc(j))* drhodx);
    
    %for energy equation
    dtdt_c(j) = -v_nc(j) * dtdx - (g -1)*t_nc(j)*(dvdx + v_nc(j)* d_log_a_dx);
    
    end   % corrector step is updated
    
    %compute average time derivative
    d_rho_dt = 0.5*(d_rho_dt_p + d_rho_dt_c);
    dvdt     = 0.5*(dvdt_p + dvdt_c);
    dtdt     = 0.5* (dtdt_p + dtdt_c);
    
        for i = 2:n-1
        %final solution update
        v_nc(i)   = v_old(i) + dvdt(i)*dt;
        t_nc(i)   = t_old(i) + dtdt(i)*dt;
        rho_nc(i) = rho_old(i) + d_rho_dt(i)*dt;
        end
        
        %Apply boundary conditions
        %inlet
        v_nc(1) = 2*v_nc(2) -v_nc(3);
        
        %outlet
        v_nc(n)   = 2*v_nc(n-1) - v_nc(n-2);
        rho_nc(n) = 2*rho_nc(n-1) - rho_nc(n-2);
        t_nc(n)   = 2*t_nc(n-1) - t_nc(n-2);
        
        m_nc = rho_nc.*v_nc.*a;
        M_nc = v_nc./sqrt(t_nc);
        p_nc = rho_nc.*t_nc;
        
        %Throat values 
        rho_nc_t(k) = rho_nc(throat);
        t_nc_t(k)   = t_nc(throat); 
        p_nc_t(k)   = p_nc(throat);
        m_nc_t(k)   = m_nc(throat);
        v_nc_t(k)   = v_nc(throat);
        M_nc_t(k)   = M_nc(throat);
    
    figure(6)
    if n == 31
    if     k == 50
        plot(x,m_nc,'LineWidth',1)
        ylabel('Mass flow rate')
        xlabel('Length along the nozzle (x)')
        title({'Variation of Mass flow rate distribution inside the nozzle at different time steps during the time-marching process','Non - conservative'})
        hold on
    elseif k == 100 
        plot(x,m_nc,'LineWidth',1)
        hold on
    elseif k == 200 
        plot(x,m_nc,'LineWidth',1)
        hold on
    elseif k == 400
        plot(x,m_nc,'LineWidth',1)
        hold on
    elseif k == 800
        plot(x,m_nc,'LineWidth',1)
        hold on
    elseif k == 1000
        plot(x,m_nc,'LineWidth',1)
        hold on
    elseif k == 1200
        plot(x,m_nc,'LineWidth',1)
        legend('k = 50','k = 100','k = 200','k = 400','k = 800','k = 1000','k = 1200')
    end 
    continue
    end
  end
end

Explanation and terms used in code :

L = length = 3

g = gamma = 1.4

n = 31, 61, 81 ( 61 and 81 are used for comparison of grid independence test)

dx = 0.1000 ( for n = 31 )

dx = 0.0500 ( for n = 61 )

dx = 0.0375 ( for n = 81 )

c = courant number = 0.5

nt = 1400

`throat = (n-1)/2 + 1`

dt = minimum of `(c*dx)/(sqrtt+v)`

inputs to the function are n,x,g

a = area = `1+2.2*(x-1.5)^2`

For conservative,

`v= 0.59/(rho.a)`

t = 1.....for `0 <= x <= 0.5`

t = `1 - 0.167*(x-0.5)`......for `0.5 <= x <= 1.5`

t = `0.833 - 0.3507 *(x - 1.5)`.....for `1.5 <= x <=3.5`

`rho`= 1.....for `0 <= x <= 0.5`

`rho` = `1 - 0.366*(x - 0.5)`.....for `0.5 <= x <= 1.5`

`rho` = `0.634 - 0.3879*(x-1.5)`.....for `1.5 <= x <=3.5`

For non- conservative

`rho` = `1-0.3146*x`

`t = 1-0.2314*x`

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

Explanation of the MacCormak's method during coding of matlab

`(delterm_1)/(delt') = -(delF1)/(delx')`

`(delterm_2)/(delt') = -(delF2)/(delx') + J_2`

`(delterm_3)/(delt') = -(delF3)/(delx')`

For predictor part ( Forward differencing)

`(delF1)/(delx') = (F1_(i+1) - F1_(i))/dx`

`(delF2)/(delx') = (F2_(i+1) - F2_(i))/dx`

`(delF3)/(delx') = (F3_(i+1) - F3_(i))/dx`

For corrector part ( Backward differencing)

`(delF1)/(delx') = (F1_(i) - F1_(i-1))/dx`

`(delF2)/(delx') = (F2_(i) - F2_(i-1))/dx`

`(delF3)/(delx') = (F3_(i) - F3_(i-1))/dx`

Here, i the number of nodes

RESULTS

1) Steady-state distribution of primitive variables inside the nozzle

Plot showing the steady-state distribution of Density, Temperature, and Pressure along the length of the nozzle.

Plot showing the steady-state distribution of Mach number, and velocity along the length of the nozzle.

Interpretation:

From the above two plots, it can be observed that the accuracy is almost similar for both the methods ( Conservative and non-conservative ) at 1400 time steps. A steady-state solution is obtained at 1400 time steps and the solutions do not change further even on increasing the number of time steps.

2) Time-wise variation of the primitive variables

Plot showing the convergence of Density, Temperature, and pressure values along the time marching process from 1 to 1400 time steps.

Plot showing the convergence of Mach number, and velocity values along the time marching process from 1 to 1400 time steps.

Interpretation:

It can be observed from the above plot that the rate of convergence of primitives and solution is faster in conservative solution as compared to the non - conservative solution.

3) Variation of Mass flow rate distribution inside the nozzle at different time steps during the time-marching process

Variation of mass flow rate at different time steps during the time marching process for the conservative solution

Variation of mass flow rate at different time steps during the time marching process for the conservative solution.

Interpretation:

It can be observed from the above two plots that the rate of convergence is faster in conservative solution as compared to non - conservative solution and hence lower number of time steps are needed for the conservative governing equation. Also, the curves obtained in the conservative solution are much smoother as compared to the non - conservative solutions.

4) Comparison of Normalized mass flow rate distributions of both forms

Plot showing the steady-state mass flow rate in conservative and non-conservative solution

Interpretation:

It can be observed from the above plot that the conservative solution is much more accurate for the mass flow rate through the nozzle as compared to the non - conservative solution for mass flow rate through the nozzle.

Grid independence test :

	Temperature Ratio	Pressure Ratio	Density Ratio	Mass flow rate Ratio	Time (Seconds)
Conservative (31)	0.8410	0.5487	0.6524	1.0001	2.2559
Non - Conservative (31)	0.8365	0.5342	0.6387	0.9932	1.2188
Conservative (61)	0.8356	0.5339	0.6390	1.0002	0.5660
Non - Conservative (61)	0.8354	0.5327	0.6377	0.9997	0.4291
Conservative (81)	0.8350	0.5322	0.6374	1.0016	0.4384
Non - Conservative (81)	0.8355	0.5329	0.6378	1.0101	0.4045
Exact Value	0.833	0.528	0.634	1	-

It can be observed that as the number of grids increases the time required to run the function decreases as shown in the last column of the table. From the calculation part conservative form is more accurate than the non-conservative form which can be clearly seen from the mass flow rate comparison plot at normalized state. But, conservative calculations are more complex as compared to the non-conservative calculations, and hence much more system resources or time is needed to solve the conservative solution.

Hence, proper selection should be done before selecting a governing equation on the basis of solution accuracy, the time required, and system resources available.