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

Simulation of a 1D Super-sonic Nozzle Flow Simulation using Macormack Method

Objective:-

To Simulate the isentropic flow through a Quasi 1D subsonic nozzle by deriving both the conservation and Non-conservation forms of the governing equations and solve them by implementing the MacCormack's technique using MATLAB.

Problem Statement:-

Determine the steady-state temprature distribution for the flow -field variables and to investigate the diffrence between the two forms of governing equations by comparing their solutions along with plots stated below.

1. Steady-State distribution of primitive variable inside the nozzle.
2. Time-wise variation of the primitive variable.
3. Variation of Mass rate distribution inside the nozzle at diffrent time steps during the time-marching process.
4. Comparasion of Normalized mass flow rate distributions of both forms.

Theory:-

Any numerical solution of isentropic flow though a quasi 1D subsonic-supersonic Nozzle is reasonable, while analytic solution exists for these kind of problems. However we are trying here to establish the broadened capability of CFD in various streams rather in representing the approximated form of analytic solution. Here, we implement the numerical solution with conservative and Non-conservative forms.

Governing Equations are marched in terms of PDE for the both forms inoredr to use Time Marching solution for the problem. These forms represent in numerical aspect only, if we analytic all the methods are same.

Non Conservative Form:-

Continuity Equation : `(partialrho')/(partialt') = rho'(partialV')/(partialx')-rho'V'(partial(lnA))/(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' (partiallnA')/(partialx')]`

Conservative Form:-

Continuity Equation : `(partialrho'A')/(partialt') + (partialrho'A'V')/(partialx') = 0`

Momentum Equation : `(partial(rho'A'V'))/(partialt') + (partial(rho'A'V^2 + (1/gamma) rho'A'))/(partialx') = (1/gamma)rho'(partialA')/(partialx')`

Energy Equation : `(partial[rho'((T')/(gamma-1) + (gamma/2V^2)A']))/(partialt') + (partial[rho'((T')/(gamma-1) + (gamma/2V^2)V'A']))/(partialx')`

Non Dimensional Quantities:-

Also as part of the problem setup we represent quantites in terms of non-dimensional quantities.

Temprature `T' = (T/T_0)`

Density `rho' = (rho/rho_0)`

Non Dimensional Distance `x' = (x/L)`

Non Dimensional Velocity `V' = (V/a_0)`

Non Dimensional Time `t' = (t/(L/a_0))`

Non Dimensional Area `A' = (A/A^*)`

Nozzle diagram with sections:-

Since the problem setup is 1D, we consider only distribution in x-axis. Here Nozzle can be split into two sections

• Convergent section
• Divergent section

MacCormack method:-

To discretize the numerical solution we use the techniques called Maccarmack method. The Maccarmack method is a variation of the two-step Lax-Wenddroff Scheme but is much simpler in application

• Predictor step
• Corrector step

MAIN PROGRAM:-

In order for the code to be modulur, program is split into three sections. They are

• Main program of Quasi 1D subsonic-supersonic nozzle
• Non Conservative form
• Conservative form

MATLAB CODE:-

Main Program of Solving Quasi 1D subsonic -supersonic nozzle:-

% Main Program of solving Quasi 1D subsonic-supersonic Nozzle
clear all
close all
clc
%Solving Quasi 1D SuperSonic Nozzle Flow Equations
% Input Parameters
%No.of grid points
n = 31;
gamma = 1.4;
x = linspace(0,3,n);
dx = x(2) - x(1);
% No. of Time step
nt = 1400;
%Courant Number
C = 0.5;
% Function for Non Concervative Form
tic;
[Mass_flow_rate_NC, Pressure_NC, Mach_number_NC, rho_NC, v_NC, T_NC, rho_throat_NC, v_throat_NC, T_throat_NC, Mass_flow_rate_throat_NC, Pressure_throat_NC, Mach_number_throat_NC] = Non_conservative_form(x, dx, n, gamma, nt, C)
Time_elapsed_NC = toc;
fprintf('Time Elapsed for Non Conservative form: %0.4f', Time_elapsed_NC);
hold on
% Plotting
%Time variation of properties at throat area in Non Conservative form
figure(2)
%Density
subplot(4,1,1)
plot(linspace(1,nt,nt),rho_throat_NC,'color','m')
ylabel('Density')
legend({'Density at throat'});
axis([0 1400 0 1.33]);
grid minor;
title('Timestep variation at throat area - NON CONSERVATIVE FORM')
%Temprature
subplot(4,1,2)
plot(linspace(1,nt,nt),T_throat_NC,'color','c')
ylabel('Temprature')
legend({'Temprature at throat'});
axis([0 1400 0.6 1]);
grid minor;
%Pressure
subplot(4,1,3)
plot(linspace(1,nt,nt),Pressure_throat_NC,'color','g')
ylabel('Pressure')
legend({'Pressure at throat'});
axis([0 1400 0.4 1.1]);
grid minor
%Mach Number
subplot(4,1,4)
plot(linspace(1,nt,nt),Mach_number_throat_NC,'color','y')
xlabel('Timesteps')
ylabel('Mach Number')
legend({'Mach Number at throat'});
axis([0 1400 0.6 1.4]);
grid minor;

% Stedy state simulation for primitive values for Non Conservative form
figure(3);
% Density
subplot(4,1,1);
plot(x,rho_NC,'color','m')
ylabel('Density')
legend({'Densiy'});
axis([0 3 0 1])
grid minor;
title('Variation of Flow Rate Distribution - NON CONSERVATIVE')
%Temprature
subplot(4,1,2);
plot(x,T_NC,'color','c')
ylabel('Temprature')
legend({'Temprature'});
axis([0 3 0 1])
grid minor;
%Pressure 
subplot(4,1,3);
plot(x,Pressure_NC,'color','g')
ylabel('Pressure')
legend({'Pressure'});
axis([0 3 0 1])
grid minor;
%Mach Number
subplot(4,1,4);
plot(x,Mach_number_NC,'color','y')
xlabel('Nozzle Length')
ylabel('Mach number')
legend({'Mach number'});
axis([0 3 0 4])
grid minor;
% Function for CONSERVATIVE FORM
tic;
[Mass_flow_rate_C, Pressure_C, Mach_number_C, rho_C, v_C, T_C, rho_throat_C, v_throat_C,T_throat_C, Mass_flow_rate_throat_C, Pressure_throat_C, Mach_number_throat_C] = Conservative_form(x, dx, n, gamma, nt, C);
Time_elasped_C = toc;
fprintf('Time Elasped for Conservative form: %0.4f', Time_elasped_C);
%Plotting
%Time step variation of properties at throat area in Conservative form
figure(5)
hold on
%Density
subplot(4,1,1)
plot(rho_throat_C,'color','m')
ylabel('Density')
legend({'Density at throat'});
axis([0 1400 0 1]);
grid minor;
title('Timestep variation at throat area - CONSERVATIVE FORM')
%Temprature
subplot(4,1,2)
plot(T_throat_C,'color','C')
ylabel('Temprature')
legend({'Temprature at throat'});
axis([0 1400 0.6 1]);
grid minor;
%Pressure
subplot(4,1,3)
plot(Pressure_throat_C,'color','g')
ylabel('Pressure')
legend({'Pressure at throat'})
axis([0 1400 0.2 0.8]);
grid minor
%Mach Number
subplot(4,1,4)
plot(Mach_number_throat_C,'color','y')
xlabel('Timesteps')
ylabel('Mach Number')
legend({'Mach Number at throat'});
axis([0 1400 0.6 1.4])
grid minor;

% Steady state simulation for primitive values for Non Conservative form
figure(6);
%Density
subplot(4,1,1);
plot(x,rho_C,'color','m')
ylabel('Density')
legend({'Density'});
axis([0 3 0 1])
grid minor;
title('Variation of Flow Rate Distribution - CONSERVATIVE FORM')
%Temprature
subplot(4,1,2);
plot(x,T_C,'color','c')
ylabel('Temprature')
legend({'Temprature'})
axis([0 3 0 1])
grid minor;
%Pressure
subplot(4,1,3);
plot(x,Pressure_C,'color','g')
ylabel('Pressure')
legend({'Pressure'});
axis([0 3 0 1])
grid minor;
%Mach number
subplot(4,1,4);
plot(x,Mach_number_C,'color','y')
xlabel('Nozzle Length')
ylabel('Mach number')
legend({'Mach number'})
axis([0 3 0 4])
grid minor;

%Comparasion Plot of Normalized Mass Flow Rate of both forms
figure(7)
hold on
plot(x,Mass_flow_rate_NC,'r','linewidth',1.25)
hold on;
plot(x,Mass_flow_rate_C,'b','linewidth',1.25)
hold on;
grid on
legend('Non Conservative','Conservative');
xlabel('Length of the Domain')
ylabel('Mass Flow Rate')
title('Comparasion of Noemalized Mass Flow Rate of both forms')

Explanation:-

• Code is started with clear all, close all, clc to enable a fresh start without the manipulation of previous values.
• Input parameters are assigned for the variables x, dx, n, gamma, nt and C.
• Non conservative and conservative forms are created as sepearte functions with the neccessary inputs which is detailed below this section.
• tic toc command is used for calculating the time elasped between the functions.
• Plotting is done for the following solution.
• Each plot computed as subplots of four primitive Variables such Density, Temprature, Pressure and Mach number at throat as well as Convergent/Divergent section.

NON-CONSERVATIVE FORM

% Function for Non-Conservative Form
function [Mass_flow_rate_NC, Pressure_NC, Mach_number_NC, rho_NC, v_NC, T_NC, rho_throat_NC, v_throat_NC, T_throat_NC, Mass_flow_rate_throat_NC, Pressure_throat_NC, Mach_number_throat_NC] = Non_conservative_form(x, dx, n, gamma, nt, C)
% Initial Boundary conditions
rho_NC = 1 - 0.3146*x;
T_NC = 1 - 0.2341*x;
v_NC = (0.1 + 1.09*x).*T_NC.^0.5;

% Area & Throat profiles
a = 1 + 2.2*(x - 1.5).^2;
throat = find(a==1);

% Outer Time Loop
for k = 1:nt
    % Time-step control
    dt = min(C.*dx./(sqrt(T_NC) + v_NC));
    
    %Time-step control
    dt = min(C.*dx./(sqrt(T_NC) + v_NC));
    % Saving a copy of old values
    rho_old = rho_NC;
    T_old = T_NC;
    v_old = v_NC;
    %Predictor Method
    for j = 2:n-1
        
        % Simplifying the Equation's spatial terms
        dv_dx = (v_NC(j+1) - v_NC(j))/dx;
        drho_dx = (rho_NC(j+1) - rho_NC(j))/dx;
        dT_dx = (T_NC(j+1) - T_NC(j))/dx;
        dlog_a_dx = (log(a(j+1)) - log(a(j)))/dx;
        
        %Continuiity Equation
        drho_dt_P(j) = -rho_NC(j)*dv_dx - rho_NC(j)*v_NC(j)*dlog_a_dx - v_NC(j)*drho_dx;
        
        %Momentum Equation
        dv_dt_P(j) = -v_NC(j)*dv_dx - (1/gamma)*(dT_dx + (T_NC(j)*drho_dx/rho_NC(j)));
        
        %Energy Equation
        dT_dt_P(j) = -v_NC(j)*dT_dx - (gamma - 1)*T_NC(j)*(dv_dx + v_NC(j)*dlog_a_dx);
        
        %Solution update
        rho_NC(j) = rho_NC(j) + drho_dt_P(j)*dt;
        v_NC(j) = v_NC(j) + dv_dt_P(j)*dt;
        T_NC(j) = T_NC(j) + dT_dt_P(j)*dt;
    end
    
    % Corrector Method
    for j = 2:n-1
        % Simplifying the Equation's spatial terms
        dv_dx = (v_NC(j) - v_NC(j-1))/dx;
        drho_dx = (rho_NC(j) - rho_NC(j-1))/dx;
        dT_dx = (T_NC(j) - T_NC(j-1))/dx;
        dlog_a_dx = (log(a(j)) - log(a(j-1)))/dx;
        
        %Continuiity Equation
        drho_dt_C(j) = -rho_NC(j)*dv_dx - rho_NC(j)*v_NC(j)*dlog_a_dx - v_NC(j)*drho_dx;
        
        %Momentum Equation
        dv_dt_C(j) = -v_NC(j)*dv_dx - (1/gamma)*(dT_dx + (T_NC(j)*drho_dx/rho_NC(j)));
        
        %Energy Equation
        dT_dt_C(j) = -v_NC(j)*dT_dx - (gamma - 1)*T_NC(j)*(dv_dx + v_NC(j)*dlog_a_dx);
    end
    
    %Compute Average Time Deivative
    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);
    
    %Final Solution Update
    for i = 2:n-1
        rho_NC(i) = rho_old(i) + drho_dt(i)*dt;
        v_NC(i) = v_old(i) + dv_dt(i)*dt;
        T_NC(i) = T_old(i) + dT_dt(i)*dt;
    end
    
    %Applying Boundary conditions
    %At Inlet
    rho_NC(1) = 1;
    v_NC(1) = 2*v_NC(2) - v_NC(3);
    T_NC(1) = 1;
    
    %At Outlet
    rho_NC(n) = 2*rho_NC(n-1) - rho_NC(n-2);
    v_NC(n) = 2*v_NC(n-1) - v_NC(n-2);
    T_NC(n) = 2*T_NC(n-1) - T_NC(n-2);
    
   %Defining Mass flow rate, Pressure and Mach number
   Mass_flow_rate_NC = rho_NC.*a.*v_NC;
   Pressure_NC = rho_NC.*T_NC;
   Mach_number_NC = (v_NC./sqrt(T_NC));
   
   %Calculating Variables at Throat section
   rho_throat_NC(k) = rho_NC(throat);
   T_throat_NC(k) = T_NC(throat);
   v_throat_NC(k) = v_NC(throat);
   Mass_flow_rate_throat_NC(k) = Mass_flow_rate_NC(throat);
   Pressure_throat_NC(k) = Pressure_NC(throat);
   Mach_number_throat_NC(k) = Mach_number_NC(throat);
   
   %Plotting
   %Comparasion of non-dimensional Mass flow rates at diffrent time steps
   figure(1);
   if k == 50
       plot(x,Mass_flow_rate_NC,'m','linewidth',1.25);
       hold on
   elseif k == 100
       plot(x, Mass_flow_rate_NC,'c','linewidth',1.25);
       hold on
   elseif k == 200
       plot(x, Mass_flow_rate_NC,'g','linewidth',1.25);
       hold on
   elseif k == 400
       plot(x, Mass_flow_rate_NC,'y','linewidth',1.25);
       hold on
   elseif k == 800
       plot(x, Mass_flow_rate_NC,'k','linewidth',1.25);
       hold on
   end
   % Labelling
   title('Mass flow rates at diffrent Timesteps - NON CONSERVATIVE FORM)')
   xlabel('Distance')
   ylabel('Mass flow rate')
   legend({'50^t^h Timestep','100^t^h Timestep','200^t^h Timestep','400^t^h Timestep'})
   axis([0 3 0 2])
   grid on 
end


end

Explanation:-

• The purpose of the function is to calculate the values of Mass flow rate, Density, velocity,Temprature,Density at throat, velocity at throat, Temprature at throat, Mass flow throat, Pressure at throat, Mach number at throat by the known equations for all the grid points.
• For this initial boundary conditions are established for Density, Temprature & Velocity and profiles for calculating area & throat.
• Now, outer Time loop is computed for loop to calculate all the grid points for various formulas inside the loop.
• MacCarmack technique is used to solve the numerical solution of governing equations such as continuity,Momentum & Energy Equation.
• The calculated values are averaged in the next step and updating the solution with old values.
• Finally plotting is done for Distance(x) & Mass flow rate for different Timesteps.
• program is ended.

CONSERVATIVE FORM:-

% Function for Conservative form
function [Mass_flow_rate_C, Pressure_C, Mach_number_C, rho_C, v_C, T_C, rho_throat_C, v_throat_C, T_throat_C, Mass_flow_rate_throat_C, Pressure_throat_C, Mach_number_throat_C] = Conservative_form(x, dx, n, gamma, nt, C)
%Determining Initial profile of Density and temprature
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)
        rho_C(i) = 0.634 - 0.3507*(x(i) - 1.5);
        T_C(i) = 0.833 - 0.3507*(x(i) - 1.5);
    end
end

%Area & Throat profiles
a = 1 + 2.2*(x - 1.5).^2;
throat = find(a==1);

%Initial profile for Velocity and Pressure
v_C = 0.59./(rho_C.*a);
Pressure_C = rho_C.*T_C;

%Calculate the Initial Solution Vectors(U)
U1 = rho_C.*a;
U2 = rho_C.*a.*v_C;
U3 = rho_C.*a.*(T_C./(gamma -1)) + ((gamma/2).*v_C.^2);

% Outer Time Loop
for k = 1:nt
    
    %Time-step control
    dt = min(C.*dx./(sqrt(T_C) + v_C));
    
    %Saving a copy of old Values
    U1_old = U1;
    U2_old = U2;
    U3_old = U3;
    
    %Calculating the Flux vector (F)
    F1 = U2;
    F2 = ((U2.^2)./U1) + ((gamma - 1)/gamma)*(U3 - (gamma/2)*((U2.^2)./U1));
    F3 = ((gamma*U2.*U3)./U1) - (gamma*(gamma - 1)*((U2.^3)./(2*U1.^2)));
    
    %Predictor Method
    for j = 2:n-1
        
        %Calculating the Source term(J)
        J2(j) = (1/gamma)*rho_C(j)*T_C(j)*((a(j + 1) - a(j))/dx);
        
        %Continuity Equation
        dU1_dt_P(j) = -((F1(j+1) - F1(j))/dx);
        
        %Momentum Equation
        dU2_dt_P(j) = -((F2(j + 1) - F2(j))/dx) + J2(j);
        
        %Energy Equation
        dU3_dt_P(j) = -((F3(j + 1) - F3(j))/dx);
        
        %Updating New values
        U1(j) = U1(j) + dU1_dt_P(j)*dt;
        U2(j) = U2(j) + dU2_dt_P(j)*dt;
        U3(j) = U3(j) + dU3_dt_P(j)*dt;
    end
    
    % Updating Flux vectors
    F1 = U2;
    F2 = ((U2.^2)./U1) + ((gamma - 1)/gamma)*(U3 - ((gamma*0.5*(U2.^2))./(U1)));
    F3 = ((gamma*U2.*U3)./U1) - ((gamma*0.5*(gamma - 1)*(U2.^3))./(U1.^2));
    
    %Corrector Method
    for j = 2:n-1
        %Updating the source term(J)
        J2(j) = (1/gamma)*rho_C(j)*T_C(j)*((a(j) - a(j - 1))/dx);
        
        %Simplyfing the Equation's spatial terms
        dU1_dt_C(j) = -((F1(j) - F1(j - 1))/dx);
        dU2_dt_C(j) = -((F2(j) - F2(j - 1))/dx) + J2(j);
        dU3_dt_C(j) = -((F3(j) - F3(j - 1))/dx);
    end
    
    %Compute Average Time Derivative
    dU1_dt = 0.5*(dU1_dt_P + dU1_dt_C);
    dU2_dt = 0.5*(dU2_dt_P + dU2_dt_C);
    dU3_dt = 0.5*(dU3_dt_P + dU3_dt_C);
    
    %Final Solution Update
    for i = 2:n-1
        U1(i) = U1_old(i) + dU1_dt(i)*dt;
        U2(i) = U2_old(i) + dU2_dt(i)*dt;
        U3(i) = U3_old(i) + dU3_dt(i)*dt;
    end
    
    %Applying Boundary Conditions
    %At Inlet
    U1(1) = rho_C(1)*a(1);
    U2(1) = 2*U2(2) - U2(3);
    
    v_C(1) = U2(1)./U1(1);
    U3(1) = U1(1)*((T_C(1)/(gamma - 1)) + ((gamma/2)*(v_C(1)).^2));
    
    %At 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);
    
    %Final primitive value Update
    rho_C = U1./a;
    v_C = U2./U1;
    T_C = (gamma - 1)*((U3./U1 - (gamma/2)*(v_C).^2));
    
    % Defining Mass flow rate, Presssure and Mach number
    Mass_flow_rate_C = rho_C.*a.*v_C;
    Pressure_C = rho_C.*T_C;
    Mach_number_C = (v_C./sqrt(T_C));
    
    %Calculating Variables at throat section
    rho_throat_C(k) = rho_C(throat);
    T_throat_C(k) = T_C(throat);
    v_throat_C(k) = v_C(throat);
    
    Mass_flow_rate_throat_C(k) = Mass_flow_rate_C(throat);
    Pressure_throat_C(k) = Pressure_C(throat);
    Mach_number_throat_C(k) = Mach_number_C(throat);
    
    % Plotting
    % Comparasion of non-dimensional Mass flow rates at diffrent time steps
    figure(4)
    if k == 50
        plot(x, Mass_flow_rate_C, 'm', 'linewidth', 1.25);
        hold on
    elseif k == 100
        plot(x, Mass_flow_rate_C, 'c', 'linewidth', 1.25);
        hold on
    elseif k == 200
        plot(x, Mass_flow_rate_C, 'g', 'linewidth', 1.25);
        hold on
    elseif k == 400
        plot(x, Mass_flow_rate_C, 'y', 'linewidth', 1.25);
        hold on
    elseif k == 800
        plot(x, Mass_flow_rate_C, 'k', 'linewidth', 1.25);
        hold on
    end
    
    %Labelling
    title('Mass floe rates at diffrent Timesteps - CONSERVATIVE FORM')
    xlabel('Distance')
    ylabel('Mass flow rate')
    legend({'50^t^h Timestep', '100^t^h Timestep', '200^t^h Timestep','400^t^h Timestep','800^t^h Timestep'});
    axis([0 3 0.4 1])
    grid on
end


end

Explanation:-

• The purpose of the function in to calculate the values of Mass flow rate, Pressure, Mach number, Density, velocity, Temprature, Density throat, velocity throat, Temprature throat, mass flow rate at throat, pressure at throat, Mach number at throat by known equations for all grid points.
• Initial values of area, throat, velocity, & Pressure are calculated using formulas.
• In order to simplify and avoid errors, equation is splitted into solution vector(U), Flux Vector(F) and source term(j).
• Now, Outer time loop computed loop is computed for loop to calculate all the grid points for various formulas inside the loop.
• Boundary conditions are applied within the for loop and subsequently used for calculating mass flow rate, Mach Number & pressure and also calculating the valuest throat section.
• Finally plotting is done for Distance(x) & Mass Flow rate for different Time steps.
• program is ended.

Outputs:-

Non-Conservative form:-

Time steps variation of primitive variable at Throat area:-

Variation Flow Rate Distribution:-

Mass Flow Rate at Diffrent Time Steps:-

CONSERVATIVE-Form:-

Time steps variation of primitive variable at Throat area:-

Variation Flow Rate Distribution:

Mass Flow Rate at Diffrent Time Steps:-

Comparasion plot of Normalized Mass flow rate of Non conservative & Conservative Forms:-

Observation:-

Eventhough the final solution of both the forms are almost similar, the Non Conservative form oscillates much when compared to Conservative form. The latter also gives almost stable solution along the length. By this we can conclude that Conservative form attains stability in lesser number of iterations with earliest time.

GRID DEPENDENCE TEST:-

Observation:-

Grid dependency is important factor to determine whether the problem can be solved faster.. So here we increased the grid points from 31 to 61 and populated them in a sheet.

It is clear from the above table that the values are almost same along all the primitive Variables for both the forms even when the grid points(n) is changed from 31 to 61 showcasing that it doesn't get affected by the increase of Grid points.

So we can safely conclude that the above problem is Grid Dependent.

CONCLUSION:-

Cumulating all the above results we can evidently conclude that Conservative form converges slightly faster than Non Conservative form offering stable solution along the simulation. yet both the forms are dependable in calculating the solutions which is almost precise.