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

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

Objective: In this problem, we are going to write a code to solve the Quasi 1D supersonic nozzle flow simulation for both conservative and Non-conservative forms of equations using the MacCormack method.

• You will perform grid dependence test

• Figure out the minimum number of cycles for which the simulation should be run in order for convergence
• Compare the normalized mass flow rate between the conservative forms and non-conservative forms - make suitable conclusions
• Which method is faster?

1D Nozzle:

Definition: Turbomachines like steam turbines, water turbines, and gas turbines produce power by utilizing the kinetic energy of the jets produced by passing high-pressure steam, water, and gas through the devices called nozzles. Correspond­ing to the fluids used, the nozzles are called steam nozzles, water nozzles, and gas nozzles.

The nozzle is a device used to give the direction to the exhaust gases coming out from the combustion chamber. It converts low velocity, high pressure, high temperature of the gas to high velocity, low pressure and low temperature of the gas in the respective combustion chamber.

We consider the steady, isentropic flow through a convergent-divergent nozzle as sketched. The flow at the inlet to the nozzle comes from a reservoir where the pressure and temperature are denoted by p0, and To, respectively. The cross-sectional area of the reservoir is large ( theoretically, A --> oo ), and hence the velocity is very small ( V --> O). Thus, Po and T0 are stagnation values. The flow expands isentropically to supersonic speeds at the nozzle exit, where the exit pressure, temperature, velocity, and Mach number are denoted by Pe, Te, Ve and Me, respectively.

The flow is locally subsonic in the convergent section of the nozzle is sonic at the throat
(minimum area) and supersonic at the divergent section. The sonic flow (M = 1) at the throat means that the local velocity at this location is equal to the speed of sound.

Using an asterisk to denote sonic flow values, we have at the throat V = V*, a = a*. Similarly, the sonic flow values of pressure and temperature are denoted by p* and T*, respectively.

Conservation vs Non-Conservation:

The conservative equation is applicable to the stationary system either fixed control volume fixed in space or infinitesimally small element fixed in space. Whereas non-conservative equations are used for the moving system. Finite control volume fixed mass moving with the flow or infinitesimally small fluid element of fixed mass (fluid particle) moving with the flow. That is why total derivative is used in non-conservative equations.

If the solution is non-smooth, i.e., there are shocks, then the conservative equations should be used and discretized in a conservative way. If the solution is smooth, even at supersonic speeds but more commonly at subsonic speeds, the two equations are equivalent. Hence it does not matter which form is chosen. The important property is provable numerical stability. If proven stable, two discretizations, conservative and non-conservative will yield very similar answers given that the orders of accuracy are the same.

MacCormack Method:

In computational fluid dynamics, the MacCormack method is a widely used discretization scheme for the numerical solution of hyperbolic partial differential equations. It is an explicit scheme that is second-order accurate and consists of two steps namely predictor and corrector. Predictor uses a forward difference scheme and the corrector uses a backward difference scheme.

It is a much simpler method compared to the other numerical methods which are dealing especially with the non-linear equations. There are some other methods powerful than this, but this method is fully student-friendly.

In MacCormack method, the value of the variable is represented in a way of

In numerical analysis, predictor-corrector methods belong to a class of algorithms designed to integrate ordinary differential equations – to find an unknown function that satisfies a given differential equation. All such algorithms proceed in two steps:

1. The initial, \"prediction\" step, starts from a function fitted to the function-values and derivative-values at a preceding set of points to extrapolate (\"anticipate\") this function\'s value at a subsequent, new point.

1. The next, the \"corrector\" step refines the initial approximation by using the predicted value of the function and another method to interpolate that unknown function\'s value at the same subsequent point.

The governing equations in the form of conservative and Non-Conservative equations are as follows:

Non-Conservative Equations:

• ContinuityEquation

`(partialrho)/(partialt) = -rho*(partialv)/(partialx)-rho*v*(partialloga)/(partialx)-v*(partialrho)/(partialx)`

• Momentum Equation

`(partialv)/(partialt) = -v*(partialv)/(partialx)-1/gamma*((partialt)/(partialx)+t/rho*(partialrho)/(partialx))`

• Energy Equation

`(partialt)/(partialt) = v*(partialt)/(partialx)-(gamma-1)*t((partialv)/(partialx)+v*(partialloga)/(partialx))`

Note: above equation is a time marching problem. So, please dont confuse as here I represent both temperature and time as \'t\'.

To calculate the number of time steps;

`dt = C*(dx)/(v+a)`

where C = Courant number as 0.5. (for the stability)

v = velocity

The initial conditions given for the problem is:

`rho = 1-0.3146*x`

`t = 1-0.2314*x`

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

`gamma = 1.4`

pressure `p = rho*t`

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

Here we need to give the boundary conditions at subsonic and supersoinic sections

• Subsonic section`U_3=`

`v_1 = 2*v_2+v_3`

• Supersonic section

`v_n = 2*v_(n-1)+v_(n-2)`

`t_n = 2*t_(n-1)+t_(n-2)`

`rho_n = 2*rho_(n-1)+rho_(n-2)`

• For the grid generation, we devided the whole length of the nozzle into 31 and 61 parts along the x-axis to study the grid independance. From previous challenges, we know that less the grid size (dx) the error is also goinbg to be less. Lets see that phenomena here. The discretization technique applied is MacCormach method to solve the equations.

Conservative Equations:

• ContinuityEquation

`(partial/(partialt)*rho*a) = -(partial/(partialx)*rho*a*v)`

• Momentum Equation

`(partial/(partialt)*rho*a*v)+(partial/(partialx)[rho*a*v^2+(1/gamma)*rho*a]= 1/gamma*rho*(partiala)/(partialx)`

• Energy Equation

`partial/(partialt)[rho*(e/(gamma-1)+(gamma/2)*v^2)*a]+partial/(partialx)[rho*(e/(gamma-1)+(gamma/2)*v^2)v*a+rho*a*v]=0`

Above equations are non-dimensional conservative form equations. These equations are quite unrealistic to see and hard to write the program with these. So, we need to write these in a generic form that can be easy and useful to write code.

The equations are devided into three parts as

• Solution vector

`U_1=rho*a`

`U_2 = rho*a*v`

`U_3=[rho*(e/(gamma-1)+(gamma/2)*v^2)*a]`

• Flux vector

`F_1=rho*a*v`

`F_2=rho*a*v^2+1/gamma*rho*a`

`F_3=[rho*(e/(gamma-1)+(gamma/2)*v^2)v*a+rho*a*v]`

• Source term

`J_2=1/gamma*rho*(partiala)/(partialx)`

From these generic equatins we can reconstruct the continuity, momentum and energy equations as;

• Continuity equation

`(partialU_1)/(partialt)=-(partialF_1)/(partialx)`

• Momentum equation

`(partialU_2)/(partialt)=-(partialF_2)/(partialx)+J_2`

• Energy equation

`(partialU_3)/(partialt)=-(partialF_3)/(partialx)`

Main Program:

• The code written below is a numerical form content that published in the John. D Anderson book.
• Density, velocity, temperature and area formulas were taken from the textbook itself.

For n = 31:

clear all
close all
clc

% inputs
n = 61;
x = linspace(0,3,n);
dx = x(2)-x(1);
gamma = 1.4;

%Time step 
nt = 1000;
c1 = 0.5;
c2 = 0.5;
tol = 1e-6;
mass_f_tol1 = 1e-3;
mass_f_tol2 = 1e-2;

% Non-conservative function

[rho_nc,v_nc,t_nc,time_taken1,net_mf_rate_nc,error_v_nc,error_rho_nc,error_t_nc,m_dot_nc,p_nc] = Nonconservative(n,x,dx,gamma,nt,c1,tol,mass_f_tol1);

% Conservative function

 [rho_c,v_c,t_c,time_taken2,net_mf_rate_c,error_v_c,error_rho_c,error_t_c,m_dot_c,p_c] = conservative(n,x,dx,gamma,nt,c2,tol,mass_f_tol2);

% Plots

% velocity
figure(1)
plot(x,v_nc,\'b\',x,v_c,\'r\')
xlabel(\'Non-dimensional Distance\')
ylabel(\'Non-dimensional Velocity\')
legend(\'Non-Conservative\',\'Conservative\')
grid on

% density
figure(2)
plot(x,rho_nc,\'b\',x,rho_c,\'r\')
xlabel(\'Non-dimensional Distance\')
ylabel(\'Non-dimensional Density\')
legend(\'Non-Conservative\',\'Conservative\')
grid on

% temperature
figure(3)
plot(x,t_nc,\'b\',x,t_c,\'r\')
xlabel(\'Non-dimensional Distance\')
ylabel(\'Non-dimensional Temperature\')
legend(\'Non-Conservative\',\'Conservative\')
grid on

% pressure
figure(4)
plot(x,p_nc,\'b\',x,p_c,\'r\')
xlabel(\'Non-dimensional Distance\')
ylabel(\'Non-dimensional Pressure\')
legend(\'Non-Conservative\',\'Conservative\')
grid on

% Net mass flow rate 
figure(5)
subplot(2,1,1)
plot(net_mf_rate_nc)
xlabel(\'Time Step\')
ylabel(\'Net Mass Flow Rate\')
title(\'Non-Conservative Form\')
grid on
subplot(2,1,2)
plot(net_mf_rate_c)
xlabel(\'Time Step\')
ylabel(\'Net Mass Flow Rate\')
title(\'Conservative Form\')
grid on

% Error in velocity 
figure(6)
subplot(2,1,1)
plot(error_v_nc)
xlabel(\'Time Step\')
ylabel(\'Error in Velocity\')
title(\'Non-Conservative Form\')
grid on
subplot(2,1,2)
plot(error_v_c)
xlabel(\'Time Step\')
ylabel(\'Error in Velocity\')
title(\'Conservative Form\')
grid on

% Error in density
figure(7)
subplot(2,1,1)
plot(error_rho_nc)
axis([0 inf -0.01 0.02])
xlabel(\'Time Step\')
ylabel(\'Error in Density\')
title(\'Non-Conservative Form\')
grid on
subplot(2,1,2)
plot(error_rho_c)
axis([0 inf -0.005 0.015])
xlabel(\'Time Step\')
ylabel(\'Error in Density\')
title(\'Conservative Form\')
grid on

% Error in temperature
figure(8)
subplot(2,1,1)
plot(error_t_nc)
axis([0 inf -0.005 0.015])
xlabel(\'Time Step\')
ylabel(\'Error in Temperature\')
title(\'Non-Conservative Form\')
grid on
subplot(2,1,2)
plot(error_t_c)
axis([0 inf -0.01 0.02])
xlabel(\'Time Step\')
ylabel(\'Error in Temperature\')
title(\'Conservative Form\')
grid on

figure(11)
hold on
plot(x,m_dot_c,\'b\')
plot(x,m_dot_nc,\'r\')
line([0 3], [0.579 0.579],\'color\',\'magenta\')
legend(\'Conservative\',\'Non-Conservative\',\'Exact Solution\')
title(\'Non-Dimensional Mass Flow Distribution\')
xlabel(\'Length of Nozzle\')
ylabel(\'Mass Flow Rate\')
grid on

Function Program for Non-Conservative form:

    function [rho,v,t,k,net_mf_rate,error_v,error_rho,error_t,m_dot_nc,p_nc] = Nonconservative(n,x,dx,gamma,nt,c,tol,mass_f_tol1)

% Intial conditions
% Density
rho = (1 - 0.3146*x);  
% Tempaerature
t = (1 - 0.2314*x);
% Velocity
v = (0.1 + 1.09*x).*(t.^0.5);
% Area
a = (1 + 2.2*(x-1.5).^2)

% outer time loop

for k = 1:nt
    
     dt = min(abs(c*dx./(realsqrt(t)+v)));

% updating initial conditions
     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;
        dloga_dx = (log(a(j+1))-log(a(j)))/dx;
        dtdx = (t(j+1)-t(j))/dx;
    
        % continuty equation
    
        drhodt_p(j) = -rho(j)*dvdx- rho(j)*v(j)*dloga_dx- 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)*dloga_dx);
        
         % solution update
         
          v(j) = v(j)+dvdt_p(j)*dt ;
        rho(j) = rho(j)+ drhodt_p(j)*dt;
        t(j) = t(j) + dtdt_p(j)*dt;

    
   end
   
      % corrctor 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)-v(j-1))/dx;
        drhodx = (rho(j)-rho(j-1))/dx;
        dloga_dx = (log(a(j))-log(a(j-1)))/dx;
        dtdx = (t(j)-t(j-1))/dx;
        
        % continuty equation
    
        drhodt_c(j) = -rho(j)*dvdx- rho(j)*v(j)*dloga_dx- 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)*dloga_dx);
        

    
   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 j = 2:n-1
        v(j) = v_old(j)+dvdt(j)*dt;
        rho(j) = rho_old(j)+drhodt(j)*dt;
        t(j) = t_old(j)+dtdt(j)*dt;
   end
   
   % applying 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);
   
   %Mass flow rate
       
       inflow = rho(1)*v(1)*a(1);
       outflow = rho(n)*v(n)*a(n);
       net_mf_rate(k) = (outflow-inflow);
       m_dot_nc = rho_old.*a.*v_old;
       p_nc = rho_old.*t_old;

       
   %Error tolerance
     error_v(k) = max(v - v_old);
     error_rho(k) = max(rho - rho_old);
     error_t(k) = max(t-t_old);
      if(error_v(k) < tol && error_rho(k) < tol && error_t(k) < tol && net_mf_rate(k))

      end
      
      if n == 61
          figure(10)
          hold on
       if k == 100
           plot(x,m_dot_nc,\'r\')
       elseif k == 200
           plot(x,m_dot_nc,\'b\')
       elseif k == 500
           plot(x,m_dot_nc,\'g\')
       elseif k == 700
           plot(x,m_dot_nc,\'black\')
           line([0 3], [0.579 0.579],\'color\',\'magenta\')
           title (\'Comparison of Mass Flow Distribution for different time steps (Non-Conservative)\')
           legend(\'100 delta(t)\',\'200 delta(t)\',\'500 delta(t)\',\'700 delta(t)\',\'Exact Solution\')
           xlabel(\'Length of the Nozzle\')
           ylabel(\'Mass Flow Rate\')
           grid on
       end
      end
     
end

Function program for Conservative form:

• Here some guidelines and instructions were followed while choosing the density, velocity, and temperature are from the textbook.
• The area is same. Area doesn\'t change even though we are changing the equation form.

function [rho,v,t,k,net_mf_rate,error_v,error_rho,error_t,m_dot_c,p_c] = conservative(n,x,dx,gamma,nt,c,tol,mass_f_tol2)

% Intial conditions
% Density
% Temerature

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

% Area
a = 1 + 2.2*(x-1.5).^2;
% Velocity
v = 0.59./(rho.*a);
% pressure
p_c = rho.*t

% Generic form of continuty, momentum and energy equations

U1 = rho.*a;
U2 = rho.*a.*v;
U3 = rho.*(t/(gamma-1)+(gamma/2)*v.^2).*a;

% outer time loop

for k = 1:nt
    
    % updating generic equations
    
    U1_old = U1;
    U2_old = U2;
    U3_old = U3;
    
    rho_old = U1_old./a;
    v_old =  U2_old./(rho.*a);
    t_old = (gamma-1)*((U3_old./U1_old) - (gamma/2)*v.^2);
  
    
    % flux vector
    
     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);
     
     % source term (S)
     
     for j = 2:n-1
         S(j) = (1/gamma)*rho(j)*t(j)*(a(j+1) - a(j))/dx;
         
     end
     
     for i = 1:n
         deltat(i) = (c*(dx/((t(i)^0.5+v(i)))));
         
     end
     
     dt = min (deltat);
     
     
     % Predictor method
     
     for j = 2:n-1
         
         % Continuty equation
          dU1dt_p(j) = -(F1(j+1)-F1(j))/dx;

         % Momentum equation
         dU2dt_p(j) = -(F2(j+1)-F2(j))/dx+S(j);
         
         % Energy equation
         dU3dt_p(j) = -(F3(j+1)-F3(j))/dx;
         
         %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
     
     
         % calculating variables
     
         rho = U1./a;
         v = U2./(rho.*a);
         t = (gamma-1)*((U3./U1) - (gamma/2)*v.^2);
         
         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);
         
         
         for j = 2:n-1
             
             S(j) = (1/gamma)*rho(j)*t(j)*(a(j) - a(j-1))/dx;
         end
         
     % Corrector method
         
          for j = 2:n-1
             
          % Continuty equation
          dU1dt_c(j) = -(F1(j)-F1(j-1))/dx;

          % Momentum equation
          dU2dt_c(j) = -(F2(j)-F2(j-1))/dx+S(j);
         
          % Energy equation
          dU3dt_c(j) = -(F3(j)-F3(j-1))/dx;
         
          end
         
         %computing average time 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
         
         % applying boundary conditions
         
         %inlet 
         U2(1) = 2*U2(2) - U2(3);
         
         % 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);

         % recalculating variables
         
         rho = U1./a;
         v = U2./(rho.*a);
         t = (gamma-1)*((U3./U1) - (gamma/2)*v.^2) ;
         
         %Mass flow rate
       
       inflow = rho(1)*v(1)*a(1);
       outflow = rho(n)*v(n)*a(n);
       net_mf_rate(k) = (outflow-inflow); 
       m_dot_c = rho_old.*a.*v_old
       p_c = rho_old.*t_old;
 
       %Error tolerance
     error_v(k) = max(v - v_old);
     error_rho(k) = max(rho - rho_old);
     error_t(k) = max(t-t_old);
     if (error_v(k)  <  tol && error_rho(k)  <  tol && error_t(k)  <  tol && net_mf_rate(k))
       
     end
     
     if n == 61
         
      figure(9)
      hold on
       if k == 100
           plot(x,m_dot_c,\'r\')
       elseif k == 200
           plot(x,m_dot_c,\'b\')
       elseif k == 500
           plot(x,m_dot_c,\'g\')
       elseif k == 700
           plot(x,m_dot_c,\'black\')
           line([0 3], [0.579 0.579],\'color\',\'magenta\')
           title (\'Comparison of Mass Flow Distribution for different time steps (Conservative)\')
           legend(\'100 delta(t)\',\'200 delta(t)\',\'500 delta(t)\',\'700 delta(t)\',\'Exact Solution\')
%            ylim([0.5 0.7])
           xlabel(\'Length of the Nozzle\')
           ylabel(\'Mass Flow Rate\')
           grid on
       end
     end
           
     
end

Results: (for n = 31)

Comparison of velocity, temperature, and density using both conservative & Non-Conservative forms:

• The flow of the fluid is assumed to be one dimensional.

• As fluid enters the smaller cross-section, it has to speed up due to the conservation of mass. To maintain a constant amount of fluid moving through the restricted portion of the nozzle, the fluid must move faster. The energy to make this fluid speed up has to come from somewhere. Some energy is in the random motion of molecules, which we observe as pressure.

• A duct of varying cross-section area in which a steadily flowing fluid can be made to accelerate by a pressure drop along the duct. So, when a fluid flows through a nozzle, its velocity increases continuously and pressure decreases continuously. The above graphs can be a replica for the statement written.
• from the mass flow rate equation, mass flow rate (m_dot) is decreasing, the area is constant and velocity is decreasing. so, theoretically and numerically the density of the medium is also decreased.
• The drop in enthalpy leads to a corresponding drop in temperature.

• Graphs are shown above represent that in conservative or in Non-Conservative forms, the error ratio of mass flow rate, velocity, pressure, density, and temperature are going to be Zero while increasing in the number of time steps.
• By comparison, we can say that the Conservative form converges speedily and not much crusts and troughs in the graph (which makes solution unstable).
• Maximum 800 timesteps are enough to converge the solution. After 800 timesteps, solutions are not going to change much.
• But for the understanding of the viewer, 1000 timesteps were calculated.

• These comparison graphs of Non-Conservative and Conservative forms by changing the number of delta(t)\'s show how the equations behaving at different conditions.

Grid Dependancy:

a) n=31

b) n=61

• Figure a) and b) represents the comparison of conservative and non-conservative results with the exact value.
• From the observance, we can say that in both 31 and 61 elements the conservative form is much closer to the exact solution and shows a decent line of solution.
• Where in non-conservative form, it shows that the solution is nowhere near to the exact solution. The reason is it has a lot of variations consisting of irregular oscillations at inlet and outlet.
• Even though we change the value for a number of elements, the solution doesn\'t look much vary. So applying n = 31 can give similar results and we can notice that computational time is also less.

Conclusion:

• The mass flow distribution ratio seems to be at a good rate in the conservative form when compared.
• When the solution starts, in non-dimensional cases, Non-Conservative equation shows some unwanted up\'s and down\'s in a graph, wherein the conservative form we have seen very less of them.
• This is because of the presence of flux vectors and the source term in the conservative form. They make a solution much smoother.
• By this we cant sat that non-cons form is useless. Both forms are the same and equal.
• The non-conservative form gives faster convergence at the primary or initial conditions. It doesn\'t happen in conservative form because of the same flux and source terms.
• To compute numerically, Non-Conservative form is quite easier, easily understandable and takes less time to compute compare to conservative form.

Note:

1) Individual code for the conservative form is attached in zip format.

2) Other useful plots of n = 31, and n = 61 is attached in zip format.