Simulation of quasi 1D Super-sonic nozzle flow using MaCcormack Method

Objective:

The objective of this project is to write code solve the 1D supersonic nozzle flow equations using the Macormack Method.

• To implementing both the conservative and non-conservative forms of the governing equations
• To perform grid dependence test

MacCormack method:

MacCormack's technique is used, which provides 2nd order accuracy without the involvement of the 2nd derivative, thus reducing the complexity of the mathematical calculation. We will be setting up a finite difference solution, using the time marching approach.

Simulation of 1D supersonic nozzle flow

Non-conservative form of 1D supersonic nozzle flow (User-defined function)

function[K,rho,t,v,m_n,M_n,p_n,rho_exact_n,t_exact_n,p_exact_n,M_exact_n,rho_th,v_th,t_th,p_th,M_th,rho_res_n,t_res_n,v_res_n]=Non_conservative_f(x,n,gamma,nt,a,dx,dt)
    tic;

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

    %Calculate initial condition
    rho = 1-0.3146*x;
    t=1-0.2314*x;
    v=(0.1+1.09*x).*t.^2;

    %Area
    a=1+2.2*(x-1.5).^2;

    %No.of time step
    nt=1400;
    C=0.5;
    dt=min(C*(dx./(sqrt(t)+v)));
    gamma =1.4;
    
    %Throat value
    th=find(a==1);
    
    %Outer time loop
    for K=1:nt
        rho_old =rho;
        v_old = v;
        t_old = t;
       
        
        %Predictor step
        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 eqation
            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 step
        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 eqation
            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

        %calulate the avg time step
        drhodt = 0.5*(drhodt_p+drhodt_c);
        dvdt = 0.5*(dvdt_p+dvdt_c);
        dtdt = 0.5*(dtdt_p+dtdt_c);

        for i=2:n-1
        %Solution update
        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

        %Boundary condtions
        
        %inlet BC
        v(1)=2*v(2)-v(3);

        %Outlet boundary condition
        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
        m_n=rho.*a.*v;

        %pressure ratio
        p_n=rho.*t;

        %Mach number
        M_n=v./sqrt(t);

        %% @ throat 

        rho_th(K) = rho(th);
        v_th(K) = v(th);
        t_th(K) = t(th);
        p_th(K) =p_n(th);
        M_th(K) = M_n(th);


        %Exact solution
        rho_exact_n=0.634*ones(1,1400);
        t_exact_n=0.833*ones(1,1400);
        p_exact_n=0.528*ones(1,1400);
        M_exact_n=1*ones(1,1400);

        %% Residual convergence
        rho_res_n(K) = drhodt(th);
        t_res_n(K) = dtdt(th);
        v_res_n(K) = dvdt(th);

        %% plot for non dimensional mass flow rate 
        
        figure(1)
        hold on
        if K==1
            plot(x,m_n,'linewidth',1.4)
         hold on
        elseif K==50
                plot(x,m_n,'g','linewidth',1.4)
        hold on
        elseif K==100
                plot(x,m_n,'r','linewidth',1.4)
        hold on
        elseif K==200
                plot(x,m_n,'c','linewidth',1.4)
        hold on
        elseif K==700
                plot(x,m_n,'y','linewidth',1.4)
        hold on
        elseif K==1400
                plot(x,m_n,'--','linewidth',1.4);
            grid on;    
            xlabel('Non - dimensional nozzle dimension','Fontsize',12);
            ylabel('Non- dimensional mass flow rate','Fontsize',12);
            legend('0dt','50dt','100dt','200dt','700dt','1400dt')
            title('Mass flow rate at different time step','Fontweight','bold','Fontsize',12)    
        end
        end
        toc;    
    end
     

Conservative form of 1D supersonic nozzle flow (User-defined function)

function [v_c,rho_c,t_c,m,M_c,pr,v_t,rho_t,t_t,m_t,mach_t,pr_t,rho_res_c,t_res_c,v_res_c] = conservative_from(n,x,nt,dx)
tic;
% input variables
L=3;
n=61;
x = linspace(0,3,n);
dx = x(2)-x(1);
c=0.5;

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

%Initial conditions 
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)
        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


% Assuming steady state flow for the initial condition
  v_c = 0.59./(rho_c.*a);
  gamma = 1.4;
  throat=find(a==1); %
 
 
% Solution vector parameters
  pr = rho_c.*t_c;
  u1 = rho_c.*a;
  u2 = rho_c.*a.*v_c;
  u3 = u1.*((t_c./(gamma-1))+((gamma/2).*(v_c.^2)));

  %calculation of timestep
  dt = min (c.*(dx./((t_c.^0.5)+v_c))); 

  %No.of time cycle 
   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;
        
        
        dadx(i) = (a(i+1)-a(i))/dx;
        j2(i) = (1/gamma)*(rho_c(i)*t_c(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_c = u1./a;
    v_c = u2./u1;
    t_c = (gamma-1)*((u3./u1)-((gamma*v_c.^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;
        
        dadx1(i) = (a(i)-a(i-1))/dx;
        j3(i) = (1/gamma)*(rho_c(i)*t_c(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_c(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_c(1) = u2(1)/u1(1);
    
    % updating u3
    u3(1) = u1(1)*((t_c(1)/(gamma-1))+((gamma/2)*v_c(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_c = u1./a; 
    v_c = u2./u1;
    t_c = (gamma-1)*((u3./u1)-((gamma*v_c.^2)/2));
    
   
    % calculating solutions 
    m = rho_c.*v_c.*a;
    M_c = v_c./((t_c).^0.5);
    p_c = rho_c.*t_c;
    
    % Throat values 
    rho_t(k) = rho_c(throat);
    t_t(k) = t_c(throat); 
    pr_t(k) = p_c(throat);
    mach_t(k) = M_c(throat);
    v_t(k) = v_c(throat);
    m_t(k) = m(throat);
    
    %Residual convergence
    
    rho_res_c(k)=du1dt(throat);
    t_res_c(k)=du2dt(throat);
    v_res_c(k)=du3dt(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('Variation of mass through nozzle at different time steps')
     legend({'50 time steps','100 time steps','150 time steps','200 time steps','550 time steps','700 time steps','1400 time steps'});
   end     
end
toc;
end

Main program of 1D super-sonic nozzle flow

clear all;close all;clc;

%Input varibales 
n = 61; %No.of grid points
x = linspace(0,3,n);
dx = x(2)-x(1);

gamma = 1.4;
nt=1400;%No.of time step


tic;
[K,rho,t,v,m_n,M_n,p_n,rho_exact_n,t_exact_n,p_exact_n,M_exact_n,rho_th,v_th,t_th,p_th,M_th,rho_res_n,t_res_n,v_res_n]=Non_conservative_f(x,n,gamma,nt,dx);
Non_conservative_execution_time=toc;

%hold on

tic;
[v_c,rho_c,t_c,m,M_c,p_c,v_t,rho_t,t_t,m_t,mach_t,pr_t,rho_res_c,t_res_c,v_res_c] = Conservative_from(x,n,nt,dx);
Conservative_execution_time=toc;

%% Conservative flow -Residual plot
figure(8)
plot(rho_res_c)
hold on
plot(t_res_c)
hold on
plot(v_res_c)
xlabel('time step')
ylabel('Density,Temp,Velocity')
title('Residual convergence - Conservative flow')
legend('Density','Temperature','Velocity')
grid on
%% Non-conservation flow - Residual plot
figure(7)
plot(rho_res_n)
hold on
plot(t_res_n)
hold on
plot(v_res_n)
xlabel('time step')
ylabel('Density,Temp,Velocity')
title('Residual convergence - Non conservative flow')
legend('Density','Temperature','Velocity')
grid on

%% Conservative form Time wise variation of density,pressure,temperature & Mach numer
figure(6)
subplot(4,1,1)
plot(rho_t)
hold on
plot(rho_exact_n)
xlabel('Time step')
ylabel('rho/{rho_0}')
title('Time wise variation of density,temperature,pressure and Mach No. at the nozzle throat')
legend('Time wise variation','Exact solution')

subplot(4,1,2)
plot(t_t)
hold on
plot(t_exact_n)
xlabel('Time step')
ylabel('T/{T_0}')
legend('Time wise variation','Exact solution')

subplot(4,1,3)
plot(pr_t)
hold on
plot(p_exact_n)
xlabel('Time step')
ylabel('p/{p_0}')
legend('Time wise variation','Exact solution')

subplot(4,1,4)
plot(mach_t)
hold on
plot(M_exact_n)
xlabel('Time step')
ylabel('Mach no')
legend('Time wise variation','Exact solution')
%% Non-Conservative form Time wise variation of density,pressure,temperature & Mach numer
figure(5)
subplot(4,1,1)
plot(rho_th)
hold on 
plot(rho_exact_n)
xlabel('Time step')
ylabel('rho/{rho_0}')
legend('Time wise variation','Exact solution')
title('Time wise variation of density,temperature,pressure and Mach No. at the nozzle throat')


subplot(4,1,2)
plot(t_th)
hold on 
plot(t_exact_n)
xlabel('Time step')
ylabel('T/{T_0}')
legend('Time wise variation','Exact solution')

subplot(4,1,3)
plot(p_th)
hold on 
plot(p_exact_n)
xlabel('Time step')
ylabel('p/{p_0}')
legend('Time wise variation','Exact solution')

subplot(4,1,4)
plot(t_th)
hold on 
plot(t_exact_n)
xlabel('Time step')
ylabel('Mach no')
legend('Time wise variation','Exact solution')


%% Non-Conservative vs Conservative form mass flow rate
figure(4)
plot(x,m_n,'r','linewidth',1.4)
hold on
plot(x,m,'k','linewidth',1.4)
xlabel('x/L')
ylabel('Non-dimensional mass flow rate')
legend('non-conservative flow','Conservative flow')
grid on;

%% Qualitative aspects of 1D nozzle flow
figure(3)
subplot(4,1,3)
plot(x,rho,'r')
hold on
plot(x,rho_c,'*')
xlabel('Length of the nozzle')
ylabel('Density')
legend('Non-conservative','conservative')
grid on;

subplot(4,1,4)
plot(x,t,'r')
hold on
plot(x,t_c,'*')
xlabel('Length of the nozzle')
ylabel('Temperature')
legend('Non-conservative','conservative')
grid on;

subplot(4,1,1)
plot(x,M_n,'r')
hold on
plot(x,M_c,'*')
xlabel('Length of the nozzle')
ylabel('Mach number')
legend('Non-conservative','conservative')
title('Qualitative aspects of quasi-1D nozzle flow')
grid on;

subplot(4,1,2)
plot(x,p_n,'r')
hold on
plot(x,p_c,'*')
xlabel('Length of the nozzle')
ylabel('Mach number')
legend('Non-conservative','conservative')
grid on;

%%

Result & Discussion

A.VARIATION FLOW THROUGHOUT THE NOZZLE

Both forms produce almost similar plots, once after the attainment of steady-state. 1400 time steps are carried out by the time marching approach, at which the steady-state is almost obtained and the values of the primitive variables are in good approximation with the exact analytical solution.

B.TIMEWISE VARIATION OF FLOW PARAMETER ALONG THE NON-DIMENSIONAL LENGTH OF NOZZLE AT THROAT

1. On a deep observation, it can be found that the solution of the non-conservative form lies closer to the exact values of rho', T', pr' and mach' when compared to the conservative form solution.

2. Tracking of the values of flow parameters at a particular region is a good way of investigating the attainment of the steady-state by the system.

3. The non-conservative form generates more accurate results for the primitive variables.

4. But even after 1400 time steps, residue (Difference between the numerical solution and the exact solution) exists. For a non-conservative form, the value of residue will be almost constant after 1400 time steps and equal to 10^-6. For a conservative form, the value of residue will be almost constant after 1400 time steps and equal to 10^-3.

5. Non-conservative form leads to the 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.

6. From the plots below, it can be seen that the steady-state is almost obtained with a good level of approximation, after 500-time steps for practical purposes. But this evidence is not enough to justify our observation.

7. For further investigation, we are going to observe the mass flow rate (non-dimensional) behavior along the non-dimensional distance 'x' at different steps.

1.NON-CONSERVATIVE FORM

2.CONSERVATIVE FORM

C.NON-DIMENSIONAL MASS FLOW RATE BEHAVIOUR ALONG NON-DIMENSIONAL LENGTH OF NOZZLE

1. From the following figure, it can be observed that the mass flow rate curve obtained for 700 & 1400 time steps are almost overlying in the same region of space.

2. So this plot substantiates our previous assumption from the observation of the "Timewise variation of flow parameters at the nozzle throat", i.e after 500 time steps almost steady state is obtained (for practical purposes) for the non-conservative form of the flow equations.

3. Theoretically, it takes infinite time steps to attain a steady-state.

4. Then what will be the case for conservative form? It will be discussed with relevant plots obtained from that case.

1.NON-CONSERVATIVE FORM

Number of grid points, n=31

Number of grid points, n=61

2.CONSERVATIVE FORM

1. When compared to the previous plot, the fluctuations are less severe in the conservative form. But comparing these two plots makes no sense because we assumed different sets of initial
conditions. In fact, we made an intelligent guess while choosing the initial conditions(density, temperature, velocity) with constant mass flow, for the conservative form because of its sensitivity of the stability behavior. This is the reason for its mild fluctuating behaviour.

2. Even from this plot, it can be observed that the mass flow rate curve obtained for 550,700 & 1400 time steps are almost overlying in the same region of space.

3. But, the reasonable settlement of the fluctuations at the 550th time step may be due to the more realistic initial conditions with constant mass flow, assumed by the intelligent guess. And also for the conservative form, the value of residue at 550th-time step will be somewhere between 10^-1(Initially) and 10^-3(after 1400 timesteps), which doesn't lead to a good approximation with the exact solution. i.e in other words, @ 550th-time step, milder fluctuations still exist, which can't be approximated with the exact solution at a good level of accuracy, due to the presence of large values of residue.

4. After 700 time steps, it almost (but not quite) becomes equal to a constant value. Moreover, it is quite close to the exact analytical solution of 0.579, and therefore the steady-state is almost obtained (for practical purposes). Theoretically, it takes infinite time steps to attain a steady-state.

Number of grid points,n=31

Number of grid points,n=61

D.COMPARISION OF MASS FLOW RATE OF NON-CONSERVATIVE FORM AND CONSERVATIVE FORM

1. The conservative form generates a better mass flow distribution with less fluctuations, when compared with the non-conservative form.

2. The reason is the mass flow itself is one of the dependent variables in the equations of the conservative form, whereas in the non-conservative form, the dependent variables are the primitive variables and the mass flow is obtained after some calculation.

3. The conservative form generates more accurate results for the flux variables.

Number of grid points,n=31

Number of grid points, n=61

RESIDUAL

WHAT ARE THE FACTORS WHICH INFLUENCE THE RESIDUE? HOW TO REDUCE THE VALUE OF RESIDUE?
There are three reasons for the small numerical inaccuracies(residue). They are
1. Small inflow boundary condition error:
We assumed that the density, temperature, and pressure at the inlet are equal to the reservoir properties. This is only true if the mach number at the inlet is equal to zero. But mach number is a function of the area ratio A/A*. Even at the inlet finite value exists for the area ratio(5.95), and hence there must be a finite mach number at inlet (x=0), both numerically and analytically. It is a built-in error and it is not viewed as serious.

2. Courant number effects:
According to the Courant number stability criteria, the CFL number must be less than or equal to 1. But we assumed C=0.5 which is smaller-than-necessary at which the solution is stable and the results are closer to the exact analytical solution (calculations are not presented, as it is out of the scope of this project). So, we can conclude that we made a safe guess. CFL criteria provide a good estimate to calculate the time step value for even nonlinear hyperbolic partial differential equations (Quasi 1D nozzle flow), but it does not hold exactly because it is derived for linear hyperbolic equations.

3. Grid dependence:
a. Even If the number of grid points is doubled, and there is only a marginal difference in the solution, then our original calculations with the initial no.of grid points is said to be GRID-INDEPENDENT

b. This factor is purely based on the circumstances, depending upon the expected order of accuracy and the allowable range of tolerance.

c. In our case we assumed n=31, so lets check whether it is grid-independent or not. Initially, to start the discussion it would be better if we plot the graphs with respect to the modified value of 'n', i.e n=61. when modifying the value of 'n', the number of time steps should also be modified, so that the elapsed non-dimensional time is the same for both cases. By doing so, the grid independence test can be carried out by comparing the curves obtained for both cases, at a particular non-dimensional time.

d. The total non-dimensional time elapsed (n=31), is equal to dt*no.of time steps. Initially 'dt' is calculated to all the grid points from i=1 to n using the CFL criteria, and then the smallest of all the values is chosen as 'dt'. By following this approach the same value of dt will be used throughout the time marching. So that there won't be artificial 'time warp' and it provides a time-accurate solution of the actual transient flow field.

e. The value of the non-dimensional time step (n=31) is 0.0201 and 0.02036 for non-conservative and conservative forms respectively. Lets choose the higher value for consideration. Therefore the total non-dimensional time elapsed, 0.02036*1400 = 28.50.

f. The value of the non-dimensional time step (n=61) is found similarly and it is equal to 0.010181. So, the simulation should run for higher no.of time steps to reach the non-dimensional time of magnitude 28.50. The reason why it takes higher no.of time steps is due to the fact that the value of the time step reduces as the no.of grid points are increased.

g. No.of time steps required for n=61 is (28.50 / 0.010181) = 2800 time steps approximately. It seems to be almost double in our case, i.e in other words, for n=31 the curves produced at 100th time step, 500th time step & 1400th time step will be similar, to the curves generated for n=61 at 200th time step, 1000th time step & 2800th time step respectively, if our initial assumption n=31 is GRID-INDEPENDENT.

h. Grid independence test can be also carried out by comparing the values obtained by conservative and non-conservative forms at the throat, for n = 31 & 61. It has to be noted that the non-conservative form results are more accurate when compared with the conservative form results. Even though by 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.

Grid dependence test

1.NON-CONSERVATIVE FORM

Number of grid points, n=31

Number of grid points, n=61

CONSERVATIVE FORM

Number of grid points, n=31

Number of grid points, n=61

Which method is faster?

when n=31, c=0.5

simulation time(s)

When n=61, c=0.5

In the case of a conservative form, primitive variables need to be decoded from the flux variables, but this is not necessary in the case of non-conservative form. And therefore solving conservative form requires more computational effort.