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

Problem definition:

• Simulate quasi-a1-d supersonic nozzle flow using maccormack method.
• To write separate functions for conservative and non conservative form.
• minimum number of cycles for which the simulation should be run in order for convergence.
• Compare the normalized mass flow rate between the conservative forms and non-conservative forms - make suitable conclusions
• Which method is faster?

schematic diagram for subsonic-supersonic isentropic nozzle flow

Introduction:

Flow through the nozzle whose area is varying wrt x is in reality atleast 2D but assuming flow properties vary only along x direction we solve the above problem, therefore it is called quasi-1d-flow.

As there is change in physics of problem by making above assumption we need to derive the governing equations for quasi-1d-flow from the integral form of control-volume for continuity, momentum and energy equations so that the resulting equations conserves mass, momentum and energy.

Integral form of continuty equation:

`partial/(partialt)intintintrhodvee + intintrhoV.ds = 0`

Integral form of momentum equation along x direction:

`partial/(partialt)intintint(rhou)dvee + intint(rhouV).ds = -intint(Pds)_x`

Integral form of energy equation along x direction:

`partial/(partialt)intintint_vrho(e+V^2/2)dvee + intint_srho(e+v^2/2)V.ds = -intint_s(PV).ds`

by solving the above equations for a nozzle in the above figure by considering a small control volume of length 'dx' we can obtain governing equations for quasi-1d-flow:

Governing equations in conservative form:

1. Continuity equation: `(partial(rhoA))/(partialt) + (partial(rhoAV))/(partialx) = 0`

2. Momentum equation: `(partial(rhoAV))/(partialt) + (partial(rhoAV^2))/(partialx) = -A(partialP)/(partialx)`

3. Energy equation: `(partial[rho(e+V^2/2)A])/(partialt) + (partial[rho(e+V^2/2)AV])/(partialx) = -(partial(PAV))/(partialx)`

Governing equations in non-conservative form:

1. Continuity equation: `(partial(rhoA))/(partialt) + rhoA(partialV)/(partialx) + partialV(partialA)/(partialx) + VA(partialrho)/(partialx) = 0`

2. Momentum equation: `rho(partialV)/(partialt) + rhoV(partialV)/(partialx) = -R(rho(partialT)/(partialx) + T(partialrho)/(partialx))`

3. Energy equation: `rhoc_v(partialT)/(partialt) +rhoVc_v(partialT)/(partialx) = -rhoRT[(partialV)/(partialx) + V(partial(lnA))/(partialx) ]`

Non-dimensionalization of primitive variables:

`T' = T/T_0`, `rho' = rho/rho_0` , `x' = x/l`, `V' = V/V_o`, `A' = A/A_0`, `P' = P/P_0`

this gives

Non-dimensional form of governing differential equations for conservative form:

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

2. Momentum equation: `(partial(rho'A'V'))/(partialt') + (partial(rho'A'v'^2 + 1/gamma(P'A')))/(partialx') = 1/gammaP'(partial(A'))/(partialx')`

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

Initial conditions:

`rho = 1 - 0.3146X`

`T = 1 - 0.2314X`

`V = (0.1+1.09X)T^(1/2)`

Area

`A = 1+2.2(X -1.5)^2` for `0<=X<=3`

Non-dimensional form of governing differential equations for non-conservative form:

1. Continuity equation: `(partialrho')/(partialt') = -rho'(partialV')/(partialx') - rho'V'(partial(lnA'))/(partialx') - V'(partialrho')/(partialx')`

2. Momentum equation: `(partialV')/(partialt') = -V'(partialV')/(partialx') - 1/gamma((partialT')/(partialx') + (T')/(rho')(partialrho')/(partialx'))`

3. Energy equation: `(partialT')/(partialt') = -V(partialT)/(partialx') - (gamma-1)T'[(partialV')/(partialx')+V'(partial(lnA'))/(partialx')]`

Initial conditions:

`rho' = 1, T' = 1` for `0<=X'<=0.5`

`rho' = 1 - 0.366(X'-0.5)` and `T' = 1 - 0.167(X' - 0.5)` for `0.5<=X'<=1.5`

`rho' = 0.634-0.3879(X' - 1.5)` and `T' = 0.833 - 0.3507(X' - 1.5)` for `1.5<=X'<=3.5`

Area

`A = 1+2.2(X -1.5)^2` for `0<=X<=3`

To obtain the solution the above quasi-1D-partial differential equation CFD techniques are used but why maccormack method?

• In previous project solving 2d heat conduction equation it is observed that explicit scheme is easy to handle, faster and accurate if CFl condition is valid.
• Here maccormack method is Explicit finite difference scheme with second order accuracy in both space and time.
• Maccormack method is easy to program inspite of crank-nicolson method which is also second order accurate in both space and time. Implicit and explicit method are only first order accurate wrt time and second order accurate wrt space.
• Quasi-1d-supersonic flow equation is hyperbolic type which is suitable to solve using maccormack method.

CFL condition:

As Maccormack method is an explicit scheme it is highly dependent on CFL condition

CFL number restricts the time step for a given number of grid point and viceversa so that solution remains stable i.e. ratio of advection distance travelled to the grid size(dx) is less than 1.

`(trianglet)_i^t = C(trianglex)/(a_i^t +V_i^t)`

where CFL number `C <=1` in this project C is taken as 0.5.

For a fixed C value time step varies for every grid point. so, minimum time step value among all of them is taken as a dt to find the variables for next time step and the same procedure is used to find dt for all time loops.

Procedure:

Main program:

clear all
close all
clc
%solving quasi-1d-nozzle flow using maccormack method both conservative and
%non conservative
gamma = 1.4;
n = 31;
throat = (n+1)/2;
x = linspace(0,3,n);
dx = 3/(length(x)-1);
A = 1 + 2.2*(x - 1.5).^2;
nt = 1400;
%function to solve non conservative form
[rho_nc,v_nc,t_nc,p_nc,m_nc,mach_nc,density_nc,temperature_nc,velocity_nc,pressure_nc,mach_number_nc,mass_flow_nc,rho_residual_nc,v_residual_nc,t_residual_nc,ct_nc] = non_c(gamma,x,dx,n,throat,A,nt);

%function to solve conservative form
[rho_c,v_c,t_c,p_c,m_c,mach_c,density_c,temperature_c,velocity_c,pressure_c,mach_number_c,mass_flow_c,u1_residual_c,u2_residual_c,u3_residual_c,ct_c] = conserv_form(gamma,x,dx,n,throat,A,nt);
time = 1:1:nt;


% plotting parameters wrt time to observe convergence
figure(3)
subplot(3,2,1)
plot(time,density_nc(time),'r',time,density_c(time),'g');
legend('non-conservative form','conservative form')
xlabel('time step')
ylabel('density')
title('density comparison wrt time')

subplot(3,2,2)
plot(time,temperature_nc(time),'r',time,temperature_c(time),'g');
legend('non-conservative form','conservative form')
xlabel('time step')
ylabel('temperature')
title('temperature comparison wrt time')

subplot(3,2,3)
plot(time,velocity_nc(time),'r',time,velocity_c(time),'g');
legend('non-conservative form','conservative form')
xlabel('time step')
ylabel('velocity')
title('velocity vs time')

subplot(3,2,4)
plot(time,pressure_nc(time),'r',time,pressure_c(time),'g');
legend('non-conservative form','conservative form')
xlabel('time step')
ylabel('pressure')
title('pressure comparison wrt time')

subplot(3,2,5)
plot(time,mach_number_nc(time),'r',time,mach_number_c(time),'g');
legend('non-conservative form','conservative form')
xlabel('time step')
ylabel('mach number')
title('mach number comparison wrt time')

subplot(3,2,6)
plot(time,mass_flow_nc(time),'r',time,mass_flow_c(time),'g');
legend('non-conservative form','conservative form')
xlabel('time step')
ylabel('mass flow rate')
title('mass flow rate comparison wrt time')

exact = 0.579*ones(1,n);
%plotting mass flow rate wrt x
figure(4)
plot(x,m_nc,'r',x,m_c,'g',x,exact,'b');
legend('non-conservative form','conservative form','exact analytical result')
xlabel('non dimensional distance along x axis')
ylabel('mass flow rate')
title('mass flow rate comparison wrt non-dimensional length of nozzle ')

%plotting all the observed parameters wrt x
figure(5)
subplot(3,2,1)
plot(x,rho_nc,'r',x,rho_c,'g');
legend('non-conservative form','conservative form')
xlabel('non dimensional length')
ylabel('density')
title('density comparison wrt length')

subplot(3,2,2)
plot(x,t_nc,'r',x,t_c,'g');
legend('non-conservative form','conservative form')
xlabel('Non dimensional length')
ylabel('temperature')
title('temperature comparison wrt length')

subplot(3,2,3)
plot(x,v_nc,'r',x,v_c,'g');
legend('non-conservative form','conservative form')
xlabel('Non dimensional length')
ylabel('velocity')
title('velocity vs length')

subplot(3,2,4)
plot(x,p_nc,'r',x,p_c,'g');
legend('non-conservative form','conservative form')
xlabel('Non dimensional length')
ylabel('pressure')
title('pressure comparison wrt length')

subplot(3,2,5)
plot(x,mach_nc,'r',x,mach_c,'g');
legend('non-conservative form','conservative form')
xlabel('Non dimensional length')
ylabel('mach number')
title('mach number comparison wrt length')

subplot(3,2,6)
plot(x,m_nc,'r',x,m_c,'g',x,exact,'b');
legend('non-conservative form','conservative form','exact')
xlabel('Non dimensional length')
ylabel('mass flow rate')
title('mass flow rate comparison wrt length')

ti = 1:10:1400;
figure(6)
subplot(2,1,1)
%plotting residuals for non-conservative form
loglog(ti,rho_residual_nc(ti),'--r',ti,v_residual_nc(ti),'--b',ti,t_residual_nc(ti),'--g')
xlabel('time step')
ylabel('residuals')
legend('drho/dt','dv/dt','dT/dt')
title('comparsion of residuals for Non-Conservative form')

subplot(2,1,2)
%plotting residuals for conservative form
loglog(ti,u1_residual_c(ti),'--r',ti,u2_residual_c(ti),'--b',ti,u3_residual_c(ti),'--g')
xlabel('time step')
ylabel('residuals')
legend('du1/dt','du2/dt','du3/dt')
title('comparsion of residuals for Conservative form')

figure(7)
%comparison of computational time
y = [ct_c ct_nc];
bar(y)
set(gca, 'XTick', 1:2, 'XTickLabels', {'conservative','non-conservative'})
title('comparison of computational time')

Explaination:

Few denotions used in above codes varying wrt non-dimensional length

	Density	velocity	temperature	pressure	mass flow rate	mach number
Non-conservative form	rho_nc	v_nc	t_nc	p_nc	m_nc	mach_nc
conservative form	rho_c	v_c	t_c	p_c	m_c	mach_c

Few other denotions varying wrt time

	Density	velocity	temperature	pressure	mass flow rate	mach number
Non-conservative form	density_nc	velocity_nc	temperature_nc	pressure_nc	mass_flow_nc	mach_number_nc
conservative form	density_c	velocity_c	temperature_c	pressure_c	mass_flow_c	mach_number_c

Also, ct_nc = computational time for non-conservative form including figure(1)

ct_c = computational time for conservative form including figure(2)

In the above code consevative form and non-conservative form are solved in the form of functions to make code user friendly.

Function for Non-conservative form:

function [rho,v,t,p,m,mach,density,temperature,velocity,pressure,mach_number,mass_flow,rho_residual,v_residual,t_residual,ct] = non_c(gamma,x,dx,n,throat,A,nt)
tic;
%initial conditions
rho = 1 - 0.3146*x;
t = 1 - 0.2314*x; %t is temperature
v = (0.1+1.09*x).*t.^0.5;

for k = 1:nt
    
    rho_old = rho;
    v_old = v;
    t_old = t;


%Maccormack method
for i = 2:n-1
    
    %technique
    dvdx = (v(i+1)-v(i))/dx;
    drhodx = (rho(i+1)-rho(i))/dx;
    dlogadx = (log(A(i+1)) - log(A(i)))/dx;
    dtdx = (t(i+1) - t(i))/dx;
    
    %predictor
    %continuity equation
    drhodt_p(i) = -rho(i)*dvdx - rho(i)*v(i)*dlogadx - v(i)*drhodx;
    %momentum equation
    dvdt_p(i) = -v(i)*dvdx - (1/gamma)*(dtdx + (t(i)/rho(i))*drhodx);
    %energy equation
    dtdt_p(i) = -v(i)*dtdx -(gamma-1)*t(i)*(dvdx + v(i)*dlogadx);
    
    %finding dt
    delt = (0.5*(dx./(t.^.5 + v)));
    delt(1) = (0.5*(dx/(t(1)^.5 + v(1))));
    dt = min(delt(2:n-1));
   
    %updating 
    rho(i) = rho(i) + drhodt_p(i)*dt;
    v(i) = v(i) + dvdt_p(i)*dt;
    t(i) = t(i) + dtdt_p(i)*dt;
end

    
for i = 2 : n-1
     %corrector
    dvdx = (v(i)-v(i-1))/dx;
    drhodx = (rho(i)-rho(i-1))/dx;
    dlogadx = (log(A(i)) - log(A(i-1)))/dx;
    dtdx = (t(i) - t(i-1))/dx;
    
    %continuity equation
    drhodt_c(i) = -rho(i)*dvdx - rho(i)*v(i)*dlogadx - v(i)*drhodx;
    %momentum equation
    dvdt_c(i) = -v(i)*dvdx - (1/gamma)*(dtdx + (t(i)/rho(i))*drhodx);
    %energy equation
    dtdt_c(i) = -v(i)*dtdx -(gamma-1)*t(i)*(dvdx + v(i)*dlogadx);
    
end

% compute average time derivatives
drhodt = 0.5*(drhodt_p + drhodt_c);
dvdt = 0.5*(dvdt_p + dvdt_c);
dtdt = 0.5*(dtdt_p + dtdt_c);

%finding new variables
for i = 2:n-1
    
    rho(i) = rho_old(i) + drhodt(i)*dt;
    v(i) = v_old(i) + dvdt(i)*dt;
    t(i) = t_old(i) + dtdt(i)*dt;
end

%applying boundary conditions
%inlet
v(1) = 2*v(2) - v(3);

%outlet
rho(n) = 2*rho(n-1) - rho(n-2);
v(n) = 2*v(n-1) - v(n-2);
t(n) = 2*t(n-1) - t(n-2);

p = rho.*t;
mach = v./sqrt(t);
m = rho.*A.*v;

%For storing primitive variables wrt time
density(k) = rho(throat);
temperature(k) = t(throat);
velocity(k) = v(throat);
pressure(k) = p(throat);
mach_number(k) = mach(throat);
mass_flow(k) = m(throat);

%for storing residuals wrt time
rho_residual(k) = abs(drhodt(throat));
v_residual(k) = abs(dvdt(throat));
t_residual(k) = abs(dtdt(throat));

z = 0*x;

%for plotting mass flow rate wrt x for different time steps
figure(1)
if k == 1
    plot(x,m,'--b')
    xlabel('Non dimensional length')
    ylabel('Mass flow rate(Non-conservative form)')
    title('Comparison of mass flow rate wrt length for various time iterations')
    hold on
end
    if k == 100
        plot(x,m,'--g')
        hold on
    end
    if k == 200
        plot(x,m,'--c')
        hold on
    end
    if k == 700
        plot(x,m,'r')
        hold on
    end
    if k == 1400
        plot(x,m,'*m')
        hold on
        legend('timestep = 1','timestep =100','timestep =200','timestep =700','timestep =1400')
    end

 ct = toc;
end

% m = rho.*A.*v;
% h(j) = plot(x,m);
% hold on

Explaination: In the above function non-conservative form of quasi-1d-equation is solved using maccormack method to obtain all the output parameters mentioned in the above table.

steps:

1. Initial values of density, velocity and temperature are to be provided

2. Predictor step

3. Finding dt using cfl condition

4. updating primitive variables

5. corrector step

6. finding average residuals

7. final update of primitive variables

8. Finding other dependent variables

9. boundary conditions

10. plotting variables wrt x for diferent time steps

Function for conservative form:

function [rho,v,t,p,m,mach,density,temperature,velocity,pressure,mach_number,mass_flow,u1_residual,u2_residual,u3_residual,ct] = conserv_form(gamma,x,dx,n,throat,A,nt)
tic;
%finding initial values of density(rho) and temperature(t)
[rho,t] = initialization(x,n,dx);
%finding v using U2 = 0.59
v = 0.59./(rho.*A);
%finding pressure
p = rho.*t;
%initial conditions for variables u1,u2 and u3
u1 = rho.*A;
u2 = rho.*A.*v;
u3 = rho.*A.*((t./(gamma-1)) + (gamma/2)*v.^2);

for k = 1:nt
    %storing old values
    u1_old = u1;
    u2_old = u2;
    u3_old = u3;

%finding dt
    delt = 0.5*(dx./(t.^.5 + v));
    delt(1) = (0.5*(dx/(t(1)^.5 + v(1))));
    dt = min(delt(2:n-1));
% predictor step
    f1 = u2;
    f2 = (u2.^2)./u1 + ((gamma -1)/gamma)*(u3 - 0.5*gamma*((u2.^2)./u1));
    f3 = gamma*u2.*u3./u1 - 0.5*gamma*(gamma-1)*((u2.^3)./u1.^2);

for j = 2:n-1
    %finding j2
    j2(j) = (1/gamma)*rho(j).*t(j)*((A(j+1)-A(j))/dx);
%finding derivatives using forward difference

    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 u1 u2 and u3 values
    u1(j) = u1(j) + du1dt_p(j)*dt;
    u2(j) = u2(j) + du2dt_p(j)*dt;
    u3(j) = u3(j) + du3dt_p(j)*dt;
end

%updating rho, v, t    
    rho = u1./A;
    v = u2./u1;
    t = (gamma - 1)*(u3./u1 - 0.5*gamma*(v.^2));
    p = rho.*t;
%updating f1, f2 and f3 for corrector step
    f1 = u2;
    f2 = ((u2.^2)./u1) + ((gamma -1)/gamma)*(u3 - 0.5*gamma*((u2.^2)./u1));
    f3 = (gamma*u2.*u3./u1) - (0.5*gamma*(gamma-1)*((u2.^3)./(u1.^2)));

for j = 2:n-1
    %finding j2
    j2(j) = (1/gamma)*rho(j).*t(j)*((A(j)-A(j-1))/dx);                
%corrector step
    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
%computing average time derivative
    du1dt = 0.5*(du1dt_p + du1dt_c);
    du2dt = 0.5*(du2dt_p + du2dt_c);
    du3dt = 0.5*(du3dt_p + du3dt_c);
    
%final corrected values of u1, u2 and u3
for i = 2:n-1
    u1(i) = u1_old(i) + du1dt(i)*dt;
    u2(i) = u2_old(i) + du2dt(i)*dt;
    u3(i) = u3_old(i) + du3dt(i)*dt;
end
    
%applying boundary conditions
% inlet
    u1(1) = A(1);
    u2(1) = 2*u2(2) - u2(3);
    u3(1) = u1(1)*((t(1)/(gamma-1)) + 0.5*gamma*(u2(1)/u1(1))^2);
    %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);
%updating primitive variables

    rho = u1./A;
    v = u2./u1;
    t = (gamma - 1)*((u3./u1) - 0.5*gamma*v.^2);
    p = rho.*t;
    mach = v./sqrt(t);
    m = rho.*A.*v;
    
  %plotting primitive variables wrt time
    density(k) = rho(throat);
    temperature(k) = t(throat);
    velocity(k) = v(throat);
    pressure(k) = p(throat);
    mach_number(k) = mach(throat);
    mass_flow(k) = m(throat);
    
   %for storing residuals wrt time
    u1_residual(k) = abs(du1dt(throat));
    u2_residual(k) = abs(du2dt(throat));
    u3_residual(k) = abs(du3dt(throat)); 
    
   %for plotting mass flow rate with varying time step
figure(2)
if k == 1
    plot(x,m,'--b')
    xlabel('Non dimensional length')
    ylabel('Mass flow rate (conservative form)')
    title('Comparison of mass flow rate wrt length for various time iterations ')
    hold on
end
    if k == 100
        plot(x,m,'--g')
        hold on
    end
    if k == 200
        plot(x,m,'--c')
        hold on
    end
    if k == 700
        plot(x,m,'r')
        hold on
    end
    if k == 1400
        plot(x,m,'*m')
        hold on
        legend('timestep = 1','timestep =100','timestep =200','timestep =700','timestep =1400')
    end
 
 ct = toc;   
end

Explaination:

In the above function conservative form of quasi-1d-equation is solved using maccormack method to obtain all the output parameters mentioned in the above table.

steps:

1. Initial values of density, velocity and temperature are to be provided using initialization function.

2. Finding initial values of u1, u2 and u3

3. Finding f1, f2, f3 and continuing with Predictor step

4. Finding dt using cfl condition

5. updating solution variables u1, u2 and u3

6. updating primitive variables and flux variables for corrector step.

7. corrector step

8. finding average residuals for solution variables

9. final update of solution variables.

10. Boundary conditions

11. Finding primitive variables.

12. Finding other dependent variables

13. plotting variables wrt x for diferent time steps

Function for initialization used in conservative form:

%initial values for conservative form
function [rho,t] = initialization(x,n,dx)
rho = zeros(1,n);
t = zeros(1,n);
for i = 1:(1 + 0.5/dx)
    rho(i) = 1;
    t(i) = 1;
end

for i = (1 + 0.5/dx):(1 + 1.5/dx)
    rho(i) = 1 - 0.366*(x(i)-0.5);
    t(i) = 1 - 0.167*(x(i) - 0.5);
end
    
for i = (1 + 1.5/dx):(1 + 3/dx)
    rho(i) = 0.634 - 0.3879*(x(i) - 1.5);
    t(i) = 0.833 - .3507*(x(i) - 1.5);
end

Explaination:

In the above code initial values for density and temperature are calculated for conservative form.

Results & Discussion:

I. For number of grid points = 31

1. Comparing time wise variation of density, velocity, temperature, pressure, mach number and mass flow rate for both conservative and non-conservative form:

Note: All the variables are non-dimensional

• From the above figure it can be observed that initially there is variation in variables (more variation in case of non-conservative form compared to conservative form).
• As the number of time iterations increases the variables tend to reach steady state at around 550 iterations.

2. comparing variation of density, velocity, temperature, pressure, mach number and mass flow rate wrt length for both conservative and non-conservative forms

Note: All the variables are non-dimensional.

• In the above figure steady state results of variation of non-dimensional variables wrt non-dimensional length is shown.
• it can be seen as expected that for subsonic-supersonic flow mach number increases from <1 to >1 .
• Flow is from high pressure at the inlet of nozzle to low pressure at the exit of nozzle, i.e, pressure energy is converted to kinetic energy. This inturn can be seen as velocity increases wrt x.
• Force is required for any gas to flow through nozzle this is usually provided by pressure ratio of inlet to exit.

3. Comparison of mass flow rate wrt length for conservative, non-conservative and exact solutions

note: mass flow rate is non dimensional mass flow rate.

From the above figure it can be clearly seen that when consevative form is used mass is conserved and only minor variations of mass flow rate can be observed which is more closer to exact solution compared to non-conservative form where more variation of mass flow rate is observed along the length.

4. comparison of mass flow rate wrt length for various time steps for non-conservative form:

Here variation of mass flow rate wrt length for different time iterations i.e 1,100,200,700,1400 shows how mass flow rate conveges to final answer and it can clearly be seen that mass flow rate converged at 700 iterations even for non-conservative form.

by using cfd technique mass flow rate found computationally is nearly same as the exact mass flow rate.

5. comparison of mass flow rate wrt length for various time steps for conservative form:

Here variation of mass flow rate wrt length for different time iterations i.e 1,100,200,700,1400 shows how mass flow rate conveges to final answer and it can clearly be seen that mass flow rate converged at 700 iterations even for conservative form but the values are closer compared to non-conservative form.

6. Residuals

• Residuals shows how fast and how much accurately the solution converges, they are called quality indexes for a given method.
• For non-conservative form residuals drop to an order of 1e-7 and for conservative form residuals drop to an order of 1e-5.
• This shows results of variable outputs from non-conservative form are more accurate than conservative form, but for practical purposes both the results are acceptable.

7. comparing computational time for both conservative and non-conservative forms

• Conservative form takes less time compared to non-conservative form.

II. Grid independency check:

Number of grid points	Form of GDE	Density	velocity	Temperature	mass flow rate	mach number	pressure
31	Non-conservative	0.6387	0.914	0.8365	0.5838	0.9994	0.5342
31	Conservative	0.6502	0.9	0.8399	0.5852	0.9892	0.5462
61	Non-conservative	0.6376	0.9137	0.8354	0.5825	0.9997	0.5326
61	conservative	0.6382	0.9133	0.8351	0.5829	0.9994	0.5330
	Exact analytical result	0.634		0.833		1.000	0.528

• From the above data it can be observed that by increasing number of grid points from 31 to 61 there is only slight change in values but solution is still grid independent for practical purpose and if high precision data is needed more number of grid points are recommended but not too high as there will also be round off errors.
• when number of grid points increased to 61 there is less difference between conservative and non-conservative variables compared to at n = 31.
• It can be clearly seen that as number of grid points increases solution marches nearer to the exact solution.

Conclusion:

1. From the above results it can be seen that conservative form does a better job in mass conservation along the length of the nozzle, where as in non-conservation form there are deviations.

2. For number of grid points = 31

Form of GDE	Density	Temperature	mach number	pressure
Non-conservative	0.6387	0.8365	0.9994	0.5342
Conservative	0.6502	0.8399	0.9892	0.5462
Exact analytical result	0.634	0.833	1.000	0.528

From the above table it can be clearly seen that the primitive variables in non-conservative case are close to exact analytical result as residuals in case of non-conservative case reach to 1e-7 for 1400 iterations compared to conservative case .

3. Non conservative form is easier to solve compared to conservative form.

4. But, Non conservative form takes more computational time compared to conservative form as observed from bar graph of result 7.

5. For practical purposes both conservative and non-conservative form results for this simulation are accepted. Both have their own advantages.

6. As the number of grid points increases more number of iterations need to be done to reach steady state. For n = 61 solution reaches steady state after 950 iterations.