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

CHALLENGE REPORT :

AIM :

To simulate the isentropic flow through a quasi 1D subsonic-supersonic nozzle and derive both the conservation and non-conservation forms of the governing equations and sovle them using the MacCormack's technique.

OBJECTIVE :

* To simulate the isentropic flow through a quasi 1D subsonic-supersonic nozzle and

* Derive both the conservation and non-conservation forms of the governing equations and sovle them using the MacCormack's technique.

* Plot the Steady-state distribution of primitive variables inside the nozzle.

* Plot the time-wise variation of the primitive variables.

* Plot Variation of Mass flow rate distribution inside the nozzle at different time steps during the time-marching process.

* Plot Comparison of Normalized mass flow rate distributions of both forms.

THEORY :

Nozzle diagram :

The nozzle has three main sections

1. convergent.

2. Throat.

3. Divergent

• At convergent section the flow velocity is less than the sound velocity therefore there mach number is less than 1.
• At throat section the flow velocity is equal to the sound velocity therefore there mach number is equal to 1.
• At divergent section the flow velocity is greater than the sound velocity therefore there mach number is greater than 1.

So here, the flow at the inlet of the nozzle comes from a reservoir where pressure and temperature are denoted as byrho_0 and T_0 respectively. The cross sectional area of the reservoir is large (theoretically, A->oo) and hence the velocity is very small (V->0). Thus, rho_0 and T_0 are the stagnation values or total pressure and temperature respectively. The flow expands isentropically to supersonic speeds at the nozzle exits, where the exit pressure, temperature, velocity, and Mach number are denoted by P_e, T_e, V_e and M_e respectively. The flow is locally subsonic in the convergent section, sonic at the throat section 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. Using an asterisk to denote sonic flow values, we have at the throat V = V* = a* Similarly, the sonic flow values of pressure and temperature are denoted by P* and T*, respectively.

Although the area of nozzle changes as a function of distance along the nozzle x and therefore in reality the flow field is two dimensional (xy), we make an assumption that the flow properties vary only with x, this tendency to assume uniform flow properties across any given cross-section. Such flow is defined as Quasi-one-dimensional flow.

In order to solve quasi-one-dimensional problem, we need PDE or Integral form of equations for CFD analysis. These PDE or Integral form of the governing equations are in both conservative and non-conservative form.

Before, writing the governing equations in conservative and non-conservaticve form, we need to non-dimensionalize the primitive variables for removing the constraints of dimensions to analyse the numerical solution. It will certainly ease the method of analysis for a given problem.

Mach number :

It is a dimensionless quantity representing the ratio of flow velocity to the local speed of the sound.

M = u/c

M - Mach number.

u - Flow velocity.

c - Velocity of sound.

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 as a_0 is the local sonic velocity (a_0 = sqrt(gammaRT_0))

Non-dimensional time t' = t/((L/a_0)) where L/a_0 is characteristic time scale.

Non dimensional area A'=(A)/(A^(**)).

The above is the generic form of non-dimensional conservation form of the continuity, momentum and energy equations for quasi-one-dimensional flow. Let us define the above equations in a simple generic form, which will be used in calculations instead of the above rigid and complex form. Elements of the solutions vector U, the flux vector F, and the source term J as follows,

Non Conservative form :

MacCormack method :

MacCormack's method is an explicit finite-difference technique which is second order accurate in both space and time. MacCormack method generally uses three step approach(Predictor, Corrector and averaging) to calculate a final numerical solution for an iteration. In predictor step, forward difference is taken for spatial derivatives whereas in corrector step, rearward (backward) difference is taken for spatial derivatives. Furthermore, the solutions obtained from these two steps are averaged and final numerical solution is calculated for a single iteration.

Initial and boundary conditions :

To start the time-marching calculations, we must stipulate initial conditions for rho, T and V as a function of x; that is, we must set up values of rho, T and V at time t = 0.

Non-conservative form :

Initial conditions :

rho = 1-0.3146x

T = 1-0.2314x

V = (0.1+1.09x).T^(1/2)

Boundary conditions :

In our case, we are having a subsonic flow in the convergent section (M<1), choked flow at the throat section (M=1) and supersonic outflow in the divergent section (M>1). So by using method of characteristics for an unsteady, one-dimensional flow we can precisely allot inflow and outflow boundary conditions and can identify which quantities to state at the boundaries.

What does these characteristics curve and streamline direction tells us?

• If any characteristic line is entering into the domain, implies that the value of one dependent flow field variable must be specified at that boundary.
• However, if any characteristic line is exiting the domain, implies that the value of the another dependent flow field variable must be allowed to float at that boundary i.e. it must be calculated in steps of time as function time wise solution of the flow field.
• Similarly, if streamline flows into the domain then specify third dependent flow field variable and if streamline flows out of the domain then third dependent flow field variable must be allowed to float.

In our case, at i = 1(inlet), we have primitive varibles like rho, A, V and T. So which variables must be taken for initial and boundary conditions? Here, we know the rho and T at node 1 are equal to the reservoir's rho and T which is assumed to be rho=1 and T=1, so it becomes,

Inflow boundary condition :

rho(1)=1and T(1) =1(Since it is Dirichlet boundary condition)

However, velocity and area tends to infinity, so for velocity, we need to calculate it from the information obtained from inner nodes,

V(1)=2⋅V(2)-V(3)">V(1) = 2V(2) - V(3) (It is Neumann boundary condition)

Outflow boundary condition :

At the outflow boundary condition (Node =n), as we can see that, both the characteristics lines are out of the domain and streamline also flows out of the domain. This is because, at the outflow condition, the flow is supersonic (M>1), which is certainly greater than the speed of sound. So in this condition, information transport is very fast. Therefore, there are no flow field variables which require their values to be stipulated at the supersonic outflow boundary; all variables must be allowed to float at this boundary.

rho(n) = 2rho(n-1) - rho(n-2)

V(n) = 2V(n-1) - V(n-2)

T(n) = 2T(n-1) - T(n-2)

The supersonic outflow boundary is purely Neumann boundary condition, which is extrapolating the boundary nodes using the information from inner nodes.

Conservative form :

In conservative form, we will initialize the nozzle profile rationally. As in conservative form, seen above from the governing equations, we have solution vectors, flux terms and a source term. So it is the prior matter of concern to state initial and boundary conditions in those terms rather than that of primitive variables.

Initial conditions :

These variations are slightly more realistic than those assumed in the nozzle shape and initial conditions in non-conservative form of governing equations discussed above. This is an anticipation that the stability behaviour of the finite-difference formulation using the conservation form of the governing equations might be slightly more sensitive, and therefore it is useful to start with more improved initial conditions than those discussed above in non-conservative form.

Boundary conditions :

In conservative form, by referring fig (b) and absorbing the same method of characteristics approach we can allot boundary conditions for inflow and outflow boundary conditions. However, in this case, the boundary conditions will be in terms of solution vectors. Similarly as discussed above, at the subsonic inflow (M<1), we need to specify two quantities at the initial node, where we have rho',A' and U2 = rho'A'V' = 0.59(Mass flow rate) specified and one quantity will be allowed to float by referring left hand characteristic line. However, at the supersonic outflow (M>1), the information is travelling exceeds the speed of sound, hence it is irrational to specify outflow boundary condition. Therefore, all the dependent variables are allowed to float and they will be calculated in steps of time as a function time-wise solution of the flow field.

Inflow boundary condition :

The inflow boundary condition can be given as,

U1 = rho'(1)A'(1) (Dirichlet boundary condition)

U2 = 2.U2(2) - U2(3)

V(1)= (U2(1))/(U1(1)

Note: This velocity at the first node will be substituted below in calculating U3

U3 = U1(1)'((T'(1))/(gamma-1)+gamma/2(V'(1))^2) (Neumann boundary condition)

Outflow boundary condition :

All the dependent variables will be allowed to float. Hence, the outflow boundary can be given as,

U1(n) = 2U1(n-1) - U1(n-2) (Neumann boundary condition)

U2(n) = 2U2(n-1) - U2(n-2) (Neumann boundary condition)

U3(n) = 2U3(n-1) - U3(n-2) (Neumann boundary condition)

Time step calculation :

The governing system of equations is hyperbolic with respect to time. So there is stability constraint on this system which is known as Courant-Friedrich-Lewy (CFL) condition. The Courant number is a simple stability analysis of a linear hyperbolic equations that gives result CFL<=1 for an explicit numerical solution to be stable.

CFL condition can be given as,

CFL condition = (any quantity)*(time-step with which information would travel)/(element step size)

So, in our case, the quantity which is going to travel is velocity, or in other sense velocity is the quantity which will be travelling with certain time from one node to other. So here velocity of the fluid at a point in the flow is coupled with local speed of sound. The quantity can be expressed to be as a+v'

CFL = ("a + v').(Deltat)/(Deltax)

We will assume the Courant number as C = 0.5, accordingly will calculate the Deltat. In addition to it, we will calculate (Deltat)_i^t

at all the grid points, i = 1 to i = n, and then we'll choose the minimum value of Deltat in our calculations.

After considering all these governing equations and their respective pre-requisites to solve the problem, then we will move towards to the MATLAB code for execution.

PROGRAM CODE IN MATLAB :

MAIN CODE :

clear all
close all
clc
%inputs
n=31;
x=linspace(0,3,n)
dx=x(2)-x(1);
gamma=1.4;

%time step
nt=1400;
% dt=0.001;

% function for non conservative form
[mass_flow_nc, pressure_nc, Mach_nc, rho, v, t, rho_throat_nc, v_throat_nc, t_throat_nc, mass_flow_throat_nc, pressure_throat_nc, Mach_throat_nc] = non_conservative_form(x,dx,n,gamma,nt)
hold on

%plot for Steady-state distribution of primitive variables inside the nozzle

figure(2)
subplot(5,1,1)
plot(x,rho,'color','r','linewidth',2)
xlabel('length of the nozzle , x')
ylabel('density')
grid on
title('Steady-state distribution of primitive variables inside the nozzle- non conservative form')
subplot(5,1,2)
plot(x,v,'color','b','linewidth',2)
xlabel('length of the nozzle , x')
ylabel('velocity')
subplot(5,1,3)
plot(x,t,'color','g','linewidth',2)
xlabel('length of the nozzle , x')
ylabel('Temperature')
subplot(5,1,4)
plot(x,pressure_nc,'color','r','linewidth',2)
xlabel('length of the nozzle , x')
ylabel('Pressure')
subplot(5,1,5)
plot(x,Mach_nc,'color','r','linewidth',2)
xlabel('length of the nozzle , x')
ylabel('Mach no')

% plot for Time-wise variation of the primitive variables

figure(3)
subplot(4,1,1)
plot(rho_throat_nc,'color','r','linewidth',2)
xlabel('length of the nozzle , x')
ylabel('density')
grid on
title(' Time-wise variation of the primitive variables- non conservative form')
subplot(4,1,2)
plot(v_throat_nc,'color','b','linewidth',2)
xlabel('length of the nozzle , x')
ylabel('velocity at throat')
subplot(4,1,3)
plot(t_throat_nc,'color','g','linewidth',2)
xlabel('length of the nozzle , x')
ylabel('Temperature at throat')
subplot(4,1,4)
plot(pressure_throat_nc,'color','r','linewidth',2)
xlabel('length of the nozzle , x')
ylabel('Pressure at throat')
% subplot(5,1,5)
% plot(Mach_nc,'color','r','linewidth',2)
% xlabel('length of the nozzle , x')
% ylabel('Mach no')
hold on



% function for conservative form
[mass_flow_c, pressure_c, Mach_c, rho_c, v_c, t_c, rho_throat_c, v_throat_c, t_throat_c, mass_flow_throat_c, pressure_throat_c, Mach_throat_c] = conservative_form(x,dx,n,gamma,nt)


%plot for Steady-state distribution of primitive variables inside the
%nozzle for conservative form

figure(5)
subplot(5,1,1)
plot(x,rho_c,'color','r','linewidth',2)
xlabel('length of the nozzle , x')
ylabel('density')
grid on
title('Steady-state distribution of primitive variables inside the nozzle- conservative form')
subplot(5,1,2)
plot(x,v_c,'color','b','linewidth',2)
xlabel('length of the nozzle , x')
ylabel('velocity')
subplot(5,1,3)
plot(x,t_c,'color','g','linewidth',2)
xlabel('length of the nozzle , x')
ylabel('Temperature')
subplot(5,1,4)
plot(x,pressure_c,'color','r','linewidth',2)
xlabel('length of the nozzle , x')
ylabel('Pressure')
subplot(5,1,5)
plot(x,Mach_c,'color','r','linewidth',2)
xlabel('length of the nozzle , x')
ylabel('Mach no')

% plot for Time-wise variation of the primitive variables

figure(6)
subplot(4,1,1)
plot(rho_throat_c,'color','r','linewidth',2)
xlabel('length of the nozzle , x')
ylabel('density at throat')
grid on
title(' Time-wise variation of the primitive variables-conservative form')
subplot(4,1,2)
plot(v_throat_c,'color','b','linewidth',2)
xlabel('length of the nozzle , x')
ylabel('velocity at throat')
subplot(4,1,3)
plot(t_throat_c,'color','g','linewidth',2)
xlabel('length of the nozzle , x')
ylabel('Temperature at throat')
subplot(4,1,4)
plot(pressure_throat_c,'color','r','linewidth',2)
xlabel('length of the nozzle , x')
ylabel('Pressure at throat')


figure(7)
plot(x,mass_flow_nc,'color','r','linewidth',2)
hold on
grid on
plot(x,mass_flow_c,'color','b','linewidth',2)
hold on
xlabel('length of the nozzle , x')
ylabel('mass flow rate of non conservative')

title('comparing the mass flow rate for both consevative & non conservative form')

FUNCTION FRO NON CONSERVATIVE FORM :

function [mass_flow_nc, pressure_nc, Mach_nc, rho, v, t, rho_throat_nc, v_throat_nc, t_throat_nc, mass_flow_throat_nc, pressure_throat_nc, Mach_throat_nc] = non_conservative_form(x,dx,n,gamma,nt)

%calculate initial profile
rho=1-0.3146*x;
t=1-0.2314*x;
v=(0.1+1.09*x).*t.^0.5;
%area
a=1+2.2*(x-1.5).^2;
throat=find(a==1);
c=0.5;



%time loop
for k=1:nt
    
    dt=min(c.*dx./(sqrt(t)+v));
    rho_old=rho;
    v_old=v;
    t_old=t;
    
    %predictor method
    for j=2:n-1
        dv_dx=(v(j+1)-v(j))/dx;
        drho_dx=(rho(j+1)-rho(j))/dx;
        dt_dx=(t(j+1)-t(j))/dx;
        dloga_dx=(log(a(j+1))-log(a(j)))/dx;
        %continuity equation
        drho_dt_p(j)=-(rho(j)*dv_dx)-(rho(j)*v(j)*dloga_dx)-(v(j)*drho_dx);
        %momentum equation
        dv_dt_p(j)=-(v(j)*dv_dx)-((1/gamma)*(dt_dx+(t(j)/rho(j))*drho_dx));
        %energy equation
        dt_dT_p(j)=-(v(j)*dt_dx)-((gamma-1)*t(j)*(dv_dx+(v(j)*dloga_dx)));
        
        %solution
        rho(j)=rho(j)+drho_dt_p(j)*dt;
        v(j)=v(j)+dv_dt_p(j)*dt;
        t(j)=t(j)+dt_dT_p(j)*dt;
    end
    
    
    %corrector method
    for j=2:n-1
        dv_dx=(v(j)-v(j-1))/dx;
        drho_dx=(rho(j)-rho(j-1))/dx;
        dt_dx=(t(j)-t(j-1))/dx;
        dloga_dx=(log(a(j))-log(a(j-1)))/dx;
        %continuity equation
        drho_dt_c(j)=-(rho(j)*dv_dx)-(rho(j)*v(j)*dloga_dx)-(v(j)*drho_dx);
        %momentum equation
        dv_dt_c(j)=-(v(j)*dv_dx)-((1/gamma)*(dt_dx+(t(j)/rho(j))*drho_dx));
        %energy equation
        dt_dT_c(j)=-(v(j)*dt_dx)-((gamma-1)*t(j)*(dv_dx+(v(j)*dloga_dx)));
        
        
    end
    
    
    %compute the 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
    for i=2:n-1
        rho(i)=rho_old(i)+drho_dt(i)*dt;
        v(i)=v_old(i)+dv_dt(i)*dt;
        t(i)=t_old(i)+dt_dT(i)*dt;
    end
    
    %Apply the boundary conditions
    %Inlet
    v(1)=2*v(2)-v(3);
    %outlet
    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);
    
    %find the mass flow rate
    mass_flow_nc=rho.*a.*v;
    pressure_nc=rho.*t;
    Mach_nc=v./(sqrt(t));
    
    %variables at throat section
    v_throat_nc(k)=v(throat);
    rho_throat_nc(k)=rho(throat);
    t_throat_nc(k)=t(throat);
    mass_flow_throat_nc(k)=mass_flow_nc(throat);
    pressure_throat_nc(k)=pressure_nc(throat);
    Mach_throat_nc(k)=Mach_nc(throat);
    
    
    %plot for mass flow rate at different time step
    figure(1)
    if k==50
        plot(x,mass_flow_nc,'color','r','linewidth',2)
        hold on
    elseif k==100
        plot(x,mass_flow_nc,'color','b','linewidth',2)
        hold on
    elseif k==300
        plot(x,mass_flow_nc,'color','g','linewidth',2)
        hold on
    elseif k==500
        plot(x,mass_flow_nc,'color','r','linewidth',2)
        hold on
    elseif k==800
        plot(x,mass_flow_nc,'color','y','linewidth',2)
        hold on
    elseif k==1000
        plot(x,mass_flow_nc,'color','m','linewidth',2)
        hold on
    end
    
    
    
    
end

xlabel('length of the nozzle')
ylabel('mass flow rate at different time step')
title('mass flow rate at different time step ,Non conservative form')
legend('50th time step','100th time step','300th time step','500th time step','800th time step','1000th time step')
axis([0 3 0 2])

end

FUNCTION FOR CONSERVATIVE FORM :

function [mass_flow_c, pressure_c, Mach_c, rho_c, v_c, t_c, rho_throat_c, v_throat_c, t_throat_c, mass_flow_throat_c, pressure_throat_c, Mach_throat_c] = conservative_form(x,dx,n,gamma,nt)

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
a = 1+2.2*(x-1.5).^2;

throat=find(a==1);
c=0.4;


    % Solution Vectors 
 
v_c = 0.59./(rho_c.*a);
U1= rho_c.*a;
U2= rho_c.*a.*v_c;
U3= rho_c.*a.*((t_c./(gamma-1))+(gamma/2).*v_c.^2);

 
 
% Solution Vectors 
%% Time loop 
 for k = 1:(nt) 
 % Storing Initial Values 
    U1_old = U1;
    U2_old = U2;
    U3_old = U3;
    
    dt =min((c.*dx)./(t_c.^(0.5)+v_c));
   
    F1= U2;
    F2= ((U2.^2)./U1 )+ ((gamma-1)/gamma).*(U3-((gamma/2) * (U2.^2)./U1));
    F3= ((gamma*U2.*U3)./U1) - ((gamma*(gamma-1)/2) * (U2.^3)./(U1.^2));
      
    for j = 2:n-1 
    %% Predictors step 
      % Source term J2
      J2_c(j) = (1/gamma)*rho_c(j)*t_c(j)*((a(j+1)-a(j))/dx);  
      % Predictors Step derivative  
      dU1_dt_p(j) = -(F1(j+1)-F1(j))./dx ;
      dU2_dt_p(j) = -(F2(j+1)-F2(j))./dx + J2_c(j) ;
      dU3_dt_p(j) = -(F3(j+1)-F3(j))./dx ;
      
      % Solution update 
      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 
    rho_c =U1./a ;
    v_c=U2./U1;
    t_c = (gamma-1).*((U3./U1) - (gamma./2).*(v_c).^2);
    
    % correctors step 
    % Flux Vector 
    F1= U2;
    F2= ((U2.^2)./U1 )+ ((gamma-1)/gamma).*(U3-((gamma/2) * (U2.^2)./U1));
    F3= ((gamma*U2.*U3)./U1) - ((gamma*(gamma-1)/2) * (U2.^3)./(U1.^2));
    
    
    for j = 2:n-1 
    % Source term J2
    J2_c(j) = (1/gamma)*rho_c(j)*t_c(j)*((a(j)-a(j-1))/dx);
    % Correctors Step derivative  
    dU1_dt_c(j) = -(F1(j)-F1(j-1))./dx ;
    dU2_dt_c(j) = -(F2(j)-F2(j-1))./dx + J2_c(j) ;
    dU3_dt_c(j) = -(F3(j)-F3(j-1))./dx ;

    end

    %Average terms
    dU1_dt_avg = 0.5.*(dU1_dt_p + dU1_dt_c);
    dU2_dt_avg = 0.5.*(dU2_dt_p + dU2_dt_c);
    dU3_dt_avg = 0.5.*(dU3_dt_p + dU3_dt_c);

    %% Final Solution update 
    for i =2:n-1 
      
        U1(i) = U1_old(i) + dU1_dt_avg(i).*dt;
        U2(i) = U2_old(i) + dU2_dt_avg(i).*dt;
        U3(i) = U3_old(i) + dU3_dt_avg(i).*dt;  
     
    end 

%% Apply Boundary Conditions 
% Inlet Boundary Conditions
    %boundary conditions

    U1(1) = rho_c(1).*a(1);
    U2(1)= 2*U2(2)-U2(3);
    U3(1)= U1(1)*((t_c(1)/(gamma-1))+(gamma/2)*v_c(1).^2);
    
    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);

    % calcultaing primitive variables
    rho_c = U1./a;
    v_c= U2./U1;
    t_c = (gamma-1)*((U3./U1)-(gamma/2)*v_c.^2);
   
    pressure_c=rho_c.*t_c;
    Mach_c=v_c./(sqrt(t_c));
    mass_flow_c=rho_c.*a.*v_c;

     %variables at throat section
    v_throat_c(k)=v_c(throat);
    rho_throat_c(k)=rho_c(throat);
    t_throat_c(k)=t_c(throat);
    mass_flow_throat_c(k)=mass_flow_c(throat);
    pressure_throat_c(k)=pressure_c(throat);
    Mach_throat_c(k)=Mach_c(throat);
    
     %plot for mass flow rate at different time step
   figure(4)
    if k==50
        plot(x,mass_flow_c,'color','r','linewidth',2)
        hold on
    elseif k==100
        plot(x,mass_flow_c,'color','b','linewidth',2)
        hold on
    elseif k==300
        plot(x,mass_flow_c,'color','g','linewidth',2)
        hold on
    elseif k==500
        plot(x,mass_flow_c,'color','r','linewidth',2)
        hold on
    elseif k==800
        plot(x,mass_flow_c,'color','y','linewidth',2)
        hold on
    elseif k==1000
        plot(x,mass_flow_c,'color','m','linewidth',2)
        hold on
    end
xlabel('length of the nozzle')
ylabel('mass flow rate at different time step')
title('mass flow rate at different time step , conservative form')
legend('50th time step','100th time step','300th time step','500th time step','800th time step','1000th time step')
axis([0 3 0 2])

end        
    

OUTPUTS :

PLOTS :

i) Mass flow rate at different time step ,Non conservative form :

* The above figure shows that the variation of mass flow rate with respect to time step for non-conservation form. The graph shows that the mass flow rate is almost stable after the 500th time step.

ii) Steady-state distribution of primitive variables inside the nozzle- non conservative form :

* The above graph shows that the variation of flow field varialbles along the length of the nozzle for the non conservative form. In the above graph it is clear that the flow field variables such as density, pressure, temperature decreases and mach number increases with a increase in length of the nozzle.

iii) Time-wise variation of the primitive variables- non conservative form :

* The variation of flow field variables at throat section with respect to time step for non-conservative form as shown above. The graph shows that the flow field variables are almost reached the steady state solution after the 400th iteration.

iv) Mass flow rate at different time step ,conservative form :

* The graph shows the mass flow rate becomes almost stable after the 500th time step.

v) Steady-state distribution of primitive variables inside the nozzle- conservative form :

* The above graph shows that the variation of flow field varialbles along the length of the nozzle for the conservative form. The above graph shows that the flow field variables such as density, pressure, temperature decreases and mach number increases with a increase in length of the nozzle.

vi) Time-wise variation of the primitive variables- conservative form :

* The variation of flow field variables at throat section with respect to time step for conservative form as shown above. The graph shows that the flow field variables are almost reached the steady state solution after the 600th iteration.

vii) Comparison of Normalized mass flow rate distributions of both conservative and non conservative forms:

* The variation of mass flow rate at the both the forms is shown in the above figure. The graph shows that the mass flow rate of conservative forms produces a stable result when compare to that of non conservative form. It can also find that the mass flow rate is minimum at the throat then it expands in the non conservative form.

Grid dependence test :

To conduct the gride dependence test, we have to consider the number of gride as 31 and 61 , and just compare the results for flow variables at throat section at the end of 1400th time step.

The table shows that the solution of 61 grid points are almost equal when compare to the solution of 31 grid points. For more number of grid points will results in more computational time but it can able to produce the accurate result than the lesser grid point.