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

1D Sub-sonic and Super-sonic Nozzle by Using McCormack Method

Objective:

To simulate the isentropic flow through a quasi 1D subsonic-supersonic nozzle by deriving both conservative and non-conservative forms of governing equation and solve them by using McCormack Method in MATLAB.

NOZZLE DIAGRAM:

The nozzle has three main sections

1. convergent.
2. Throat.
3. Divergent

In this diagram The Mach Number is describes the level/section of the Nozzle on the flow of sound velocity.

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.

• 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, the fig. shows, the flow at the inlet of the nozzle comes from a reservoir where pressure and temperature are denoted as byρ0ρ0 and T0T0 respectively. The cross-sectional area of the reservoir is large (theoretically, A) and hence the velocity is very small (V0). Thus and 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 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-conservative 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.

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 ρ'= ρ/ρ_0

Non-dimensional distance X' = X/L

Non-dimensional velocity V' = v/a_0

a_0 is local sonic velocity ( a_0 = √γRT_0 )

Non-dimensional time t' = t/((L/a_0 ))

L/a_0 is characteristic time scale

Non-dimensional Area A' = A/A^*

Non-Conservative form:

The non-conservative form has its differentiated variables as coefficients of some derivatives in the equation.

Continuity Form:

∂ρ'/∂t'= -ρ'(∂v^')/(∂x^' )-ρ'v'(∂(lnA^'))/(∂x^' ) -v'∂ρ'/∂x'

Momentum Form:

∂v'/∂t' =-v'∂v'/∂x'-1/γ(∂T'/∂x'+T'/ρ' ∂ρ'/∂x')

Energy Form:

∂T'/∂t' = -v'∂T'/∂x'-(γ-1)T'[∂v'/∂x'+v'(∂(lnA'))/∂x']

Where

T' = T/T_0 ; ρ' = ρ/ρ_0 ; x' = x/L ;

v' = v/a_0 ; A' = A/A^* ; t' = t/(L/a_0 ) ;

a_0 = Speed of sound

Conservative Form:

In the conservative form of equations, the variables that are being differentiated does not exist as coefficients of any derivative in the equation.

Continuity Equation:

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

Momentum Equation:

(∂(ρ'A'v'))/∂t' + (∂[ρ'A'〖v'〗^2+(1/γ)ρ'A'])/∂x' = 1/γ p'∂A'/∂x'

Energy Equation:

(∂[ρ'(e'/(γ-1)+γ/2 〖v'〗^2)A'])/∂t' + (∂[ρ'(e'/(γ-1)+γ/2)v'A'+ρ'A'v'])/∂x' = 0

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,

U_1 = ρ'A'

U_2 = ρ'A'v'

U_3 = ρ'(e'/(γ-1)+γ/2 〖v'〗^2)A'

F_1 = ρ'A'v'

F_2 = ρ'A'〖v'〗^2+1/γ p'A'

F_3 = ρ'(e'/(γ-1)+γ/2 〖v'〗^2)v'A'+p'A'v'

J_2 = 1/γ p'∂A'/∂x'

Would be written as

(∂U_1)/∂t' = -(∂F_1)/∂x' ; (∂U_2)/∂t' = -(∂F_2)/∂x'+J_2 ; (∂U_3)/∂t' = -(∂F_2)/∂x' ;

ρ' = U_1/A' ; V' = U_2/U_1 ; T'=e'= (γ-1)(U_3/U_1 -γ/2 〖v'〗^2) ; p' = ρ'T' ;

Adaptive time step method

∆t = minimum(∆〖t'〗_1,∆〖t'〗_2_ _ _ _ ∆〖t'〗_i_ _ _ _ _ ∆〖t'〗_n)

(〖∆t)〗_i^t = c ∆x/(a_i^t+v_i^t )

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, 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.

The Predictor step:

Initialization of flux:

(〖(∂F_1)/∂x)〗_i^t = (F_1^t (i+1)-F_1^t (i))/dx

(〖(∂F_2)/∂x)〗_i^t = (F_2^t (i+1)-F_2^t (i))/dx

(〖(∂F_3)/∂x)〗_i^t = (F_3^t (i+1)-F_3^t (i))/dx

Time step solution using space(flux) derivatives

(〖(∂U_1)/∂t)〗_i^t = -(〖(∂F_1)/∂x)〗_i^t

(〖(∂U_2)/∂t)〗_i^t = -(〖(∂F_2)/∂x)〗_i^t+J_2^t (i)

(〖(∂U_3)/∂t)〗_i^t = -(〖(∂F_3)/∂x)〗_i^t

Obtaining the intermediate value by using predictor/old value of time derivative

U_(1_i)^(-t+∆t) = U_(1_i)^t+(〖(∂U_1)/∂t)〗_i^t.∆t

U_(2_i)^(-t+∆t) = U_(2_i)^t+(〖(∂U_2)/∂t)〗_i^t.∆t

U_(3_i)^(-t+∆t) = U_(3_i)^t+(〖(∂U_3)/∂t)〗_i^t.∆t

Primitive values:

ρ' = U_1/A' T'=e'= (γ-1)(U_3/U_1 -γ/2 〖v'〗^2)

V' = U_2/U_1 p' = ρ' T'

Corrector Method:

(〖(∂F_1)/∂x)〗_i^(t+∆t) = F_1^(t+∆t) (i)-F_1^(t+∆t) (i-1)

(〖(∂F_2)/∂x)〗_i^(t+∆t) = F_2^(t+∆t) (i)-F_2^(t+∆t) (i-1)

(〖(∂F_3)/∂x)〗_i^(t+∆t) = F_3^(t+∆t) (i)-F_3^(t+∆t) (i-1)

Time step solution using space(flux) derivatives

(〖(∂U_1)/∂t)〗_i^(t+∆t) = -(〖(∂F_1)/∂x)〗_i^(t+∆t)

(〖(∂U_2)/∂t)〗_i^(t+∆t) = -(〖(∂F_2)/∂x)〗_i^(t+∆t)+J_2^(t+∆t) (i)

(〖(∂U_3)/∂t)〗_i^(t+∆t) = -(〖(∂F_3)/∂x)〗_i^(t+∆t)

Average value of predictor and corrector method

We get

({\frac{\partial\mathbit{U}_\mathbf{1}}{\partial t})}_{\mathbit{average}}\ =\ \mathbf{0}.\mathbf{5}[({\frac{\partial\mathbit{U}_\mathbf{1}}{\partial t})}_\mathbit{i}^\mathbit{t}+(U1∂t)it+∆t]

({\frac{\partial\mathbit{U}_\mathbf{2}}{\partial t})}_{\mathbit{average}}\ =\ \mathbf{0}.\mathbf{5}[({\frac{\partial\mathbit{U}_\mathbf{2}}{\partial t})}_\mathbit{i}^\mathbit{t}+(U2∂t)it+∆t]

({\frac{\partial\mathbit{U}_\mathbf{3}}{\partial t})}_{\mathbit{average}}\ =\ \mathbf{0}.\mathbf{5}[({\frac{\partial\mathbit{U}_\mathbf{3}}{\partial t})}_\mathbit{i}^\mathbit{t}+(U3∂t)it+∆t]

Final Solution:

Boundary condition:

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 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.

Initial and Boundary condition:

we must stipulate initial conditions for , T and V as a function of x , i.e. we must set up values of T and V at time t = 0.

Non-conservative Form:

Nozzle area A = 1+2.2(x+1.5) ^2 at ()

Initialize rho, v, & T

nt = 1400 n = 31 i = 1

Initial conditions:

= 1 – 0.3146x

T = 1 – 0.2314x

V = (0.1 + 1.09x)

Inflow boundary condition :

In our case, at i = 1(inlet), we have primitive variables like ρ, A, V and T.

we know the ρ and T at node 1 are equal to the reservoir's ρ and T which is assumed to be ρ = 1 and T = 1, so it becomes,

ρ (1) = 1 and 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.

ρ(n) = 2ρ(n-1) - ρ(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:

Initial conditions: Assumptions for initial profile of the nozzle

For 0 ≤ x' ≤ 0.5

ρ' = 1

T' = 1

For 0.5 ≤ x' ≤ 1.5

ρ' = 1.0 - 0.366(x'-0.5)

T' = 1.0 - 0.167(x'-0.5)

For 1.5 ≤ x' ≤ 3.5

ρ' = 0.634 - 0.3879(x'-1.5)

T' = 0.833 - 0.3507(x'-1.5)

Before stating the initial conditions in terms of solution vectors, we need to initialize velocity(V'). We have one of the dependent variables in our governing equation as , which is physically local mass flow; that is = ρ'A'V'. Therefore, for the initial conditions only; let us assume a constant mass flow through the nozzle and calculate V' as,

v\prime\ =\ \frac{U_2}{\rho\prime\ A\prime}\ \ =\ \frac{0.59}{\rho\prime\ A\prime}\

Note: The value 0.59 is chosen for , because it is close to exact analytical value of the steady-state mass flow (which is for this case is 0.579)

Therefore, the initial condition for V' as a function of x' is obtained by substituting the ρ′>`rho'` and A' variation. Finally, the initial conditions for U_1, U_2 and U_3 and are obtained by substituting the above variations for ρ' T' and V' into the definitions discussed above.

and the source term J as follows,

U_1\ =\ \rho\prime A\prime

U_2\ =\ \rho\prime A\prime v\prime

U_3\ =\ \rho\prime(\frac{e\prime}{\gamma-1}+\frac{\gamma}{2}{v\prime}^2)A\prime

where, e' = T'. Certainly, for the initial conditions described above, V' is calculated such that = ρ' A' V' = 0.59

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 ρ', A' and = ρ' 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.

Subsonic(inflow) boundary condition:

\mathbit{U}_\mathbf{1}(\mathbf{1})\ =\ \mathbit{\rho}\prime(\mathbf{1})\ \mathbit{A}\prime(\mathbf{1})\ \ =\ \mathbit{A}_{\mathbit{i}=\mathbf{1}}^\mathbit{t} Fixed value (Dirichlet boundary condition)

\mathbit{U}_\mathbf{2}(\mathbf{1})\ =\ \mathbf{2}\mathbit{U}_\mathbf{2}(\mathbf{2})-\mathbit{U}_\mathbf{2}(\mathbf{3})

\mathbit{U}_\mathbf{3}\ =\ \mathbit{U}_\mathbf{1}(\mathbf{1})(\frac{\mathbit{T}\prime(\mathbf{1})}{\mathbit{\gamma}-\mathbf{1}}+\frac{\mathbit{\gamma}}{\mathbf{2}}\ {\mathbit{v}\prime(\mathbf{1})}^\mathbf{2}) (Neumann boundary condition)

Outflow boundary condition :

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

\mathbit{U}_\mathbf{1}(\mathbit{N})\ =\ \mathbf{2}\mathbit{U}_\mathbf{1}(\mathbit{N}-\mathbf{1})\ -\ \mathbit{U}_\mathbf{1}(\mathbit{N}-\mathbf{2}) (Neumann boundary condition)

\mathbit{U}_\mathbf{2}(\mathbit{N})\ =\ \mathbf{2}\mathbit{U}_\mathbf{2}(\mathbit{N}-\mathbf{1})\ -\ \mathbit{U}_\mathbf{2}(\mathbit{N}-\mathbf{2}) (Neumann boundary condition)

\mathbit{U}_\mathbf{3}(\mathbit{N})\ =\ \mathbf{2}\mathbit{U}_\mathbf{3}(\mathbit{N}-\mathbf{1})\ -\ \mathbit{U}_\mathbf{3}(\mathbit{N}-\mathbf{2}) (Neumann boundary condition)

MAIN PROGRAM

close all
clear all
clc
%input variables
%valve of n can be varied by our own as it is grid size
n = 31;
L = 3;
x = linspace(0,L,n);
%Courent no.
cfl = 0.5;
%nodes at throat
th = ((n-1)/2)+2;
% time steps
nt = 1400;

tic
%Function for non conservative form
[v_nc,rho_nc,T_nc,mass_nc,M_nc,p_nc,v_throat_nc,rho_throat_nc,T_throat_nc,mass_throat_nc,M_throat_nc,p_throat_nc] = non_conservative_approach(n,L,cfl,nt,th);
% Calculating the simulation time
Non_conservative_simulation_time = toc;

tic
% Function for conservative form
[v_c,rho_c,T_c,mass_c,M_c,p_c,v_throat_c,rho_throat_c,T_throat_c,mass_throat_c,M_throat_c,p_throat_c]=conservative_approach(n,L,cfl,nt,th);
%Calculating the simulation time
conservative_simulation_time = toc;

%Analatical Solution Steady State
rho_an = 0.634*ones(900,1400);
T_an = 0.833*ones(900,1400);
p_an = 0.528*ones(900,1400);
M_an = 1*ones(900,1400);

figure(5)
a1 = subplot(4,1,1);
plot(linspace(1,nt,nt),rho_throat_nc);
hold on 
plot(linspace(1,nt,nt),rho_throat_c,'g');
hold on
plot(linspace(1,nt,nt),rho_an,'k');
hold on
grid on
ylabel('Density')
legend('Non-conservative approch','Conseravative approach','Analatical approach')
title('Variation in parameter at throat w.r.t time step')

a2 = subplot(4,1,2);
plot(linspace(1,nt,nt),T_throat_nc);
hold on 
plot(linspace(1,nt,nt),T_throat_c);
hold on
plot(linspace(1,nt,nt),T_an);
hold on
grid on
ylabel('Temperature')
legend('Non-conservative approch','Conseravative approach','Analatical approach')
title('Variation in parameter at throat w.r.t time step')

a3 = subplot(4,1,3);
plot(linspace(1,nt,nt),p_throat_nc);
hold on 
plot(linspace(1,nt,nt),p_throat_c);
hold on
plot(linspace(1,nt,nt),p_an);
hold on
grid on
ylabel('Pressure')
legend('Non-conservative approch','Conseravative approach','Analatical approach')
title('Variation in parameter at throat w.r.t time step')

a4 = subplot(4,1,4);
plot(linspace(1,nt,nt),M_throat_nc);
hold on 
plot(linspace(1,nt,nt),M_throat_c);
hold on
plot(linspace(1,nt,nt),M_an);
hold on
grid on
ylabel('Mach number')
legend('Non-conservative approch','Conseravative approach','Analatical approach')
title('Variation in parameter at throat w.r.t time step')

figure(6)
b1 = subplot(4,1,1);
plot(x,rho_nc);
hold on 
plot(x,rho_nc,'*');
hold on
grid on
ylabel('Density')
legend('Non-conservative approch','Conseravative approach','Analatical approach')
title('Variation in properties along X')

b2 = subplot(4,1,2);
plot(x,T_nc);
hold on 
plot(x,T_nc,'*');
hold on
grid on
ylabel('Temperature')
legend('Non-conservative approch','Conseravative Solution')
title('Variation in properties along X')

b3 = subplot(4,1,3);
plot(x,p_nc);
hold on 
plot(x,p_nc,'*');
hold on
grid on
ylabel('Pressure')
legend('Non-conservative approch','Conseravative Solution')

b4 = subplot(4,1,4);
plot(x,M_nc);
hold on 
plot(x,M_nc,'*');
hold on
grid on
ylabel('Mach Number')
legend('Non-conservative approch','Conseravative Solution')

%mass flow rate comparission between two forms of equations
figure(7)
plot(x,mass_nc)
hold on 
plot(x,mass_c,'g')
hold on 
plot(x,linspace(0.579,0.579,n),'k')
hold on 
ylabel('non-dimentionless mass flow rate')
legend('Non-conservative approch','Conseravative approach','Analatical approach')
title('non-dimentionless mass flow rate comparission')
axis([0 3 0.55 0.6])
grid on

y = [rho_nc(th) rho_c(th) 0.634; T_nc(th) T_c(th) 0.833; p_nc(th) p_c(th) 0.528; M_nc(th) M_c(th) 1];
x_axis = categorical({'1.denisty','2.Temperature''3.Pressure','4.mach number'});

figure(10)
labels = {'Non-conservative approch','Conseravative approach','Analatical approach'};
b = bar(x_axis,y,'group');
legend(labels,'location','northoutside','Orientation','Vertical')
a = (1:size(y,1)).';
x_ax = [a-0.25 a a+0.25];
for k = 1:size(y,1)
    for m = 1:size(y,2)
        text(x_ax(k,m),y(k,m),num2str(y(k,m),'%0.2f'),'HorizontalAlignment','center','VerticalAlignment','bottom')
    end
end





NON-CONSERVATIVE FORM

function[v,rho,T,mass,M,p,v_throat,rho_throat,T_throat,mass_throat,M_throat,p_throat]=non_conservative_approach(n,L,cfl,nt,th)
%inputs
n=31;
x=linspace(0,3,n);
dx=x(2)-x(1);

%conidtions
rho=1-0.3146*x;
T=1-0.2314*x;
v=(0.1+1.09*x).*T.^0.5;
a=1+2.2*(x-1.5).^2;
th=find(a==1);
gama=1.4;

%time steps
%Time Step calculation based on Courant number.
cfl=0.5;
dt = min(cfl.*(dx./((sqrt(T))+v)));
%Number of iterations
nt=1400;
%rho_initial=rho;
%T_initial=T;
%v_initial=v;

for k=1:nt
%copy of the variables
rho_old=rho;
T_old=T;
v_old=v;
    
%predictor step
    for j=2:n-1
        dv_dx= (v(j+1)-v(j))/dx;
        drho_dx =(rho(j+1)-rho(j))/dx;
        dTdx = (T(j+1)-T(j))/dx;
        da_dx =(log(a(j+1))-log(a(j)))/dx;
        %continuity
        drho_dt_p(j) = -rho(j)*dv_dx - rho(j)*v(j)*da_dx - v(j)*drho_dx;
        %momentum
        dv_dt_p(j) = -v(j)*dv_dx - (1/gama)*(dTdx+(T(j)/rho(j))*drho_dx);
        %energy
        dT_dt_p(j) = -v(j)*dTdx -(gama-1)*T(j)*(dv_dx+v(j)*da_dx);
       %solution updates
       v(j)= v(j)+dv_dt_p(j)*dt;
       T(j)= T(j)+dT_dt_p(j)*dt;
       rho(j)= rho(j)+drho_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;
        dTdx = (T(j)-T(j-1))/dx;
        da_dx =(log(a(j))-log(a(j-1)))/dx;
        %continuity
        drho_dt_c(j) = -rho(j)*dv_dx-rho(j)*v(j)*da_dx-v(j)*drho_dx;
        %momentum
        dv_dt_c(j) = -v(j)*dv_dx-(gama^-1)*(dTdx+(T(j)/rho(j))*drho_dx);
        %energy
        dT_dt_c(j) = -v(j)*dTdx -(gama-1)*T(j)*(dv_dx+v(j)*da_dx);
    end
    
    %average
    dv_dt_avg=0.5*(dv_dt_c+dv_dt_p);
    drho_dt_avg = 0.5*(drho_dt_c+drho_dt_p);
    dT_dt_avg = 0.5*(dT_dt_p+dT_dt_c);
    
    % Solution update
    for i =2:n-1
     v(i)= v_old(i)+dv_dt_avg(i)*dt;
     T(i)= T_old(i)+dT_dt_avg(i)*dt;
     rho(i)= rho_old(i)+drho_dt_avg(i)*dt;
    end
    
    %Inlet condition by interpolation
    v(1)=2*v(2)-v(3);
    %outlet condition by interpolation
    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);
     
   %properties
   mass = rho.*a.*v;
   M = v./((T).^0.5);
   p= rho.*T;
   
   %throat values
   rho_throat(k)=rho(th);
   v_throat(k)=v(th);
   T_throat(k)=T(th);
   M_throat(k)=M(th);
   p_throat(k)=p(th);
   mass_throat(k)=mass(th);
   
   %residues at throat.
   v_res(k)=dv_dt_avg(th);
   rho_res(k)=drho_dt_avg(th);
   T_res(k)=dT_dt_avg(th);
   
%plotting the results
% plotting mass flow rate at different time steps
figure(1)
if k == 50
plot(x, mass, 'b', 'linewidth', 1.5)
hold on
elseif k == 100
plot(x, mass, 'r', 'linewidth', 1.5)
hold on
elseif k == 150
plot(x, mass, 'c', 'linewidth', 1.5)
hold on
elseif k == 200
plot(x, mass, 'g', 'linewidth', 1.5)
hold on
elseif k == 550
plot(x, mass, 'y', 'linewidth', 1.5)
hold on
elseif k == 700
plot(x, mass, 'bl', 'linewidth', 1.5)
hold on
elseif k == 1000
plot(x, mass, 'k', 'linewidth', 1.5)
hold on
elseif k == 1400
plot(x, mass,'m', 'linewidth', 1.5)
grid on
ylabel('mass flow rate');
xlabel('nozzle length (x)');
title('mass flow rate at different time steps Non conservative approach');
legend('50 time steps','100 time steps','150 time steps','200 time steps','550 time steps','700 time steps','1000 time steps','1400 time steps');
end
end

figure(3)
plot(linspace(1,nt,nt),v_res);
hold on
plot(linspace(1,nt,nt),rho_res,'g');
hold on
plot(linspace(1,nt,nt),T_res,'k');
grid on
ylabel('residues');
xlabel('timestep');
legend('vplot','rhoplot','Tplot');
title('residual plot');
end






CONSERVATIVE FORM

function[v,rho,T,mass,M,p,v_throat,rho_throat,T_throat,mass_throat,M_throat,p_throat]=conservative_approach(n,L,cfl,nt,th)
% geometry
L=3;
n=31;
x=linspace(0,L,n);
dx= x(2)-x(1);

%Initial Conditions 
a=1+2.2*(x-1.5).^2;
th=find(a==1);  %thraot node.
rho_initial=ones(1,n);
T_initial=ones(1,n);


for i= 1:n
    if (x(i)<=0.5)
        rho_initial(i)=1;
        T_initial(i)=1;
    elseif (x(i)>=0.5 && x(i)<=1.5)
            rho_initial(i)=1-0.366*(x(i)-0.5);
            T_initial(i)=1-0.167*(x(i)-0.5);
    elseif (x(i)>=1.5 && x(i)<=3)
           rho_initial(i)=0.634-0.3879*(x(i)-1.5);
            T_initial(i)=0.833-0.3507*(x(i)-1.5);
       
  
    end
end

v = 0.59./(rho_initial.*a);
gamma=1.4;

%Solution term
p= rho_initial.*T_initial;
u1 = rho_initial.*a;
u2 = rho_initial.*a.*v;
u3 = u1.*((T_initial/(gamma-1))+(0.5*gamma*v.^2));


%Time Step calculation based on Courant number.
cfl=0.5;
dt = min(cfl.*(dx./((sqrt(T_initial))+v)));
%Number of iterations
nt=1400;

rho=rho_initial;
T= T_initial;


%time loop

for k=1:nt
    
   %copy of solution terms
    u1old=u1;
    u2old=u2;
    u3old=u3;
    
    
   %Flux terms 
   
   f1= u2;
   term=(((gamma-1)/gamma)*(u3-((gamma/2)*((u2.^2)./u1))));
   f2=((u2.^2)./u1)+ term;
   f3 = ((gamma*u2.*u3)./u1)-(gamma*(gamma-1)*(u2.^3))./(2*u1.^2);
   
   %predictor step
   for i =2:n-1
       df1dx(i)=(f1(i+1)-f1(i))/dx;
       df2dx(i)=(f2(i+1)-f2(i))/dx;
       df3dx(i)=(f3(i+1)-f3(i))/dx;
       dlogadx(i) = (log(a(i+1))-log(a(i)))/dx;
       
       dadx(i) = (a(i+1)-a(i))/dx;
       j2(i)=(1/gamma)*rho(i)*T(i)*dadx(i);
       
       %continuity equation
       du1dt_p(i)= -df1dx(i);
       %momentum equation
       du2dt_p(i)= -df2dx(i)+j2(i);
       %energy equation.
       du3dt_p(i)= -df3dx(i);
       
       %solution update.
       u1(i)=u1(i)+du1dt_p(i)*dt;
       u2(i)=u2(i)+du2dt_p(i)*dt;
       u3(i)=u3(i)+du3dt_p(i)*dt;  
   end
   
   %primitive variables
   rho = u1./a;
   v = u2./u1;
   T= (gamma-1)*((u3./u1)-((gamma*v.^2)/2));
   
   %flux terms before corrector step
   f1= u2;
   term=(((gamma-1)/gamma)*(u3-(gamma*u2.^2)./(2*u1)));
   f2=((u2.^2)./u1)+ term;
   f3 = ((gamma*u2.*u3)./u1)-((gamma*(gamma-1)/2)*((u2.^3)./(u1.^2)));
   
   
   %corrector method
   for i=2:n-1
       df1dx(i)=(f1(i)-f1(i-1))/dx;
       df2dx(i)=(f2(i)-f2(i-1))/dx;
       df3dx(i)=(f3(i)-f3(i-1))/dx;
       dlogadx(i) = (log(a(i+1))-log(a(i)))/dx;
       dadx(i) = (a(i)-a(i-1))/dx;
       j2(i)=(1/gamma)*rho(i)*T(i)*dadx(i);
       
       %continuity equation
       du1dt_c(i)= -df1dx(i);
       %momentum equation
       du2dt_c(i)= -df2dx(i)+j2(i);
       %energy equation
       du3dt_c(i)= -df3dx(i);      
   end
   
   %avergae derivative
   du1dt = 0.5*(du1dt_c+du1dt_p);
   du2dt = 0.5*(du2dt_c+du2dt_p);
   du3dt = 0.5*(du3dt_c+du3dt_p);
   
   %final solution update.
   for j =2:n-1
   u1(j)=u1old(j)+du1dt(j)*dt;
   u2(j)=u2old(j)+du2dt(j)*dt;
   u3(j)=u3old(j)+du3dt(j)*dt;
   end
   
  %boundary conditions
  %inlet condition
        u1(1) = rho(1)*a(1);
        u2(1)= 2*u2(2)-u2(3);
        v(1) = u2(1)/u1(1);
        u3(1)= u1(1)*((T(1)/(gamma-1))+(0.5*gamma*v(1)^2));
        
        % outle Boundary
        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);
  
        %primitive terms
        rho = u1./a;
        v = u2./u1;
        T= (gamma-1)*((u3./u1)-(gamma/2)*(v.^2));
   
        %properties
        mass = rho.*a.*v;
        M = v./((T).^0.5);
        p= rho.*T;
   
       %throat values
       rho_throat(k)=rho(th);
       v_throat(k)=v(th);
       T_throat(k)=T(th);
       M_throat(k)=M(th);
       p_throat(k)=p(th);
       mass_throat(k)=mass(th);
   
       %residues at throat.
       u1_res(k)=du1dt(th);
       u2_res(k)=du2dt(th);
       u3_res(k)=du3dt(th);
 
%plotting the results 
% plotting mass flow rate at different time steps
figure(2);
if k == 50
plot(x, mass, 'b','linewidth', 1.5)
hold on
elseif k == 100
plot(x, mass, 'r','linewidth', 1.5)
hold on
elseif k == 150
plot(x, mass, 'c','linewidth', 1.5)
hold on
elseif k == 200
plot(x, mass, 'g','linewidth', 1.5)
hold on
elseif k == 550
plot(x, mass, 'y', 'linewidth', 1.5)
hold on
elseif k == 700
plot(x, mass, 'k', 'linewidth', 1.5)
hold on
elseif k == 1400
plot(x, mass,'m','linewidth', 1.5)
grid on
ylabel('mass flow rate');
xlabel('nozzle length (x)');
title('mass flow rate at different time steps conservative appraoch');
legend('50 time steps','100 time steps','150 time steps','200 time steps','550 time steps','700 time steps','1400 time steps');
  
 end
end


figure(4)
plot(linspace(1,nt,nt),u1_res);
hold on
plot(linspace(1,nt,nt),u2_res,'g');
hold on
plot(linspace(1,nt,nt),u3_res,'k');
grid on
ylabel('residues');
xlabel('timestep');
legend('u1plot','u2plot','u3plot');
title('residual plot');
end

Figure 1.For non-conservative form the mass flow rate v/s nozzle length the time step is varies at 50 and it can be almost constant at higher time step(1400).

Figure 3.In this plot the velocity plot may increasing with residual and the temperature plot may decreases with residual and for density at certain time step varies the residual at it be constant .

Figure 2.For conservative in mass flow rate v/s nozzle length the lesser time step may have the high mass flow rate and the highr time step may heve the constant mass flow rate w.r.t. nozzele length.

Figure 4.In this plot the U3 plot varies more and it is unstable at starting condition and other two U1 & U2 plot may be stable condition.

Figure 5. In above fig. the non-conservative approach is unstable for density, temperature and pressure. For mach number the approach is higher to lesser. The conservative and the analytical approaches are stable.

Figure 6. The fig. shows the density, temperature and pressure are decrease at throat and the mach number increase at outlet.

Figure 7. In above plot the non-conservative approach is unstable. The conservatuve and analytical approaches are in stable condition.

Figure 10. The bar graph shows the density, temperature, pressure, and mach number to non-conservative, conservative and analytical approaches.

Primitive Grid test values:

Flow variable at throat at n=31	Density	Temperature	Pressure	Mach Number
Non-conservative form	0.6387	0.8365	0.5342	0.9994
conservative form	0.6524	0.8410	0.5487	0.9780
Flow variable at throat at n=61
Non-conservative form	0.6377	0.8354	0.5327	0.9998
conservative form	0.6391	0.8355	0.5340	0.9979
Analytical solution	0.6340	0.8330	0.5280	1

From the table it is clear that the solution of 61 grid points almost equal to the exact solution when compare to that of the solution of 31 grid points. More the number of grid points more will the computational time but at the same time we can able to produce the accurate result than the lesser grid point.

Conclusion:

Simulation of isentropic flow through a quasi 1D subsonic-supersonic nozzle by deriving both conservative and non-conservative forms of governing equation and solve them by using MacCormack Method in MATLAB is successfully compiled and the results are discussed at the output. Between conservative and non-conservative forms, in order to get a smooth solution conservative form is preferred because non-conservative form may leads to shocks at the shock flow region.