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

Super-sonic nozzle flow

We consider the steady, isentropic flow through a convergent-divergent nozzle. The pressure, temp, area, and velocity are denoted `P0,T0,V,A` at the nozzle's intel. The flow expands isentropically to supersonic speeds at the nozzle exit, where the exit pressure, temp, velocity, and Mach no are denoted by`Pe,Te,Ve,Me`.

The flow is locally subsonic in the convergent section of the nozzle, sonic at the throat, and supersonic at the divergent section. The Sonic flow(M=1) at the throat means that the local velocity at this location is equal to the local speed of sound (V=V*=a).similarly, the sonic flow values of pressure and temp are denoted by P* and T*. 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.

The governing equations in terms of nondimensional variables for this quasi-one-dimensional steady and isentropic flow can be expressed respectively as:

Non- Conservative form:-

Continuty equation:

`(∂ ρ ' )/(∂ t ') = − ρ ' ( (∂ V ') /(∂ x ') ) − ρ ' V ' (( ∂ ( ln A ' ) )/(∂ x ') ) − V ' ( (∂ ρ ' )/(∂ x ' ))`

Momentum equation:

`(∂ V ')/( ∂ t ') = − V ' ( (∂ V ') /(∂ x ') ) − ( 1 /gamma) ( (∂ T ')/( ∂ x ') + T ' ρ ' (( ∂ ρ ')/ ∂ x ') )`

Energy equation:

`(∂ T ')/( ∂ t ') = − V ' ( (∂ T ')/( ∂ x ') ) − ( gamma − 1 ) T ' [ (∂ V ')/ (∂ x ') + V ' (∂ (ln A ' ))/( ∂ x ' )]`

Conservative form:-

Continuity equation:

`(∂ ( ρ ' A ' )) /(∂ t ') +( ∂ ( ρ ' A ' V ' ) )/(∂ x ') = 0`

Momentum equation:

`(∂ ( ρ ' A ' V ' ))/( ∂ t ') + (∂ [ ρ ' A ' ( V ' ) 2 + ( 1 /gamma) ρ ' A ' ])/( ∂ x ') = ( 1 /gamma ) ρ ' (( ∂ A ' )/(∂ x ' ))`

Energy equation:

`(∂ [ ρ ' ( (e ') /(gamma − 1) + ( gamma/ 2 ) ( V ' ) ^2 ) A ' ])/( ∂ t ') + (∂ [ ρ ' ( (e ') /(gamma − 1) + ( gamma/ 2 ) ( V ' ) ^2 ) V ' A ' + ρ ' A ' V ' ]) /(∂ x ') = 0`

For air at standard condition, `gamma=1.4`
These equations are solved by the Macormack method using the time marching process

Non-Dimensional Quantities:

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

Temperature `T'=T/(To)`

Density, `rho'=rho/(rho0)`

Velocity, `V'= V/(a0)` where a0 is local sonic velocity

Time, `t'=t/(L/(a0))`

Area, `A'= A/(A ^ast)`

Macormack 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 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 of the problem.

ii. predictor step

`((∂ f )/(d t ))underset(p) =( f ( i + 1 ) − f ( i ) ) /(Deltax)`

iii. solution update

`f p = f I n i t i a l + ( (∂ f )/(∂ t )) p ⋅ Deltat`

`funderset(p)`= Flow field variables calculated after the predictor step

`funderset(Inital)` = low field variable in the initial step

`Deltat=` time step size

iv. Corrector step

`(∂ f)/(d t ) c = (f ( i ) − f ( i − 1 ))/(Delta x)`

`( (∂ f )/(d t ))underset(avg) = 0.5 ( ( (∂ f )/(d t )) p + ( (∂ f )/(d t ))c )`

v. Final solution update

`f underset(sol) = funderset( Initial)+ ( (∂ f)/( ∂ t) )underset(avg) ⋅ Δ t`

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

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

Program for Non-conservation from

Explanation

• First, we define the function of the Non-Conservation form, and we pass the Destiny, Pressure, and velocity.Temperature.Mach number and Mass parameters at initial and throat for the nozzle.
• initial boundary conditions are established for Density, Temperature & Velocity, and profiles for calculating area & throat.
• Now, we start to run the loop for timesteps.
• MacCormack technique is used to solve the numerical solution of governing equations such as Continuity, Momentum & Energy Equation.
• Predictor & Corrector methods are assigned within the for loop individually for calculating values.
• The Avg for Predictor & Corrector method computes and updates the solution with old values.
• Boundary Conditions are computed and next 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.

%Function for Non-Conservative form

function [Mass_nc,Pressure_nc, Machno_nc, rho_nc, v_nc, t_nc, rho_throat_nc, v_throat_nc, t_throat_nc, Mass_throat_nc,  Pressure_throat_nc, Machno_throat_nc] = Non_conservative_form(x, dx, n, gamma, nt, c)

%Intial 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;

if n==31
throat= find(a==1);
else
    throat= find(a==min(a));
end



%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_nc = rho_nc.*a.*v_nc;
    Pressure_nc = rho_nc.*t_nc;
    Machno_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_throat_nc(k) = Mass_nc(throat);
    Pressure_throat_nc(k) = Pressure_nc(throat);
    Machno_throat_nc(k) = Machno_nc(throat);
    
    
    
  
    %PLOTTING
    %Comparison of non-dimensional Mass flow rates at different time steps
    if n==31
    figure(1);
    if k == 50
        plot(x, Mass_nc, 'm', 'linewidth', 1.25);
        hold on
    elseif k == 100
        plot(x, Mass_nc, 'b', 'linewidth', 1.25);
        hold on
    elseif k == 200
        plot(x, Mass_nc, 'g', 'linewidth', 1.25);
        hold on
    elseif k == 400
        plot(x, Mass_nc, 'y', 'linewidth', 1.25);
        hold on
    elseif k == 800
        plot(x, Mass_nc, 'k', 'linewidth', 1.25);
        hold on
    end
    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])


end

Program for conservation from

Explanation

• First, we define the function of the Non-Conservation form, and we pass the Destiny, Pressure, and velocity.Temperature.Mach number and Mass parameters at initial and throat for the nozzle.
• initial boundary conditions are established for Density, Temperature & Velocity, and profiles for calculating area & throat.
• To simplify and avoid errors, the equation is split into Solution Vector(U), Flux Vector(F), and Source Term(J).
• Now, we start to run the loop for timesteps.
• MacCormack technique is used to solve the numerical solution of governing equations such as Continuity, Momentum & Energy Equation.
• Predictor & Corrector methods are assigned within the for loop individually for calculating values.
• The Avg for Predictor & Corrector method computes and updates the solution with old values.
• Boundary Conditions are computed and next 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.

%Function for consevation form

function [Mass_c, Pressure_c, Machno_c,  rho_c, v_c, t_c, rho_throat_c, v_throat_c, t_throat_c, Mass_throat_c, Pressure_throat_c,Machno_throat_c] = Conservative(x, dx, n, gamma, nt, c)


%intial condition for the Densiity 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 and throat
a= 1+2.2*(x-1.5).^2;

if n==31
throat= find(a==1);
else
    throat= find(a==min(a));
end


%initial velocity and Pressure
v_c = 0.59./(rho_c.*a);
Pressure_c= rho_c.*t_c;

%intial solutions of vector(u)
u1=rho_c.*a;
u2=rho_c.*a.*v_c;
u3=rho_c.*a.*((t_c./(gamma-1))+ ((gamma/2).*v_c.^2));

for k=1:nt

    

    %time step
    dt = min(c.*dx./(sqrt(t_c) + v_c));

    %update the values
    u1_old=u1;
    u2_old=u2;
    u3_old=u3;

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

        %calculate sourse term(s)
        s(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)+s(j);

        %Energy Equation
        du3_dt_p(j) = -((f3(j+1)-f3(j))/dx);


        %updating new vector 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 new flux values
    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)));

    rho_c = u1./a;
    v_c = u2./u1;
    t_c = (gamma - 1)*((u3./u1 - (gamma/2)*(v_c).^2));

    %correction method
    for j = 2:n-1

        %calculate sourse term(s)
        s2(j) = (1/gamma)*rho_c(j)*t_c(j)*((a(j) - a(j - 1))/dx);

        %Continuity Equation
        du1_dt_c(j) = -((f1(j)-f1(j-1))/dx);


        %Momentum Equation
        du2_dt_c(j) = -((f2(j)-f2(j-1))/dx)+s2(j);

        %Energy Equation
        du3_dt_c(j) = -((f3(j)-f3(j-1))/dx);


    end

    %Avg time derivatives
    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


    %Boundary conidtion

    %intlet

    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_c = rho_c.*a.*v_c;
    Pressure_c = rho_c.*t_c;
    Machno_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_throat_c(k) = Mass_c(throat);
    Pressure_throat_c(k) = Pressure_c(throat);
    Machno_throat_c(k) = Machno_c(throat);


   
    %PLOTTING
    %Comparison of non-dimensional Mass flow rates at different time steps
 if n==31 
    figure(4)
    if k == 50
        plot(x, Mass_c, 'm', 'linewidth', 1.25);
        hold on
    elseif k == 100
        plot(x, Mass_c, 'b', 'linewidth', 1.25);
    elseif k == 200
        plot(x, Mass_c, 'g', 'linewidth', 1.25);
        hold on
    elseif k == 400
        plot(x, Mass_c, 'y', 'linewidth', 1.25);
        hold on
    elseif k == 800
        plot(x, Mass_c, 'k', 'linewidth', 1.25);
        hold on
    end
end    

end
%Labelling
title('Mass flow rates at different Timesteps - CONSERVATIVE FORM')
xlabel('Distance')
ylabel('Mass flow rate')
legend({'50th timestep','100th timestep','200th timestep','400th timestep','800th timestep'})
axis([0 3 0.4 1])

end

Main function

Explanation

• First, we define the input x, dx, n, gamma, nt, and C parameters and assign the values to them.
• Now, we call the conservation and non-conversation form function and pass the input parameters.
• tic toc command is used for calculating the time elapsed between the functions.
• Plotting is done for the following solutions
  1)Timestep variation in the throat area of both forms.
  2)Variation of Flow Rate Distribution of both forms.
  3)Comparison of Normalized Mass Flow Rate of both forms.
• Each plot is computed as subplots of four primitive variables such as Density, Temperature, Pressure, and Mach Number at the throat as well as the Convergent/Divergent section.

%Main Function
close all
clear all
clc

% Initializing the domain
L=3;                        %length of Nozzle
n = 31;
x= linspace(0,L,n);
dx= x(2)-x(1);              %each grid size

nt=1400;                    %time step
c=0.5;                      %Courant Number
gamma=1.4;                  %constant number
a = 1 + 2.2*(x-1.5).^2;     %Nozzle area profile

%calling Non_conservation form function
tic;
[Mass_nc, Pressure_nc, Machno_nc, rho_nc, v_nc, t_nc, rho_throat_nc, v_throat_nc, t_throat_nc, Mass_throat_nc, Pressure_throat_nc, Machno_throat_nc] =...
    Non_conservative_form(x, dx, n, gamma, nt, c);
time_nc=toc;
fprintf('time elasped for non conservation=%d',time_nc);
hold on

%density
figure(2)
subplot(4,1,1)
plot(linspace(1,nt,nt),rho_throat_nc,'m');
ylabel('density');
legend('density at throat');
axis([0 1400 0 1.3])
title('timestep variation at throat for non conversation')

%Temperature
subplot(4,1,2)
plot(linspace(1,nt,nt),t_throat_nc,'color','b')
ylabel('Temperature')
legend({'Temperature at throat'});
axis([0 1400 0.6 1]);


%pressure
subplot(4,1,3)
plot(linspace(1,nt,nt),Pressure_throat_nc,'g');
ylabel('pressure');
legend('pressure at throat');
axis([0 1400 0.4 1.1])

%Mach number
subplot(4,1,4)
plot(linspace(1,nt,nt),Machno_throat_nc,'k');
ylabel('Mach number');
legend('Mach number at throat');
axis([0 1400 0.6 1.4])



%Steady state simulation for primitive values for non conversative

%density
figure(3)
subplot(4,1,1)
plot(x,rho_nc,'m');
ylabel('density');
legend('density ');
axis([0 3 0 1])
title(' variation of flow rate distribution for non conversation')

%Temperature
subplot(4,1,2)
plot(x,t_nc,'b');
ylabel('Temperature');
legend('Temperature');
axis([0 3 0 1])


%pressure
subplot(4,1,3)
plot(x,Pressure_nc,'g');
ylabel('pressure');
legend('pressure');
axis([0 3 0 1])

%Mach number
subplot(4,1,4)
plot(x,Machno_nc,'k');
ylabel('Mach number');
legend('Mach number');
axis([0 3 0 1])

%calling Conservation form function
tic;
[Mass_c, Pressure_c, Machno_c,  rho_c, v_c, t_c, rho_throat_c, v_throat_c, t_throat_c, Mass_throat_c, Pressure_throat_c, Machno_throat_c] = ... 
    Conservative(x, dx, n, gamma, nt, c);
time_c=toc;
fprintf('time elasped for  conservation=%d',time_c);
hold on

figure(5)
%Density
subplot(4,1,1)
plot(rho_throat_c,'color','m')
ylabel('Density')
legend({'Density at throat'});
axis([0 1400 0 1]);
title('Timestep variation at throat area - CONSERVATIVE FORM')

%Temperature
subplot(4,1,2)
plot(t_throat_c,'color','b')
ylabel('Temperature')
legend({'Temperature at throat'});
axis([0 1400 0.6 1]);

%Pressure
subplot(4,1,3)
plot(Pressure_throat_c,'color','g')
ylabel('Pressure')
legend({'Pressure at throat'});
axis([0 1400 0.2 0.8]);

%Mach Number
subplot(4,1,4)
plot(Machno_throat_c,'color','k')
xlabel('Timesteps')
ylabel('Mach Number')
legend({'Mach Number at throat'});
axis([0 1400 0.6 1.4]);


%Steady state simulation for primitve 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])
title('Variation of Flow Rate Distribution - CONSERVATIVE FORM')

%Temperature
subplot(4,1,2);
plot(x,t_c,'color','b')
ylabel('Temperature')
legend({'Temperature'});
axis([0 3 0 1])

%Pressure
subplot(4,1,3);
plot(x,Pressure_c,'color','g')
ylabel('Pressure')
legend({'Pressure'});
axis([0 3 0 1])

%Mach number
subplot(4,1,4);
plot(x,Machno_c,'color','k')
xlabel('Nozzle Length')
ylabel('Mach number')
legend({'Mach number'});
axis([0 3 0 4])


%COMPARISON PLOT of Normalized Mass Flow Rate of both forms
figure(7)
hold on
plot(x,Mass_nc,'r','linewidth',1.25)
hold on
plot(x,Mass_c,'b','linewidth',1.25)
hold on;

legend('Non Conservative','Conservative');
xlabel('Length of the Domain')
ylabel('Mass Flow Rate')
title('Comparison of Normalized Mass Flow Rate of both forms')

code for Grid independence test

close all
clear all
clc

L=3;                    %length of nozzle
n=[31 63 93];           %number of grid points
nt=1400;                %number of time steps
c = 0.5;                %CFL number
gamma=1.4;              %constant value

for i = 1:length(n)

    x=linspace(0,L,n(i));
    dx =x(2)-x(1);

    [Mass_nc, Pressure_nc, Machno_nc, rho_nc, v_nc, t_nc, rho_throat_nc, v_throat_nc, t_throat_nc, Mass_throat_nc,  Pressure_throat_nc, Machno_throat_nc] = Non_conservative_form(x, dx, n(i), gamma, nt, c);
    [Mass_c, Pressure_c, Machno_c,  rho_c, v_c, t_c, rho_throat_c, v_throat_c, t_throat_c,  Mass_throat_c, Pressure_throat_c,Machno_throat_c] = Conservative(x, dx, n(i), gamma, nt, c);



    %Grid independency test at throat
    figure(8)
    subplot (2,1,1)
    plot(Mass_throat_c,'linewidth',1.5)
    hold on
    xlabel('Time steps (nt)')
    ylabel('Mass flow rate')
    legend('@31 nodes','@63 nodes','@93 nodes')
    title('Conservative form: Grid independency test at throat point ')

    subplot (2,1,2)
    plot(Mass_throat_nc,'linewidth',1.5)
    hold on
    xlabel('Time steps (nt)')
    ylabel('Mass flow rate')
    legend('@31 nodes','@63 nodes','@93 nodes')
    title('Non-Conservative form: Grid independency test at throat point ')


    %Grid independency test along the nozzle
    figure(9)
    subplot (2,1,1)
    plot(x,Mass_c,'linewidth',1.5)
    hold on
    xlabel('Length of nozzle')
    ylabel('Mass flow rate')
    legend('@31 nodes','@63 nodes','@93 nodes')
    title('Conservative form: Grid independency test along the nozzle')

    subplot (2,1,2)
    plot(x,Mass_nc,'linewidth',1.5)
    hold on
    xlabel('Length of nozzle')
    ylabel('Mass flow rate')
    legend('@31 nodes','@63 nodes','@93 nodes')
    title('Non-Conservative form: Grid independency test along the nozzle')
    axis([0 3 0.580 0.595])
end

Output Graphs:

Time steps variation of primitive variable at Throat area:

Non-Conservative Form

Conservative Form

Observation:

When we compare both graphs we can conclude that Non-conservation from initial variations is observed. That means it consumes a higher amount of time to achieve a stable solution when compared to the Conservative form.

Variation of Flow Rate Distribution

Non Conservative form

Conservative form

Observation:

When we compare both graphs we can conclude that, when the Mach number increases the other's parameters pressure, density, and temperature decrease in both forms.

Mass flow rates at different Timesteps:

Non Conservative form:

Conservative form:

Observation:

When we compare both graphs we can conclude 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 Normalized Mass Flow Rate of both NON-CONSERVATIVE & CONSERVATIVE FORMS:

Observation:

The Non-Conservative form oscillates much when compared to the Conservative form and also gives an almost stable solution along the length. And we can say that the Conservative form attains stability in less number of iterations and time.

Grid independency test for both Non-conservative form & Conservative form:

Observation:

From figure 8 we can conclude that Non-conservation from initial variations is observed. That means it consumes a higher amount of time to achieve a stable solution when compared to the Conservative form. And when grid size increases initial variations in a mass flow rate decrease in both form.

From fig(8) and fig(9), we can conclude that the plot is a slight change in the throat area at a different grid Size.