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

The main purpose of this project is to compare the different numerical schemes such as \"Conservative Form\" and \"Non-Conservative Form\" for 1D Subsonic-Supersonic Nozzle flow simulations.

1. Assumption:

• Steady State & Isentropic flow
• Quasi 1D flow
• Inlet is a reserviour @ `P_{0}, T_{0}` Area of reserviour is `Ararroo`

In actual the flow will be 2D i,e in XY Plane. For convient of writing program we consider the variables remain constant along the Y-axis (at particular length of nozzle the flow variables remain constant along the cross sectional area of nozzle).

Nozzle area will vary parabolically along the length of nozzle as shown in fig.

2. The Setup of the problem contain:

1. The Governing Flow Equations:

1. Conservation of mass : `(del(rhoA))/(delt) + (del(rhoAV))/(delx) = 0`
2. Conservation of momentum: `rho(delV)/(delt) + rhoV(delV)/(delx) = -A(delp)/(delx)`
3. Conservation of Energy: `rho(dele)/(delt) + rhoV(dele)/(delx) = -p(delV)/(delx) - pV(del(lnA))/(delx)`

The above equcations are called Navier-stokes equcations. The equcation can be return in Non-dimensional form w.r.t to values @ reserviour. Aerodynamicist prefare both dimensionalised and Non-dimensional form however many will go for the Non-dimensional form because it easy to study the result where as to get the physical intution of magnitude of flow variables we can go for the dimensional form. One can find the numerical form of the governing equcation in 7th chapter of book by John D Anderson \"Computational Fluid Dynamics: Basic and Applications\"

2. The Finite-Difference Equations: We will follow the MacCormack time marching technique for finite difference expressions i,e Predicator Corrector Method

3. Calculation of time step: We will consider CFL criteria to find the time step value at grid points

`Deltat = C. (Deltax)/(a + V)`

Where: `Deltat rarr`time step @ a grid point

`C rarr Courant Number`

`a rarr`Acoustic Velocity

`V rarr`Local Velocity

There are two methods to calculate the time step

1. Local time step method: At each grid point we will calculate the time and use it calculate the flow variables.
   
2. Minimum time step: We select the lowest possible value among all time step at each nodes. We will follow this method because it provides common time step for all nodes and variables will be stable and approaches steady state fast.
   

3. Inital Boundry Conditions: There are 2 methods to provide the intial boundry conditions

1. Methods of characteristic & stream line direction: In this method of characteristic tells us that at a boundary where one chracteristic propagates into the domain then the values of one dependent flow field variable must be specified at that boundary and if one characteristic line propagates out of the domain then the values of another dependent flow field variable must be allowed to float at the boundary. At outflow boundary all variables are allowed to float because of supersonic velocity it points all charecteristic out of domain.
2. Boundary Conditions for Numerical Method Form: Boundary conditions are set based on linear extrapolation.Intial boundary conditions in following program @ t=0 are:

`rho = 1 - 0.3146x`

`T = 1 - 0.2314x`

`V = (0.1+1.09x)T^(1/2)`

The following Main program call for function \"Conservative Form of Solver\" & \"Non-Conservative Form of Solver\"

Main Program:

clear all
close all
clc

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

% Area profile

A = 1 + 2.2*((x-1.5).^2);
n_throat = find(A-1<0.001);

% Intial profile

C = 0.5;
error_noncon = 1e9;
error_con = 1e9;
tol = 1e-5;


% Non Conservative solver

tic

[rho,v,T,p,M,massflow,rho_throat,T_throat,M_throat,p_throat,count,drho_dt_throat,dv_dt_throat,dt_dt_throat] = Quasi_nonconservative_nozzleflow(n,x,dx,gamma,A,n_throat,C,error_noncon,tol);

time_noncon = toc;

fprintf(\'Number of Itteration for convergence in Non-conservative form: %f\\n\',count);
fprintf(\'Computation time for Non-conservative form: %f s\\n\',time_noncon);

% Transient variation of the throat properties
Transient_Throat=figure(2);
subplot(4,1,1);
plot(rho_throat,\'m\',\'linewidth\',1.0);
title(\'Transient variation of flow variables @ throat (Non-conservative part)\');
legend(\'Density ratio\');
ylabel(\'\\rho/\\rho_{0}\');
axis([0 count 0.4 1]);
grid minor;

subplot(4,1,2);
plot(T_throat,\'r\',\'linewidth\',1.0);
legend(\'Temperature ratio\');
ylabel(\'T/T_{0}\');
axis([0 count 0.6 1]);
grid minor;

subplot(4,1,3);
plot(p_throat,\'g\',\'linewidth\',1.0);
legend(\'Pressure ratio\');
ylabel(\'p/p_{0}\');
axis([0 count 0.1 1]);
grid minor;

subplot(4,1,4);
plot(p_throat,\'b\',\'linewidth\',1.0);
legend(\'Mach Number\');
ylabel(\'M\');
axis([0 count 0.2 1]);
xlabel(\'Number of time steps\');
grid minor;

% Flow variables along the length of nozzle
figure(3)
subplot(4,1,1);
plot(x,rho,\'m\',\'linewidth\',1.0);
title(\'Flow variables along the nozzle length Non-conservative part\');
ylabel(\'\\rho/\\rho_{0}\');
legend(\'Density ratio\');
axis([0 3 0 1]);
grid minor;

subplot(4,1,2);
plot(x,T,\'r\',\'linewidth\',1.0);
ylabel(\'T/T_{0}\');
legend(\'Temperature ratio\');
axis([0 3 0 1]);
grid minor;

subplot(4,1,3);
plot(x,p,\'g\',\'linewidth\',1.0);
ylabel(\'p/p_{0}\');
legend(\'Pressure ratio\');
axis([0 3 0 1]);
grid minor;

subplot(4,1,4);
plot(x,M,\'b\',\'linewidth\',1.0);
ylabel(\'M\');
legend(\'Mach Number\',\'location\',\'northwest\');
axis([0 3 0 4]);
xlabel(\'x/L\');
grid minor;

figure(4);
plot(abs(drho_dt_throat),\'r\');
hold on;
plot(abs(dv_dt_throat),\'b\');
hold on;
plot(abs(dt_dt_throat),\'k\');
title(\'Timewise variation of the absolute values of the time derivatives @ Throat (Non-Conservative Form)\');
legend(\'drho__dt\',\'dv__dt\',\'dT__dt\');
axis([0 count -0.1 0.8]);
ylabel(\'Residual\');
xlabel(\'Number of time step\');
grid minor;

%Conservative Solver

tic 

[rho_con,v_con,T_con,p_con,M_con,massflow_con,rho_throat_con,T_throat_con,M_throat_con,p_throat_con,count_con] = Quasi_Conservative_nozzleflow(n,x,dx,gamma,A,n_throat,C,error_con,tol);

time_con = toc;

fprintf(\'Number of Itteration for convergence in Conservative form: %f\\n\',count_con);
fprintf(\'Computation time for Conservative form: %f s\\n\',time_con);

% Transient variation of the throat properties
Transient_Throat=figure(6);
subplot(4,1,1);
plot(rho_throat_con,\'m\',\'linewidth\',1.0);
title(\'Transient variation of flow variables @ throat (conservative part)\');
legend(\'Density ratio\');
ylabel(\'\\rho/\\rho_{0}\');
axis([0 count 0.5 0.8]);
grid minor;

subplot(4,1,2);
plot(T_throat_con,\'r\',\'linewidth\',1.0);
legend(\'Temperature ratio\');
ylabel(\'T/T_{0}\');
axis([0 count 0.7 1]);
grid minor;

subplot(4,1,3);
plot(p_throat_con,\'g\',\'linewidth\',1.0);
legend(\'Pressure ratio\');
ylabel(\'p/p_{0}\');
axis([0 count 0.4 .8]);
grid minor;

subplot(4,1,4);
plot(p_throat_con,\'b\',\'linewidth\',1.0);
legend(\'Mach Number\');
ylabel(\'M\');
axis([0 count 0.4 0.8]);
xlabel(\'Number of time steps\');
grid minor;

% Flow variables along the length of nozzle
figure(7)
subplot(4,1,1);
plot(x,rho_con,\'m\',\'linewidth\',1.0);
title(\'Flow variables along the nozzle length Conservative part\');
ylabel(\'\\rho/\\rho_{0}\');
legend(\'Density ratio\');
axis([0 3 0 1.2]);
grid minor;

subplot(4,1,2);
plot(x,T_con,\'r\',\'linewidth\',1.0);
ylabel(\'T/T_{0}\');
legend(\'Temperature ratio\');
axis([0 3 0 1.2]);
grid minor;

subplot(4,1,3);
plot(x,p_con,\'g\',\'linewidth\',1.0);
ylabel(\'p/p_{0}\');
legend(\'Pressure ratio\');
axis([0 3 0 1.2]);
grid minor;

subplot(4,1,4);
plot(x,M_con,\'b\',\'linewidth\',1.0);
ylabel(\'M\');
legend(\'Mach Number\',\'location\',\'northwest\');
axis([0 3 0 4.0]);
xlabel(\'x/L\');
grid minor;

figure(8)
plot(x,massflow_con,\'r\',\'linewidth\',1.25);
hold on;
plot(x,massflow,\'b\',\'linewidth\',1.25);
title(\'Comparison of Conservative form and Non-Conservative form Massflow rate\');
legend(\'Conservative form\',\'Non-Conservative form\');
xlabel(\'x/L\');
ylabel(\'\\rhoVA/\\rho_{0}V_{0}A_{0}\');
axis([0 3 0.58 0.61]);
grid minor;

Function for Non-Conservation form of solver:

function [rho,v,T,p,M,massflow,rho_throat,T_throat,M_throat,p_throat,count,drho_dt_throat,dv_dt_throat,dt_dt_throat] = Quasi_nonconservative_nozzleflow(n,x,dx,gamma,A,n_throat,C,error_noncon,tol)

rho = 1-0.3146.*x;
T = 1-0.2314.*x;
v = (0.1 + 1.09*x).*T.^0.5;
count = 0;

while error_noncon>tol
    
    count=count+1;
    
    rho_old = rho;
    v_old = v;
    T_old = T;
    a_old = real(sqrt(T_old));
    massflow_old = rho_old.*v_old.*A;    
    
    % CFL Number
    
    time_step = C*dx./(a_old+v_old);
    dt = min(time_step);
    
    %Predicator method 
    
    for j = 2:n-1
        
        dvdx = (v(j+1)-v(j))/dx;
        dlogadx = (log(A(j+1))-log(A(j)))/dx;
        drhodx = (rho(j+1)-rho(j))/dx;
        dtdx = (T(j+1)-T(j))/dx;
        
        % Continutiy Eqn
       
        drho_dt_p(j) = -rho(j)*dvdx - rho(j)*v(j)*dlogadx - v(j)*drhodx;
        
        % Momentum Eqn
        
        dv_dt_p(j) = -v(j)*dvdx - (1/gamma)*(dtdx + T(j)/rho(j)*drhodx);
        
        % Energy Eqn
        
        dt_dt_p(j) = -v(j)*dtdx - (gamma-1)*T(j)*(dvdx + v(j)*dlogadx);
        
        %Solution update 
         
        rho(j) = rho(j) + drho_dt_p(j)*dt;
        v(j) = v(j) + dv_dt_p(j)*dt;
        T(j) = T(j) + dt_dt_p(j)*dt;
 
    end
    
    %Corrector method 
    
    for j = 2:n-1
        
        dvdx = (v(j)-v(j-1))/dx;
        dlogadx = (log(A(j))-log(A(j-1)))/dx;
        drhodx = (rho(j)-rho(j-1))/dx;
        dtdx = (T(j)-T(j-1))/dx;
        
        % Continutiy Eqn
       
        drho_dt_c(j) = -rho(j)*dvdx - rho(j)*v(j)*dlogadx - v(j)*drhodx;
        
        % Momentum Eqn
        
        dv_dt_c(j) = -v(j)*dvdx - (1/gamma)*(dtdx + (T(j)/rho(j))*drhodx);
        
        % Energy Eqn
        
        dt_dt_c(j) = -v(j)*dtdx - (gamma-1)*T(j)*(dvdx + v(j)*dlogadx);
 
    end
    
    % Avg time derivative 
       
       drho_dt = 0.5*(drho_dt_p + drho_dt_c);
       dv_dt = 0.5*(dv_dt_p + dv_dt_c);
       dt_dt = 0.5*(dt_dt_p + dt_dt_c);
       
   for i = 2:n-1
       
       rho(i) = rho_old(i) + drho_dt(i)*dt;
       v(i) = v_old(i) + dv_dt(i)*dt;
       T(i) = T_old(i) + dt_dt(i)*dt;
       
   end   
   
   % Apply the BC
   % Inlet BC
   rho(1) = 1;
   v(1) = 2*v(2) - v(3);
   T(1) = 1;
   
   % Outlet BC
   rho(n) = 2*rho(n-1) - rho(n-2);
   v(n) = 2*v(n-1) - v(n-2);
   T(n) = 2*T(n-1) - T(n-2);
    
   % Calculation of other variables 
   
   M = v./(real(sqrt(T)));
   massflow = rho.*v.*A;
   p = rho.*T;
   
   % Transient Properties @ Throat 
   
   rho_throat(count) = rho(n_throat);
   T_throat(count) = T(n_throat);
   p_throat(count) = p(n_throat);
   M_throat(count) = M(n_throat);
   drho_dt_throat(count) = drho_dt(n_throat);
   dv_dt_throat(count) = dv_dt(n_throat);
   dt_dt_throat(count) = dt_dt(n_throat);
   
   error_noncon = max(abs(massflow-massflow_old));
   
   TransientMassFlow = figure(1);
   hold on;
   
   if count == 1
       plot(x,massflow,\'--m\',\'linewidth\',1.0);
   else if count == 50
           plot(x,massflow,\'k\',\'linewidth\',1.0);
   else if count == 100
           plot(x,massflow,\'g\',\'linewidth\',1.0);
   else if count == 150
           plot(x,massflow,\'b\',\'linewidth\',1.0);
   else if count == 200
           plot(x,massflow,\'y\',\'linewidth\',1.0);
   else if count == 700
           plot(x,massflow,\'-rp\',\'linewidth\',1.5);
           plot(x(end),massflow(end),\'ks\');
   end
   end
   end
   end
   end
   end
   
end

title(\'Transient Non-dimensional mass flow rates for Non-Conservative Form\');
xlabel(\'x/L\');
ylabel(\'\\rhoVA/\\rho_{0}V_{0}A_{0}\');
legend(\'0\\Deltat\',\'50\\Deltat\',\'100\\Deltat\',\'150\\Deltat\',\'200\\Deltat\',\'700\\Deltat\',\'Steady-flow analytic result\');
axis([0 3 0 2]);
grid minor;

end

Function for Conservative form of solver:

function [rho_con,v_con,T_con,p_con,M_con,massflow_con,rho_throat_con,T_throat_con,M_throat_con,p_throat_con,count_con] = Quasi_Conservative_nozzleflow(n,x,dx,gamma,A,n_throat,C,error_con,tol)

%Intial profile of Area of nozzle

x_1 = find(x==0.5);
x_2 = find(x==1.5);
x_3 = find(x==3.0);

for i=1:x_1
    rho_con(i) = 1.0;
    T_con(i) = 1.0;
end

for i=x_1:x_2
    rho_con(i) = 1.0 - 0.366*(x(i)-0.5);
    T_con(i) = 1.0 - 0.167*(x(i)-0.5);
end

for i=x_2:x_3
    rho_con(i) = 0.634 - 0.3879*(x(i)-1.5);
    T_con(i) = 0.833 - 0.3507*(x(i)-1.5);
end

% Solution vector
U1 = rho_con.*A;
v_con = 0.59./(rho_con.*A);
U2 = rho_con.*v_con.*A;
U3 = rho_con.*A.*((T_con./(gamma-1))+(gamma/2)*(v_con.^2));
p_con = rho_con.*T_con;

count_con = 0;

while error_con>tol
    
    count_con = count_con + 1;
    
    v_old_con = v_con;
    rho_old_con = rho_con;
    T_old_con = T_con;
    a_old_con = real(sqrt(T_con));
    massflow_old_con = rho_con.*v_con.*A;

    U1_old = U1;
    U2_old = U2;
    U3_old = U3;
    
    % 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);
    
    % CFL Number 
    
    timestep_con = C*dx./(a_old_con+v_con);
    dt = min(timestep_con); 
    
    % Predicator Step

for i=2:n-1

    
    J2(i) = (1/gamma)*rho_con(i)*T_con(i)*(A(i+1)-A(i))/(dx);
        
    dU1_dt_P(i) = -(F1(i+1)-F1(i))/dx;
    dU2_dt_P(i) = (-(F2(i+1)-F2(i))/dx)+J2(i);
    dU3_dt_P(i) = -(F3(i+1)-F3(i))/dx;    
    
    U1(i) = U1(i) + dU1_dt_P(i)*dt;
    U2(i) = U2(i) + dU2_dt_P(i)*dt;
    U3(i) = U3(i) + dU3_dt_P(i)*dt;
    
end
    % Predicated values of flow variable
    rho_con = U1./A;
    v_con = U2./U1;
    T_con = (gamma-1)*((U3./U1)-(gamma/2)*(U2./U1).^2);
    p_con = rho_con.*T_con;
    
    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 method

for i=2:n-1
    
    J2_1(i) = ((1/gamma)*rho_con(i)*T_con(i)*(A(i)-A(i-1))/dx);
    dU1_dt_C(i) = -(F1(i)-F1(i-1))/dx;
    dU2_dt_C(i) = (-(F2(i)-F2(i-1))/dx)+J2_1(i);
    dU3_dt_C(i) = -(F3(i)-F3(i-1))/dx;    
    
end
    
    dU1_dt_avg = 0.5*(dU1_dt_P+dU1_dt_C);
    dU2_dt_avg = 0.5*(dU2_dt_P+dU2_dt_C);
    dU3_dt_avg = 0.5*(dU3_dt_P+dU3_dt_C);

for i=2:n-1
    
    U1(i) = U1_old(i) + dU1_dt_avg(i)*dt;
    U2(i) = U2_old(i) + dU2_dt_avg(i)*dt;
    U3(i) = U3_old(i) + dU3_dt_avg(i)*dt;
end
    
    % Inlet BC
    U1(1) = rho_con(1)*A(1);
    U2(1) = 2*U2(2) - U2(3);
    U3(1) = U1(1)*(T_con(1)/(gamma-1) + (gamma/2)*(v_con(1).^2));

    % Outlet BC
    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_con = U1./A;
    v_con = U2./U1;
    T_con = (gamma-1)*((U3./U1)-(gamma/2)*v_con.^2);   
    p_con = rho_con.*T_con;
    
    % Other Variables
    
    M_con = v_con./(real(sqrt(T_con)));
    massflow_con = rho_con.*v_con.*A;
    p_con = rho_con.*T_con;
    
    % Flow variables @ throat section
    
    rho_throat_con(count_con) = rho_con(n_throat);
    M_throat_con(count_con) = M_con(n_throat);
    p_throat_con(count_con) = p_con(n_throat);
    T_throat_con(count_con) = T_con(n_throat);
    
    error_con = max(abs(massflow_con - massflow_old_con));
    
    TransientMassFlow_con=figure(5);
    hold on;
    
    if count_con == 1
        plot(x,massflow_con,\'--r\',\'linewidth\',1.0);
    else if count_con == 100
            plot(x,massflow_con,\'k\',\'linewidth\',1.0);
    else if count_con == 200
            plot(x,massflow_con,\'m\',\'linewidth\',1.0);
    else if count_con == 700
            plot(x,massflow_con,\'.-b\',\'linewidth\',1.5);
            plot(x(end),massflow_con(end),\'rs\');
    end
    end
    end
    end
    
end

title(\'Transient Non-dimensional mass flow rates for Conservative Form\');
xlabel(\'x/L\');
ylabel(\'\\rhoVA/\\rho_{0}V_{0}A_{0}\');
legend(\'0\\Deltat\',\'100\\Deltat\',\'200\\Deltat\',\'700\\Deltat\',\'Steady-flow analytic result\');
axis([0 3 0.5 0.7]);
grid minor;

end

4. The Steady State Solutions:

Figure(1):

Figure(2):

Figure(3):

Figure(4):

• Tempreature (T), Pressure (p), Density (`rho`) and Mach number (M) are ploted at throat vs time step as shown in fig. (1) & (3) for both Non-conservative & Conservative form respectively.
• The largest variation take place at early times, after which the final steady-state value is approached almost asymptotically.
• Here one can observe the \"rubber band effect\", at early times the rubber band is \"stretched\" tightly because of which flow field variables are driven by a stronger potential and changes rapidly where as after few time step the error get reduced because of which rubber band \"relaxed\" causing the flow variables changes much smaller with time.
• The non-dimensional flow variables are `rho`, `T`, `M` and `p`changes along the length of nozzle as shown in fig. (2) & (4) for both Non-conservative & Conservative form respectively.
• As nozzle designed for the subsonic-supersonic flow, non-dimensional flow variables `rho`, `T`, and `p`decrease along the length of nozzle and `M` increases along the length for both Non-conservative & Conservative form.
  
• The Conservation form & Non-conservation form are theortically suitable representation of fundemental physical principles of mass, momentum and energy conserved.
• However there are different advantage of both form such as conservation form can capture the shock were as Non-conservation form leads to poor results.

5. Non-dimensional Massflow variation with timestep:

1. Non-conservative form:

Figure(5):

ii. Conservative form:

Figure(6):

• The distribution of the nondimensional mass flow through the nozzle at different times is shown for both Non-Conservative and Conservative form as shown in fig. (5) & (6) respectively.
• The curve labeled `0Deltat` corresponds to the initial condition.
• After some time steps the massflow starts to converge towards steady state which can be seen as horizontal line in graph.
• After `50Deltat` there are still large peaks are present they will gradually start to reduce to following time steps `100Deltat`, `150Deltat`, `200Deltat`, `700Deltat` and massflow will get converged to steady state after 767 time steps for Non-conservative form.
• Conservative form take 723 time steps to converge to steady state.
• The steady state reached by conservative form is shown by square dot which is 0.586 (Analytical massflow is 0.590).
• Conservative form leads to much accurate massflow in comparison to Non-conservative form.
• However in my part due to style of code writing there is 4.3% error in Numerical and Analytical massflow, this can be corrected by improving the code style.
• Below Fig.(7) shows the comparison of Non-dimensional Non-Conservative & Conservative mass flow.

iii. Comparison of Non-dimensional massflow:

Figure(7):

6. Timewise variations of the absolute values of the time derivatives of non-dimensional density, velocity and tempreature

Figure(8):

• The following graph is obtain for timestep of 5000 (tolerence = 1e-16) where as other plots are obtained for around 700 timesteps.
• The plot is obtained for Residual vs Number of time step.
• As early times, the time derivative are large and they oscillate in value. These oscillations are occurs because of various compressions and expenssions in nozzle.
• At later times, the time derivatives rapidly grow small & steady state is achieved after 3000 time steps.
• Steady flow variables (Primitive & Flux) are achieved after the 500 time steps (i,e almost 3 decimal accuracy).

7. Demonstration of Grid Independence

• The grid independence refers to the flow variables whose values will not change with N.
• In genral we will get some values for the quntity for N number of grid points but this values may get change for 2N grid points hence we should carryout the test for different N values.
• After particular N number the variables doesnot change this condition is called grid independence.
• In our case both conservative & Non-conservative form are grid independent as shown in table 1.

Table 1:

	Condition @ the nozzle throat
	Density ratio	Temperature ratio	Pressure ratio	Mach Number
Exact analytical solutions	0.634	0.833	0.528	1
Case 1: 31 points Non-conservative form	0.639	0.836	0.534	0.999
Case 2: 61 points Non-conservative form	0.638	0.835	0.533	0.999
Case 1: 31 points Conservative form	0.645	0.788	0.508	1.057
Case 2: 61 points Conservative form	0.645	0.793	0.511	1.046

8. Courant Number Effect:

• As we know the CFL criteria C`<=`1 does not hold good for non-linear equcation still we got the stable result till C=1.1 as shown in table 2.
• The program will become unstable and blow up for C=1.2.
• However from the below results we can see that CFL criterion is certainly a good estimate for values of timestep (`Deltat`).
• It is the most reliable estimate for `Deltat`that we can use, even though the governing equations are nonlinear.

Table 2:

	Courant Number Effects (Non-conservative form)
Courant Number (C)	Density ratio	Temperature ratio	Pressure ratio	Mach Number
Exact analytical solution	0.634	0.833	0.528	1
0.5	0.639	0.836	0.528	0.999
0.7	0.639	0.837	0.535	0.999
0.9	0.639	0.837	0.535	0.999
1.1	0.639	0.837	0.535	0.999
1.2	Program went unstable

9. Results

1. In reality the conservation form leads to better massflow distribution but In my case its non-conservative form.
2. Non-conservative form leads to the smaller residuals. The amount by which residual decay is often used as an indedx of quality of the numerical algorithm. In this sense the nonconservative form does a better job.
3. There is no clear superiority of either form in terms of accuracy of the results. The non-conservative form produce the accurate result for primitive variables and conservative form produces accurate results for flux variables.
4. The conservative form required more effort while writing the algorithm compare to Non-conservative form because it have to decode the flux variables to primitive variables.