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

AIM:

• Simulate the isentropic flow through a quasi 1D subsonic-supersonic nozzle.
• Derive both the conservation and non-conservation forms of the governing equations and solve them using the MacCormack’s technique.
• Determine the steady-state temperature distribution for the flow-field variables and investigate the difference between the two forms of governing equations by comparing their solutions.
• Show the following plots inside your report
• Steady-state distribution of primitive variables inside the nozzle.
• Time-wise variation of the primitive variables.
• Variation of Mass flow rate distribution inside the nozzle at different time steps during the time-marching process.
• Comparison of Normalised mass flow rate distributions of both forms.
  
  

Expected Results:

• Perform a grid independence test on the solution and explain all your findings.
• The code must be modular. You need to write separate functions to solve both forms.
• Discuss all the plots you are providing with observed inferences.  

Theory:

Consider the study of isentropic flow through a convergent divergent nozzle. The flow at the inlet to the nozzle comes from the reservoir. The cross-section area of the reservoir is very large, and hence the velocity is very small. The flow is locally subsonic at the convergent section of the nozzle, sonic at the throat, and supersonic at the divergent section.

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

Let us consider a small control volume of the nozzle as shown below.

We assume that at a given section, where the cross-section area is A, the flow properties are uniform across the section. Hence, although the area of the nozzle changes as a function of distance along the nozzle, X, and therefore in reality the flow field is two-dimensional (the flow varies in the two-dimensional x-y space), we make the assumption that the flow properties vary only with X; this is tantamount to assuming uniform flow properties across any given cross-section. Such flow is defined as quasi-one-dimensional flow.

The governing continuity, momentum, and energy equations in terms of non-dimensional variables for this quasi-one-dimensional steady and isentropic flow can be expressed respectively as:

These equations are solved by MacCormack method using time marching process.

Non Dimensional Quantities:

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

MacCormack method:

      Let f(x) is the equation to be solved.Then the following steps are followed

I. Initial profile

      Before the calculation, the initial profile for all the flow field variables which are velocity, density and temperature has to be determined. These are guess values but should be specified with more awareness on the problem.

II. Predictor step

III. Solution update

IV. Corrector step

V. Final solution update

f_(sol) = Final solution of the flow field variables for the first time step.

Main code for solving Quasi 1D subsonic-supersonic nozzle:

% SOLVING QUASI 1D SUPERSONIC NOZZLE FLOW EQUATIONS
clear all
close all
clc

% INPUT PARAMETERS
% No. of grid points    
L = 3; %length of nozzle
n = 31;             
x = linspace(0,L,n);
dx = x(2) - x(1);
gamma = 1.4;

% No. of time step
nt = 900;  

% Courant Number(CFL)
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('\n Grid number n=%d,Time Elapsed for Non Conservative form: %0.4f', n,Time_elapsed_NC);
hold on

%PLOTTING 
%Timestep 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','linewidth',1)
ylabel('Density')
legend({'Density at throat'});
axis([0 900 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','r','linewidth',1)
ylabel('Temperature')
legend({'Temperature at THROAT'});
axis([0 900 0.6 1]);
grid minor;

%Pressure
subplot(4,1,3)
plot(linspace(1,nt,nt),Pressure_throat_NC,'color','g','linewidth',1)
ylabel('Pressure')
legend({'Pressure at THROAT'});
axis([0 900 0.4 1.1]);
grid minor

%Mach Number
subplot(4,1,4)
plot(linspace(1,nt,nt),Mach_number_throat_NC,'color','y','linewidth',1)
xlabel('Timesteps')
ylabel('Mach Number')
legend({'Mach Number at THROAT'});
axis([0 900 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','linewidth',1)
ylabel('Density')
legend({'Density'});
axis([0 3 0 1])
grid minor;
title('Variation of Flow Rate Distribution - NON CONSERVATIVE')
fprintf('\n Grid number n=%d,Density at Non Conservative form: %0.4f', n,rho_NC(n));

%Temperature
subplot(4,1,2);
plot(x,T_NC,'color','r','linewidth',1)
ylabel('Temperature')
legend({'Temperature'});
axis([0 3 0 1])
grid minor;
fprintf('\n Grid number n=%d,Temperature at Non Conservative form: %0.4f', n,T_NC(n));

%Pressure
subplot(4,1,3);
plot(x,Pressure_NC,'color','g','linewidth',1)
ylabel('Pressure')
legend({'Pressure'});
axis([0 3 0 1])
grid minor;
fprintf('\n Grid number n=%d,Pressure at Non Conservative form: %0.4f', n,Pressure_NC(n));

%Mach number
subplot(4,1,4);
plot(x,Mach_number_NC,'color','y','linewidth',1)
xlabel('Nozzle Length')
ylabel('Mach number')
legend({'Mach number'});
axis([0 3 0 4])
grid minor;
fprintf('\n Grid number n=%d,Mach number at Non Conservative form: %0.4f', n,Mach_number_NC(n));

%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('\n Grid number n=%d,Time Elapsed for Conservative form: %0.4f',n,Time_elapsed_C);

%PLOTTING 
%Timestep variation of properties at throat area in Conservative form
figure(5)

%Density
subplot(4,1,1)
plot(rho_throat_C,'color','m','linewidth',1)
ylabel('Density')
legend({'Density at throat'});
axis([0 900 0 1]);
grid minor;
title('Timestep variation at THROAT AREA-CONSERVATIVE FORM')

%Temperature
subplot(4,1,2)
plot(T_throat_C,'color','r','linewidth',1)
ylabel('Temperature')
legend({'Temperature at THROAT'});
axis([0 900 0.6 1]);
grid minor;

%Pressure
subplot(4,1,3)
plot(Pressure_throat_C,'color','g','linewidth',1)
ylabel('Pressure')
legend({'Pressure at THROAT'});
axis([0 900 0.2 0.8]);
grid minor

%Mach Number
subplot(4,1,4)
plot(Mach_number_throat_C,'color','y','linewidth',1)
xlabel('Timesteps')
ylabel('Mach Number')
legend({'Mach Number at THROAT'});
axis([0 900 0.6 1.4]);
grid minor;


%Steady state simulation for primitive values for Conservative form
figure(6);

%Density 
subplot(4,1,1);
plot(x,rho_C,'color','m','linewidth',1)
ylabel('Density')
legend({'Density'});
axis([0 3 0 1])
grid minor;
title('Variation of Flow Rate Distribution - CONSERVATIVE FORM')
fprintf('\n Grid number n=%d,Density at Conservative form: %0.4f',n,rho_C(n));


%Temperature 
subplot(4,1,2);
plot(x,T_C,'color','r','linewidth',1)
ylabel('Temperature')
legend({'Temperature'});
axis([0 3 0 1])
grid minor;
fprintf('\n Grid number n=%d,Temperature at Conservative form: %0.4f',n,T_C(n));

%Pressure
subplot(4,1,3);
plot(x,Pressure_C,'color','g','linewidth',1)
ylabel('Pressure')
legend({'Pressure'});
axis([0 3 0 1])
grid minor;
fprintf('\n Grid number n=%d,Pressure at Conservative form: %0.4f',n,Pressure_C(n));

%Mach number
subplot(4,1,4);
plot(x,Mach_number_C,'color','y','linewidth',1)
xlabel('Nozzle Length')
ylabel('Mach number')
legend({'Mach number'});
axis([0 3 0 4])
grid minor;
fprintf('\n Grid number n=%d,Mach number at Conservative form: %0.4f \n',n,Mach_number_C(n));

%COMPARISON PLOT of Normalised Mass Flow Rate of both forms
figure(7)
hold on
plot(x,Mass_flow_rate_NC,'b' ,'linewidth',1.5)
hold on;
plot(x,Mass_flow_rate_C, 'm','linewidth',1.5)
hold on;
grid on

legend('Non Conservative','Conservative');
xlabel('Length of the Domain')
ylabel('Mass Flow Rate')
title('Comparison of Normalised Mass Flow Rate of Both Forms')

Code Explanation:

• Code starts with clear all, close all and clc to enable a fresh start without the manipulation of previous values.
• Input parameters are assigned for the variables L, x, dx, n, gamma, nt and C.
• Non Conservative and Conservative forms are created as separate functions with the necessary inputs which is detailed below this section.
• tic and toc command is used for calculating the time elapsed between the functions.
• Each plot is computed as subplots of four primitive variables such as Density, Temperature, Pressure and Mach Number at throat as well as Convergent/Divergent section.
• Plotting is done for the following solutions
  1)Timestep variation at throat area of both forms
  2)Variation of Flow Rate Distribution of both forms
  3)Comparison of Normalised Mass Flow Rate of both forms

FUNCTION CODE FOR 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)

           

% calculate initial profile 
rho_NC = 1-0.3146 * x;
T_NC = 1 - 0.2341*x; % T = temperature
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 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;
        
        %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
    
    
    %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;
        
        %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 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
    %Comparison of non-dimensional Mass flow rates at different time steps
    figure(1);
    if k == 50
       plot(x, Mass_flow_rate_NC, 'm', 'linewidth', 2);
       hold on
    elseif k == 100
       plot(x, Mass_flow_rate_NC, 'c', 'linewidth', 2);
       hold on
    elseif k == 200
       plot(x, Mass_flow_rate_NC, 'g', 'linewidth', 2);
       hold on
    elseif k == 400
       plot(x, Mass_flow_rate_NC, 'y', 'linewidth', 2);
       hold on
    elseif k == 800
       plot(x, Mass_flow_rate_NC, 'k', 'linewidth', 2);
       hold on
    end
        
end
    %Labelling
    title('Mass flow rates 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

Function code 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 this initial boundary conditions are established for Density, Temperature & 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.
• MacCormack technique is used to solve the numerical solution of governing equations such as Continuity, Momentum & Energy Equation.
• The technique implies two methods for solving namely Predictor method followed by Corrector method.
• Predictor & 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 for loop and subsequently used for calculating the Mass Flow Rate, Mach Number & Pressure and also calculating the values throat section.
• Finally plotting is done for Distance(x) & Mass Flow Rate for different Timesteps.
• Program is ended.

FUNCTION CODE FOR 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 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(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);

        %Simplifying 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, 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 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
    %Comparison of non-dimensional Mass flow rates at different time steps
    figure(4)
    if k == 50
       plot(x, Mass_flow_rate_C, 'm', 'linewidth', 2);
       hold on
    elseif k == 100
       plot(x, Mass_flow_rate_C, 'c', 'linewidth', 2);
       hold on
    elseif k == 200
       plot(x, Mass_flow_rate_C, 'g', 'linewidth', 2);
       hold on
    elseif k == 400
       plot(x, Mass_flow_rate_C, 'y', 'linewidth', 2);
       hold on
    elseif k == 800
       plot(x, Mass_flow_rate_C, 'k', 'linewidth', 2);
       hold on
    end
               
end
    %Labelling
    title('Mass flow rates 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

Function code 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.
  To determine the initial profile of Density & Temperature, for loop is used inputing the known conditions within an else-if command.
• Initial values of area, throat, velocity & Pressure are calculated using formulas.
• In order to simplify and avoid errors, equation is split into Solution Vector(U), Flux Vector(F) and Source Term(J).
• Now Outer Time loop is computed for loop to calculate all the grid points for various formulas inside the loop.
• MacCormack technique is used to solve the numerical solution of governing equations such as Continuity, Momentum & Energy Equation.
• The technique implies two methods for solving namely Predictor method followed by Corrector method.
• Predictor & 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 for loop and subsequently used for calculating the Mass Flow Rate, Mach Number & Pressure and also calculating the values throat section.
• Finally plotting is done for Distance(x) & Mass Flow Rate for different Timesteps.
• Program is ended.

Results:

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 exact same solution for all the values.

Mass flow rates at different Timesteps:

Non Conservative form:

Conservative Form:

Observation:

Above plots shows that Mass Flow Rate at different Time steps along the Distance oscillates initially, however the solution becomes stable after the 400th Timestep for both the forms attaining a constant value.

COMPARISON PLOT of Normalised Mass Flow Rate of both NON CONSERVATIVE & CONSERVATIVE FORMS:

Observation:

Even though 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 minimum number of iterations with earliest time.

GRID DEPENDENCE TEST:

Observation:

Grid dependency is an 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 Dependant.

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.

we can conclude that the values are not much changing when the number of grid points are doubled up. Although there is a slight change at the throat area alone in the conservative form, which is not considerably significant.