Week 7 - 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 MacCormark's technique using Mathlab.

Problem statement :

Determine the steady stste temperature distribution for the flow field variables and to investigate the difference between the two forms of governing equation by comparing their solutions along with the plots stated below.

1. Steady state distribution of primitive variable inside the nozzle.

1. Time wise variation of the primitive variable.

1. Variation of mass rate distribution inside the nozzle at different time steps during the time marching process.

1. Comparision of Normalised mass flow rate distribution of both forms.

Theory :

Any numerical solution of isentropic flow through a quasi 1D subsonic-supersonic nozzle is unreasonable, while analytic solution exits for these kind of problems. However we are trying here to establish the broadened capability of CFD in various streams rather in representing the approximation form of analytic solution. Here, we implement the numerical solution with Conservative and Non-Conservative forms.

Governing equations are marches in terms of PDE for both forms in order to use time marching solution for the problrm. These forms represent in numerical aspect only, if we analytic all the methods are same.

Conservative form :

Consider an infinitesimal small fluid element inside nozzle, and the fluid moving that elements, such type of condition is conservative form and the equation obtained in such cases are called as conservative form of governing equation.

Continuity equation : `partial((rho^'A'))/(partialt^')+partial((rho^'A^'V^'))/(partialx^') = 0`

Momemtum equation : `partial((rho^'A^'V^'))/(partialt^')+partial((rho^'A^'V^('2)+1/gamma rho^'A^'))/(partialx^') = 1/gammap^'(partialA^')/(partialx^')`

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

Non conservative form :

Consider an infinitesimal small fluid element inside the nozzle, the fluid same fluid element moving along streamline, such type of condition is non conservative form and the equation obtained in such cases are called as non conservative form of governing equation.

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

Non dimensional quantities :

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

Temperature : `T^' = T/T_0`

Density : `rho^' = rho/rho_0`

Non dimensional distance : `x^' = x/L`

Non dimensional velocity : `V^' = V/a_0` where `a_0` is the local sonic velocity.

Non dimensional time : `t^' = t/L/a_0` were `L/a_0` is characterstic time scale.

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.

Macormack method :

To discretize the numerical solution we use the technique called MacCarmack method. The MacCormack method is a variation of the two step Lax-Wendroff schemebut is much simpler in application.

The application of MacCormack method proceeds in two steps, a predictor step which is followed by a corrector step.

• Predictor step : In the predictor step, a 'provisional value `u` at time level `n+1`(denoted by `u_i^(n+1)`) is estimated obtained by replacing the spatial and temporal derivatives in the previous first order equation using forward difference scheme `Deltat` is used for timesteps.
• Corrector step : In the corrector step, the predicted value `u_i^(n+1)` is corrected according to the equation. Note that the corrector step uses backword differencing scheme for spatial derivative. The time step used in the corrector step is `Deltat/2`in contrast to the `Deltat`used in the predictor step.
  
  

Main program of solving Quasi 1D subsonic-supersonic nozzle :

clear all
close all
clc

% solving the QUASI 1D supersonic nozzle flow equation
% Input parameters
n = 31;
gamma = 1.4;
x = linspace(0,3,n);
dx = x(2) - x(1);

% No. of time steps
nt = 1400;

% Courant number
C = 0.5;

% Function for non conservative 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 Elasped for Non Conservative form: %0.4f',Time_elapsed_NC);
hold on

% Plotting
% Timestep variation of properties at throat area in non comservative 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.3]);
grid minor;
title('Timestep variation at throat area - NON CONSERVATIVE FORM')

% Temperature
subplot(4,1,2)
plot(linspace(1,nt,nt),T_throat_NC,'color','c')
ylabel('Temperature')
legend('Temperature 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','r')
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(3);

% Density
subplot(4,1,1);
plot(x,rho_NC,'color','m')
ylabel('Density')
legend('Density');
axis([0 3 0 1]);
grid minor;
title('Variation of flow rate distribution - NON CONSERVATIVE FORM')

% Temperature
subplot(4,1,2);
plot(x,T_NC,'color','c')
ylabel('Temperature')
legend('Temperature');
axis([0 3 0 1]);
grid minor;

% Pressuer
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','r')
xlabel('Nozzle lenght')
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_elapsed_C = toc;
fprintf('Time Elasped for conservative form:%0.4f', Time_elapsed_C);

% Plotting
% Time steps variation of properties at throat area in Conservative form
figure(5)

% 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('Timesteps variation at throat area - CONSERVATIVE FORM')

% Temperature
subplot(4,1,2)
plot(T_throat_C,'color','c')
ylabel('Temperature')
legend('Temperature 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','r')
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')

% Temperature
subplot(4,1,2)
plot(x,T_C,'color','c')
ylabel('Temperature')
legend('Temperature');
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','r')
ylabel('Nozzle length')
legend('Mach number');
axis([0 3 0 4]);
grid minor;

% Comparision plot of normalised 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('comparision of normalised mass flow rate of both forms')

Explanation :

• Code is started with clear all, close, 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 conservation and conservation forms are created as seperate functions with the neccessary inputs which is detailed below this section.
• Plotting is done for the following solutions

1. Timesteps variation at throat are of both forms.

1. Variation of flow rate distribution of both forms.

1. Comparision of normalized mass flow rate of both forms.

• Each plot is computed as subplots of four primitives variables such as density, temperature, pressure and mach number at throat as well as convergent/divergent section.

Non conservative form code :

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));

    % 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 eqution spatialterms
        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;

        % Continuity 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

    % Correct method
    for j = 2:n-1

        % Simplifying the equations 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;

        % Continuity 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 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);

    % 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 sections
    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 
    % Comparision of non dimensional mass flow rates at different timesteps
    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 rate at different 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','800^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, pressure, mach number, density, velocity, temperature, density at throat, velocity at throat, temperature at throat, mass flow rate at throat, pressure at throat, mach number at throat by the known equations for all the grid points.

• For the initial boundary conditions are established for density, temperature and velocity and profiles for calculating area and throat.
• Now outer timeloop is computed for the 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 and energy equation.
• The technique impliestwo methods for solving namely predictor method followed by corrector method.
• Predictor and corrector method are assigned within the foe loop individually for calculating values.

• The calculated values are averaged in the next step and updating the solution with old values.
• Boundary conditions are applied within the for loop and subsequently used for calculating the mass flow rate, mach number and pressure and also calculating the values throat section.
• Finally plotting is done for distance(x) and mass flow rate different timesteps.
• Then the program is ended.

Conservative form code :

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 temperature
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.3879*(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 
        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
        J2(j)  = (1/gamma)*rho_C(j)*T_C(j)*((a(j) - a(j - 1))/dx);

        % Simplifying the equation 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, pressure 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 throats 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
    % Comparision of non-dimensional mass flow rates at different 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 flow rate at different 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 is to calculate the values of mass flow rate, pressure, mach number, density, velocity, temperature, density at throat, velocity at throat, temperature at throat, mass flow rate at throat, pressure at throat, mach number at throat by the known equations foe all the grid points. To determine the initial profile of density and temperature, for loop is used inputing the knows conditions within an else-if command.
• Initial values of area, throat, velocity and pressue 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 is computed for the forloop to calculated all the grid pointd for various formulas inside the loop.

• MacCarmack techniques is used to solve the numerical solution of governing equation such s continuity, momentum and energy equation.
• The technique implies two methods namely predictor method followed by corrector method.
• Predictor and corrector method are assigned within the for loop individually for calculating values.
• The calculated values are averaged in the next step and updating the solution with old values.
• Boundary conditions are applied within the forloop and subsequently used for calculating the mass flow rate, mach number and pressuer and also calculating the values throat section.

• Finally plotting is done for distance(x) and mass flow rate for different timesteps.
• Then the program is ended.

Outputs :

Time steps variation of primitive variable at throat area :

Non-conservative form :

Conservative form :

Observation :

From the above plots we can conclude that initial variations observed in Non conservative form for a certain period of time makes it consume higher amount of time to achieve stable solution when compared to Conservative form.

Variation of flow rate distribution :

Non conservative form :

 Conservative form :

Observation :

Here by observing the above plots we could see that both the forms are in same phase giving the extra same solution for all the values.

Mass flow rate at different timesteps :

Non conservative form :

Conservative form :

Observation :

Above plots shows that mass flow rate at different timesteps along th distance oscillates initally, however the solution becomes stable after the 400th timestep for both the forms attaining a constant value.

Comparision plot of normalised mass flow rate of both non conservative and conservative forms :

Observation :

Eventhough the final solution of both the forms are almost similar, the Non conservation form oscillates much when compared to conservation form. Then latter it also gives almost stabke solution along the length. By this we can say that Conservative form attains stability in less number of iterations with earliest time.

Workspace :

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 simulatio. Yet both forms are dependable in calculating the solutions which is almost precies.