NUMERICAL SOLUTIONS OF QUASI-ONE-DIMENSIONAL NOZZLE FLOWS USING MACCORMACK METHOD

Aim: To simulate the isentropic flow through a quasi 1D subsonic-supersonic nozzle using MacCormack's technique.

Abstract: The flow is steady, isentropic, through a convergent-divergent nozzle as shown in the figure below. The flow at the inlet of the nozzle comes from the reservoir where the pressure and temperature are denoted by `p_0` and `T_0` respectively. The cross-sectional of the reservoir is very large (theoretically `Atobowtie`) and hence the velocity is very small (`V==0`). The flow expands isentropically at the supersonic speed at the exit, where the exit pressure, temperature, velocity, and Mach number are denoted by `p_e`, `T_e`, `V_e`, `M_e`. The flow is subsonic (`M<1`) at the inlet of the nozzle, sonic at the throat section (`M=1`).

Please note that in the energy equation e has been replaced with T

Non-conservative equation

The integral form of the continuity equation is represented below

`frac{del}{delt}intint_Vintrho**dV+intint_Srho**v**dS=0`

The equation is applied to the shaded control volume

The left side of the figure consist of `p`, `v`,`rho`,`T` and the right side of the variable is represented as `p+dp`, `v+dv`, `rho+drho`,`T+dT`. They are uniform over the area `A+dA`.

Continuity equation

The volume integral is represented as

`frac{del}{delt}intint_Vintrho**dV=frac{del}{delt}(rho**A**dx)`

The surface integral is represented as

`intint_Srho**v**dS=-rho**v**A+(rho+drho)(v+dv)(A+dA)`

The `-` sign in the equation is because dS always point outside the control volume.

Thus expanding the above equation

`intint_Srho**v**dS=-rho**v**A+rho**v**A+rho**v**dA`

`+rho**A**dv+rho**dv**dA+A**v**drho+v**dA**drho+A**dv**drho+drho**dv**dA`

The term consisting of product of differential approaches to zero. Thus the equation becomes

`intint_Srho**v**dS=rho**v**dA+rho**A**dv+A**v**drho=d(rho**A**v)`

Thus the equation becomes

`frac{del}{delt}(rho**A**dx)+d(rho**A**v)=0`

Thus dividing by dx throughout

`frac{del}{delt}(rho**A)+frac{del}{delx}(rho**A**v)=0`

This equation represent the continuity equation for unsteady,quasi one-dimensional flow.

Momentum equation

Writing the integral form of the equation as

`frac{del}{delt}intint_Vint(rho**u)**dV+intint_s(rho**u**V)**dS=-intint_s(p**dS)_x`

`frac{del}{delt}intint_Vint(rho**u)**dV=frac{del}{delt}(rho**v**A**dx)`

`intint_s(rho**u**V)**dS=-rho**v^2**A+(rho+drho)(v+dv)^2(A+dA)`

`intint_s(p**dS)_x=-p**A+(p+dp)**(A+dA)-2**p**(frac{dA}{2})`

Solving for each term and cancelling the product of differential terms yield

`frac{del}{delt}(rho**v**A**dx)+d(rho**v^2**A)=-Adp`

Dividing by dx and taking limit `dxto0` yields partial differential equation of the form

`frac{del}{delt}(rho**v**A)+frac{del}{delx}(rho**v^2**A)=-A**frac{delp}{delx}`

This equation represent conservation form of the momentum equation.

Taking the `v` term out of the differential equation gives

`frac{del}{delt}(rho**v**A)+frac{del}{delx}(rho**v^2**A)=-A**frac{delp}{delx}`

Multiplying continuity equation by v

`v**frac{del}{delt}(rho**A)+v**frac{del}{delx}(rho**A**v)=0`

Substracting the equation

`frac{del}{delt}(rho**v**A)-v**frac{del}{delt}(rho**A)+frac{del}{delx}(rho**v^2**A)-v**frac{del}{delx}(rho**A**v)=-A**frac{delp}{delx}`

Expanding the derivatives and canceling the like terms

`rho**A**frac{delv}{delt}+rho**A**v**frac{delv}{delx}=-A**frac{delp}{delx}`

Dividing by A

`rho**frac{delv}{delt}+rho**v**frac{delv}{delx}=-frac{delp}{delx}`

The equation represent momentum equation for quasi-1d flow in non conservative form.

Energy equation

For an adiabatic flow (`q=0`) with no body forces and no viscous effects, the integral form of the equation is represented as

`frac{del}{delt}intint_Vint(e+frac{v^2}{2})dV+intint_S(rho+frac{v^2}{2})**V**dS=-intint_s(rho**V)dS`

Expanding the above equation

`frac{del}{delt}[rho**(e+frac{v^2}{2})**A**dx]-rho**(e+frac{v^2}{2})**V**A+(rho+drho)**[e+de+frac{(v+dv)^2}{2}]**(v+dv)**(A+dA)`

`=-[-rho**v**A+(p+dp)**(v+dv)**(A+dA)-2**(frac{p**v**dA}{2})]`

Neglecting the product of differentiation and cancelling the like terms we have

`frac{del}{delt}[rho(e+frac{v^2}{2})**A**dx]+d(rho**e**v**A)+frac{d**(rho**v^3**A)}{2}=-d(rho**A**v)`

or

`frac{del}{delt}[rho(e+frac{v^2}{2})**A**dx]+d[rho(e+frac{v^2}{2})**v**A]=-d(rho**A**v)`

Dividing by dx throughout and taking limit `dxto0` the partial differential equation is obtained

`frac{del(rho**(e+frac{v^2}{2})**A)}{delt}+frac{del(rho**(e+frac{v^2}{2}))}{delx}*=-frac{del(rho**A**v)}{delx}`

The above equation is the conservation form of the energy equation.

In order to achieve non conservative form of energy equation

Multiplying v in conservation form of the momentum equation

`frac{del}{delt}(rho**frac{v^2}{2}**A)+frac{del}{delx}(rho**frac{v^3}{3}**A)=-A**v**frac{delp}{delx}`

Substracting from conservation form of energy equation

`frac{del(rho**e**A)}{delt}+frac{del(rho**e**v**A)}{delx}=-p**frac{del(A**v)}{delx}`

Multiplying the above equation by e gives

`e**frac{del(rho**A)}{delt}+e**frac{del(rho**A**v)}{delx}=0`

`rho**A**frac{dele}{delt}+(rho**A**v)**frac{dele}{delx}=-pfrac{del(A**v)}{delx}`

Dividing by A throughout

`rho**frac{dele}{delt}+rho**v**frac{dele}{delx}=-rho**frac{delv}{delx}-rho**v**frac{del(lnA)}{delx}`

The above equation is non-conservative form of the energy equation.

Non dimensional form of the equation

Please note that while using the Maccormack method we will be using the non-dimensional form of the above equation for ease in the calculation.

Non-dimensional temperature `T'=frac{T}{T_0}` `T_0` is the reservoir temperature

Non-dimensional density `rho'=frac{rho}{rho_0}` `rho_0` is the reservoir density.

Non-dimensional length `x'=frac{x}{L}` `L` is the length of nozzle.

Non-dimensional velocity `v'=frac{v}{a_0}` `a_0` is the speed of the sound in the reservoir.

Non-dimensional time `t'=frac{t}{frac{L}{a_0}}`

Non-dimensional area `A'=frac{A}{A^**}` `A^**` is the sonic throat area.

Replacing the above equation with non-dimensional variables we get

Continuity equation = `frac{delrho^'}{delt^'}=-rho^'**frac{delv^'}{delx^'}-rho^'**v^'**frac{del(lnA^')}{delx^'}-v^'**frac{delrho^'}{delx^'}`

Momentum equation= `frac{delv^'}{delt^'}=-v^'**frac{delv^'}{delx^'}-frac{1}{gamma}**(frac{delT^'}{delx^'}+frac{T^'}{rho^'}**frac{delrho^'}{delx^'})`

Energy equation `frac{delT^'}{delt^'}=-v^'**frac{delT^'}{delx^'}-(gamma-1)**T^'**[frac{delv^'}{delx^'}+v^'**frac{del(lnA^')}{delx^'}]`

Conservative equation

The conservative form of equation can be expressed in terms of solution vector `U` flux vector `F` and the source term `J` which are as follows

`U_1=rho^'**A^'`

`U_2=rho^'**A^'**v^'`

`U_3=rho^'(frac{T^'}{gamma-1}+frac{gamma}{2}**v^('2))**A^'`

`F_1=rho^'**A^'**v^'`

`F_2=frac{U_2^(2)}{U_1}+frac{gamma-1}{gamma}(frac{U_3}{U_1}-frac{gamma}{2}**(frac{U_2}{U_1})^2)`

`F_3=gamma**(frac{U_2**U_3}{U_1})-frac{gamma(gamma-1)}{2}**frac{U_2^3}{U_1^2}`

`J_2=frac{rho^'**T^'}{gamma}**frac{delA}{delx}`

MacCormack method

It is used for discretising hyperbolic partial differential equation. It involves predictor and corrector method.

Predictor method: It uses forward differencing for space derivative. In this the provisional values of variables (`rho` `T``v`) at (`n+1`) timesteps.

Corrector method: It uses backward differencing for space derivative. In this the predicted values are corrected. The timesteps are used `frac{trianglet}{2}`

The method is used for non linear equation (inviscid Burgers equation, Euler equation). For non linear equation it provides the best result. For linear equation Lax-Wendroff method is used.

The above concept have been applied to produce the Matlab code as shown below

• Matlab code for the main program.

% main program
clear
close all
clc

% defining one dimensional mesh size
L=3;% length of nozzle
n=31;% number of grid points
x=linspace(0,L,n);
dx=x(2)-x(1);


% defining other input parameters
gamma=1.4;
C=0.5;
nt=linspace(1,1400,1400);

% writing function for non conservative form
[mfr_analytical,mfr_e]=nozzle(x,L,n,gamma,nt);
[mfr_con]=conservative(x,L,n,gamma,nt);

figure(7);
plot(x,mfr_analytical)
hold on
plot(x,mfr_con,'k')
hold on
plot(x,mfr_e,'r')
xlabel('Length of nozzle')
ylabel('Non dimensional mass flow rate')
legend('Analytical mass flow rate','Mass flow rate for conservative form','Mass flow rate for non conservative form')
title('Comparison of Normalised mass flow rate distributions of both the forms')

• Creating user-defined function for the conservative equation.

function [mfr_con]=conservative(x,L,n,gamma,nt)

%% Conservative form
clear
close all
clc

% defining one dimensional mesh size
L=3;% length of nozzle
n=31;% number of grid points
x=linspace(0,L,n);
dx=x(2)-x(1);
a=1+2.2*(x-1.5).^2;% area
% throat
throat=find(a==1);
%CFL number
C=0.5;



% defining other input parameters
gamma=1.4;
% intialising temperature and density
for i=1:length(x)
if x(i)>=0 && x(i)<0.5
    rho_con(i)=1;
    temp_con(i)=1;
elseif x(i)>=0.5 && x(i)<1.5
    rho_con(i)=1-(0.366.*(x(i)-0.5));
    temp_con(i)=1-(0.167.*(x(i)-0.5));
elseif x(i)>=1.5 && x(i)<=3
    rho_con(i)=0.634-(0.3879.*(x(i)-1.5));
    temp_con(i)=0.833-(0.3507.*(x(i)-1.5));
end
end

% calculating primitive variables
v_con=0.59./(rho_con.*a);
p_con=rho_con.*temp_con;

% calculating solution vector
U1=rho_con.*a;
U2=rho_con.*a.*v_con;
U3=rho_con.*a.*((temp_con)/(gamma-1)+(0.5*gamma*v_con.^2));



nt=linspace(1,1400,1400);
tic;
iteration=0;
for k=1:length(nt)
    % predictor step
    % creating copy of solution vector
     U1_old=U1;
     U2_old=U2;
     U3_old=U3;
    % calculating dt
    dt_con=min((C*dx)./((sqrt(temp_con)+v_con))); 
    % calculating 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);
    % calculating source term
    for j=2:n-1
        J2(j)=(rho_con(j).*temp_con(j).*(a(j+1)-a(j)))/(gamma*dx);
        dU1dt_p(j)=-(F1(j+1)-F1(j))/dx;
        dU2dt_p(j)=-(F2(j+1)-F2(j))/(dx)+J2(j);
        dU3dt_p(j)=-(F3(j+1)-F3(j))/dx;
        % updating values 
        U1(j)=U1(j)+(dU1dt_p(j)*dt_con);
        U2(j)=U2(j)+(dU2dt_p(j)*dt_con);
        U3(j)=U3(j)+(dU3dt_p(j)*dt_con);
    end
    % updating primitive values
    rho_con=U1./a;
    v_con=U2./U1;
    temp_con=(gamma-1)*(U3./U1-(gamma/2)*v_con.^2);
    p_con=rho_con.*temp_con;
    % updating 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);
    
    %corrector step
    for j=2:n-1
        J2(j)=(rho_con(j).*temp_con(j).*(a(j)-a(j-1)))/(gamma*dx);
        dU1dt_c(j)=-(F1(j)-F1(j-1))/dx;
        dU2dt_c(j)=-(F2(j)-F2(j-1))/(dx)+J2(j);
        dU3dt_c(j)=-(F3(j)-F3(j-1))/dx;
    end
    % average time step value
    dU1dt_a=0.5*(dU1dt_p+dU1dt_c);
    dU2dt_a=0.5*(dU2dt_p+dU2dt_c);
    dU3dt_a=0.5*(dU3dt_p+dU3dt_c);
    for j=2:n-1
        U1(j)=U1_old(j)+(dU1dt_a(j)*dt_con);
        U2(j)=U2_old(j)+(dU2dt_a(j)*dt_con);
        U3(j)=U3_old(j)+(dU3dt_a(j)*dt_con);
    end
    
    iteration=iteration+1;
    
    % left boundary condition
    U1(1)=rho_con(1).*a(1);
    U2(1)=2*U2(2)-U2(3);
    U3(1)=U1(1)*(temp_con(1)./(gamma-1)+(gamma/2)*v_con(1).^2);
    % right boundary condition
    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);
    % calculating the primitive variables at exit
    rho_con=U1./a;
    v_con=U2./U1;
    temp_con=(gamma-1)*(U3./U1-(gamma/2)*v_con.^2);
    p_con=rho_con.*temp_con;
    mfr_con=rho_con.*v_con.*a;
    M_con=v_con./temp_con.^0.5;
    % calculating throat values
    rho_con_throat(k)=rho_con(throat);
    v_con_throat(k)=v_con(throat);
    temp_con_throat(k)=temp_con(throat);
    mfr_con_throat(k)=mfr_con(throat);
    M_con_throat(k)=M_con(throat);
    p_throat(k)=p_con(throat);
    mfr(iteration,:)=rho_con.*a.*v_con;
end
computational_time=toc;
figure(4);
    plot(x,mfr(1,:))
    hold on
    plot(x,mfr(50,:))
    hold on
    plot(x,mfr(100,:))
    hold on
    plot(x,mfr(150,:))
    hold on
    plot(x,mfr(200,:))
    hold on
    plot(x,mfr(700,:))
    hold on
    plot(x,mfr(1400,:))
    hold on
    xlabel('Length of nozzle')
    ylabel('Non dimensional mass flow rate')
    legend('iteration=1','iteration=50','iteration=100''iteration=150','iteration=200','iteration=700','iteration=1400')
    title('Variation of non dimensional mass flow rate in conservative form with time steps=1400.')
    figure(5);
    subplot(5,1,1)
    plot(nt,rho_con_throat)
    xlabel('Number of time steps')
    ylabel('Density')
    legend('Density')
    title('Timewise variation of primitive variables at throat for timesteps 1400 in conservative form')
    subplot(5,1,2)
    plot(nt,v_con_throat)
    xlabel('Number of time steps')
    ylabel('Velocity')
    legend('Velocity')
    subplot(5,1,3)
    plot(nt,temp_con_throat)
    xlabel('Number of time steps')
    ylabel('Temperature')
    legend('Temperature')
    subplot(5,1,4)
    plot(nt,M_con_throat)
    xlabel('Number of time steps')
    ylabel('Mach number')
    legend('Mach number')
    subplot(5,1,5)
    plot(nt,p_throat)
    xlabel('Number of time steps')
    ylabel('Pressure')
    legend('Pressure')
    figure(6);
    subplot(5,1,1)
    plot(x,rho_con)
    xlabel('Length of nozzle')
    ylabel('Density')
    legend('Density')
    title('Length wise variation of primitive variables for conservative form')
    subplot(5,1,2)
    plot(x,v_con)
    xlabel('Length of nozzle')
    ylabel('Velocity')
    legend('Velocity')
    subplot(5,1,3)
    plot(x,temp_con)
    xlabel('Length of nozzle')
    ylabel('Temperature')
    legend('Temperature')
    subplot(5,1,4)
    plot(x,p_con)
    xlabel('Length of nozzle')
    ylabel('Pressure')
    legend('Pressure')
    subplot(5,1,5)
    plot(x,M_con)
    xlabel('Length of nozzle')
    ylabel('Mach number')
    legend('Mach number')
end

• Creating user-defined function for the non-conservative equation.

function [mfr_analytical,mfr_e]=nozzle(x,L,n,gamma,nt)

%Non conservative form

clear
close all
clc

% defining one dimensional mesh size
L=3;% length of nozzle
n=31;% number of grid points
x=linspace(0,L,n);
dx=x(2)-x(1);


% defining other input parameters
gamma=1.4;
C=0.5;
nt=linspace(1,1400,1400);


% nozzle properties at initial condition

mfr_analytical=0.579*ones(1,n);
a=1+2.2*(x-1.5).^2;% area
rho=1-0.3146.*x;% density
t=1-0.2314.*x;% temperature
v=(0.1+1.09*x).*t.^0.5;% velocity
mfr_initial=rho.*a.*v; % mass flow rate initial
M_initial=v./t.^0.5; % mach number initial
p_initial=rho.*t; % pressure initial

% calculating throat value
throat=find(a==1);
% time interval value
dt=min((C*dx)./((sqrt(t)+v))); 
number_iter=0;

% starting time loop
tic;
for k=1:length(nt)
% creating a copy of profile for comparison
rho_old=rho;
v_old=v;
t_old=t;

    
    % space loop
    
    for j=2:n-1
        
        % defining terms for calculation convenience
        % predictor step
        dvdx_p=(v(j+1)-v(j))/dx;
        dlogadx_p=(log(a(j+1))-log(a(j)))/dx;
        drhodx_p=(rho(j+1)-rho(j))/dx;
        dtdx_p=(t(j+1)-t(j))/dx;
        
        % continuity equation
        drhodt_p(j)=(-rho(j).*dvdx_p)+(-rho(j).*v(j).*dlogadx_p)+(-v(j).*drhodx_p);
        
        % momentum equation
        dvdt_p(j)=(-v(j).*dvdx_p)+((-1/gamma)*(dtdx_p+(t(j)./rho(j)).*(drhodx_p)));
        
        % energy equation
        dtdt_p(j)=(-v(j).*dtdx_p)-(gamma-1).*t(j).*(dvdx_p+(v(j).*dlogadx_p));
        
        % updating v,rho and t
        v(j)=v(j)+dvdt_p(j).*dt;
        rho(j)=rho(j)+drhodt_p(j).*dt;
        t(j)=t(j)+dtdt_p(j).*dt;
    end
    
    for j=2:n-1
        
        % corrector step
        dvdx_c=(v(j)-v(j-1))/dx;
        dlogadx_c=(log(a(j))-log(a(j-1)))/dx;
        drhodx_c=(rho(j)-rho(j-1))/dx;
        dtdx_c=(t(j)-t(j-1))/dx;
        
        % continuity equation
        drhodt_c(j)=(-rho(j).*dvdx_c)+(-rho(j).*v(j).*dlogadx_c)+(-v(j).*drhodx_c);
        
        % momentum equation
        dvdt_c(j)=(-v(j).*dvdx_c)+((-1/gamma)*(dtdx_c+(t(j)./rho(j)).*(drhodx_c)));
        
        % energy equation
        dtdt_c(j)=(-v(j).*dtdx_c)-(gamma-1).*t(j).*(dvdx_c+(v(j).*dlogadx_c));
        
        
    end
    
    % average time step
    dvdt_a=0.5*(dvdt_c+dvdt_p);
    dtdt_a=0.5*(dtdt_c+dtdt_p);
    drhodt_a=0.5*(drhodt_c+drhodt_p);
    
    % updating final values
    for j=2:n-1
       v(j)=v_old(j)+dvdt_a(j).*dt;
       t(j)=t_old(j)+dtdt_a(j).*dt;
       rho(j)=rho_old(j)+drhodt_a(j).*dt;
    end
   
    number_iter=number_iter+1;
    
    % applying boundary condition
    v(1)=2*v(2)-v(3);
    rho(1)=1;
    t(1)=1;
    v(n)=2*v(n-1)-v(n-2);
    t(n)=2*t(n-1)-t(n-2);
    rho(n)=2*rho(n-1)-rho(n-2);
    
    % primitive variables at exit
    p_e=rho.*t;
    M_e=v./t.^0.5;
    mfr_e=rho.*a.*v;
    
    % value at throat values
    mfr_throat(k)=mfr_e(throat);
    p_throat(k)=p_e(throat);
    M_throat(k)=M_e(throat);
    v_throat(k)=v(throat);
    rho_throat(k)=rho(throat);
    t_throat(k)=t(throat);
    figure(1);
    if (number_iter==1)
        plot(x,mfr_e)
        hold on
    elseif (number_iter==50)
        plot(x,mfr_e)
        hold on
    elseif (number_iter==100)
        plot(x,mfr_e)
        hold on
    elseif (number_iter==150)
        plot(x,mfr_e)
        hold on
    elseif (number_iter==200)
        plot(x,mfr_e)
        hold on
    elseif (number_iter==700)
        plot(x,mfr_e)
        hold on
    elseif (number_iter==1400)
        plot(x,mfr_e)
        hold on
        xlabel('Length of nozzle')
        ylabel('Non dimensional mass flow rate')
        legend('number iter==1','number iter==50','number iter==100','number iter==200','number iter==700','number iter==1400','mfr analytical','mfr initial')
        title('Variation of non dimensional mass flow rate along the length of nozzle with varying time steps in NC form.')
    end
end
     computational_time=toc;

figure(2);
subplot(5,1,1)
plot(nt,p_throat)
xlabel('Number of time steps')
ylabel('Pressure')
legend('Pressure at throat')
title('Variation of primitive variables at throat for time steps 1400 in NC form')
subplot(5,1,2)
plot(nt,rho_throat)
xlabel('Number of time steps')
ylabel('Density')
legend('Density at throat')
subplot(5,1,3)
plot(nt,t_throat)
xlabel('Number of time steps')
ylabel('Temperature')
legend('Temperature at throat')
subplot(5,1,4)
plot(nt,v_throat)
xlabel('Number of time steps')
ylabel('Velocity')
legend('Velocity at throat')
subplot(5,1,5)
plot(nt,M_throat)
xlabel('Number of time steps')
ylabel('Mach number')
legend('Mach number at throat')
figure(3);
subplot(5,1,1)
plot(x,p_e)
xlabel('Length of nozzle')
ylabel('Pressure')
legend('Pressure at exit')
title('Variation of primitive variables at exit along the length of nozzle')
subplot(5,1,2)
plot(x,rho)
xlabel('Length of nozzle')
ylabel('Density')
legend('Density at exit')
subplot(5,1,3)
plot(x,t)
xlabel('Length of nozzle')
ylabel('Temperature')
legend('Temperature at exit')
subplot(5,1,4)
plot(x,v)
xlabel('Length of nozzle')
ylabel('Velocity')
legend('Velocity at exit')
subplot(5,1,5)
plot(x,M_e)
xlabel('Length of nozzle')
ylabel('Mach number')
legend('Mach number at exit')
end

Demonstration of grid independence.

S No	Number of grid points.	Condition at the nozzle throat
		Conservative form				Non Conservative form
		Density `rho`	Temperature `T`	Pressure `p`	Mach number `M`	Density `rho`	Temperature `T`	Pressure `p`	Mach number `M`
1	31	`0.6504`	`0.8400`	`0.5464`	`0.9817`	`0.6387`	`0.8365`	`0.5342`	`0.9994`
2	93	`0.6373`	`0.8359`	`0.5327`	`1.0001`	`0.6399`	`0.8365`	`0.5353`	`1.0035`
3	155	`0.6415`	`0.8372`	`0.5371`	`1.0034`	`0.6666`	`0.8507`	`0.5671`	`1.0197`
Exact analytical solution		`0.634`	`0.833`	`0.528`	`1.000`	`0.634`	`0.833`	`0.528`	`1.000`

The above table represents the steady state solution at the throat section of nozzle for both the cases. The same is true for all the location within the nozzle. The deviation from the exact solution is very marginal for grid number 93 and 155 respectively. Thus increasing the grid points will increase the numerical accuracy of the solution, but it will also increase the execution time which is to be kept in mind.

Thus we can say that instead of increasing the grid points we can execute the simulation for 31 grid points which will save our time.

Courant number effect

S No	Courant Number	Condition at the nozzle throat
		Conservative form				Non Conservative form
		Density `rho`	Temperature `T`	Pressure `p`	Mach number `M`	Density `rho`	Temperature `T`	Pressure `p`	Mach number `M`
1	`0.5`	`0.6504`	`0.8400`	`0.5464`	`0.9817`	`0.6387`	`0.8365`	`0.5342`	`0.9994`
2	`0.6`	`0.6478`	`0.8387`	`0.5433`	`0.9867`	`0.6389`	`0.8365`	`0.5344`	`0.9993`
3	`0.7`	`0.6457`	`0.8376`	`0.5408`	`0.9906`	`0.6391`	`0.8366`	`0.5346`	`0.9992`
4	`0.8`	`0.6441`	`0.8368`	`0.5389`	`0.9936`	`0.6393`	`0.8367`	`0.5349`	`0.9991`
5	`0.9`	`0.6430`	`0.8362`	`0.5376`	`0.9957`	`0.6395`	`0.8368`	`0.5351`	`0.9989`
6	`1.2`	Solution did not converge/Program is unstable and it blew up.
Exact analytical sloution		`0.634`	`0.833`	`0.528`	`1.000`	`0.634`	`0.833`	`0.528`	`1.000`

From the above table it can be concluded that when courant number is increased to 1.2 instabilities do occur and the program blows up. The CFL number holds good for flow involving linear equations, but for non-linear parabolic partial differential equation it can be useful for providing the values of `trianglet`.

Expected graphs

• Steady-state distribution of primitive variables inside the nozzle for the conservative equation.

• Steady-state distribution of primitive variables inside the nozzle for the non-conservative equation.

• Time-wise variation of the primitive variables for the conservative equation.

• Time-wise variation of the primitive variables for the non-conservative equation.

• Variation of Mass flow rate distribution inside the nozzle at different time steps during the time-marching process for the conservative equation.

• Variation of Mass flow rate distribution inside the nozzle at different time steps during the time-marching process for the non-conservative equation.

• Comparison of Normalized mass flow rate distributions of both forms

Conclusion

• The program was written for both the forms and was executed in matlab.
• It can be inferred from the "steady state distribution of primitive variables inside the nozzle for both the forms'' that there is no much difference of the values in both the form as it can be observed in grid independence test.
• It can be inferred from the "time wise variation of the primitive variables for both the forms" that the largest changes takes place at the earliest time-steps, after 400 timesteps the value of the primitive variables reaches to the steady state.
• It can be inferred from the "variation of the non dimensional mass flow rate along the nozzle for different time steps in conservative form" that the result after 100 iteration shows a humped distribution. The mass flow rate has started to converge after 200 iteration and after 700 iteration it has started to converge and becomes equal to 0.579.
• It can be inferred from the "variation of the non dimensional mass flow rate along the nozzle for different time steps in non conservative form" that the first iteration value is represented by the blue line graph. After 50 iteration the mass flow rate has changed drastically, and after 100 iteration the mass flow rate has changed radically. After 700 iteration the mass flow rate begins to settle down, and in between 700 to 1400 iteration the mass flow rate has converged and becomes steady.
• It can be inferred from the " comparison of normalised mass flow rate distribution of both the forms" that the mass flow rate predicted by the conservation form of the equation is much more satisfactory than the one obtained from the nonconservation form of the equation on the basis of two accounts: that the conservation form is much closer to 0.579, on the other hand in non-conservative equations some spurious oscillations can be found at the inlet and outlet of the nozzle. The conservation form of the equation is closer to analytical mass flow rate.
• The execution time for the conservation form of equation is `1.3661` and for the non-conservative form is `2.4725`.

Source: CFD book by John Anderson, Google