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

Aim:

To simulate 1D Super-sonic nozzle flow using the MacCormack Method

Procedure:

In order to simulate the nozzle flow, we must know the governing equation of the solving method. In this project, I have divided into different sections based on the process of the work. To simulate the equation, first, we need to know the governing equation, followed by making a suitable equation for the solution. Identify the initial and boundary condition for solving the equation. By solving the equation, using the analytical method to identify the equation obtained and values are correct. Further, we going to convert the equation to the code form to simulate the answer. In final we post-process the results for the conclusion. Finally, present the results.

Formulation of governing Equation

The first step to understanding the physical aspects of the flow field.

Let us consider`P_o, T_o, rho_o, (V=0)`are the initial condition of the flow from the reservoir. The area at the reservoir is considered as infinite large so that the velocity of the flow at the initial condition is approaching to be zero.

`P_o, T_o, rho_o` is the stagnation properties of the flow that are pressure, Temperature, and density respectively.

The initial condition of the flow is

• We consider the flow is steady, isentropic through the Convergent-Divergent Nozzle.
• The flow expands isentropically to the supersonic speeds at the nozzle exit. Where the properties are given by `P_e, T_e, V_e, A_e and M_e`
• The flow reaches the sonic speed at the throat of the nozzle (`M=1`), Where the area is minimum. The properties at the throat are given by`A^**, P^**, T^**, and V^**=a^**`

Assumption: In general the flow changes with two-dimension(`x and y`). For this, we make the assumption the properties are changes only vary with the `x`. this is tantamount to assuming uniform flow properties across any given cross-section. Such flow is called Quasi-dimensional flow.

The relation between the area and velocity of the C-D nozzle is given by

In the base of combine with mass

`(1-M^2)(dV)/V=-(dA)/A`

For subsonic flow (M<1), increase in area (dA>0), Cause the velocity to decrease (dV<0).

For supersonic flow (M>1), a decrease in the area (dA<0), Cause the flow velocity to increase (dV>0).

Governing Equation:

The equation for this quasi-one-dimensional, steady, isentropic flow expressed as

Continuity Equation:

`rho_1*V_1*A_1=rho_2*V_2*A_2`

Momentum equation:

`p_1*A_1+rho_1*V_1^2*A_1+int_(A1)^(A2) p.dA=p_2*A_2+rho_2*V_2^2*A_2`

Energy Equation:

`h_1+(V_1^2)/2=h_2+(V_2^2)/2`

Perfect gas Equation of state:

`p=rho*R*T`

The relation for a calorically perfect gas,

`h=c_pT`

The above equations are solved analytically for the flow through the nozzle.The area ratio `A/(A^**)` given the Mach Number Variation through the relation as follow

`A/(A^**)=1/M^2[2/(gamma+1)(1+((gamma-1)/2)M^2)]^((gamma+1)/(gamma-1))`

The Other properties are given by

`p/p_o=(1+((gamma-1)/2)M^2)^-((gamma)/(gamma-1))` Pressure Ratio

`rho/rho_o=(1+(gamma-1)/2M^2)^(-1/(gamma-1))` Density Ratio

`T/T_o=(1+(gamma-1)/(2)M^2)^-1` Temperature Ratio

The physical principles are

• Mass conservation
• Newton's Second law
• Energy conservation

To ensure that these physical principles are satisfied, we must return to the integral forms of governing equation and apply these integral forms to a control volume consistent with quasi-one-dimensional assumption.

The integral form of the continuity equation is given by

`del/(delt)intintint_v rhodv+intint_srhoV.dS=0`

On simplifying above equation

`(del(rhoA))/(delt)+(del(rhoAV))/(delx)=0` (1a)

The integral form of the `x`component of the momentum equation, the equation for inviscid flow(neglecting the viscous stress terms) with no body forces,

`del/(delt)intintint_v(rhou)dv+intint_S(rhouV).dS=-intint_S(p*dS)_x`

On simplifying the above equation

`rho(delV)/(delt)+rhoV(delV)/(delx)=-(delp)/(delx)` (1b)

The integral form of the energy equation for an adiabatic flow(`dot q=0`) with no body forces and no viscous effects is given by,

`del/(delt)intintint_vrho(e+V^2/2)dv+intint_Srho(e+V^2/2)V.dS=-intint_s(p*V)*dS`

On simplifying the above equation

`rho(dele)/(delt)+rhoV(dele)/(delx)=-p(delV)/(delx)-pV(del(lnA))/(delx)` (1c)

On further simplifying the 1a, 1b and 1c equation using the calorically perfect gas equation and the pressure is eliminated from the equation of the state

`e=c_vT`

`p=rhoRT`

Continuity:

`(del(rhoA))/(delt)+rhoA(delV)/(delx)+rhoV(delA)/(delx)+VA(delrho)/(delx)=0`

Momentum:

`rho(delV)/(delt)+rhoV(delV)/(delx)=-R(rho(delT)/(delx)+T(delrho)/(delx))`

Energy:

`rhoc_v(delT)/(delt)+rhoVc_v(delT)/(delx)=-rhoRT[(delV)/(delx)+V(del(ln A))/(delx)]`

The nozzle flow and flow-field variables are frequently expressed in terms of non-dimensional variables,with reference to the reservoir values. The non-dimensional variables `p/(p_o),rho/rho_o, and T/T_o`

We define the non-dimensional temperature and density are

`T'=T/T_o`

`rho'=rho/rho_o`

We define a dimensionless length as

`x'=x/L`

Denoting the speed of the sound in the reservoir as `a_o`,

`a_o=sqrt(gammaRT_o)`

Dimensionless velocity as

`V'=V/a_o`

Dimensionless time and the area with throat area `A^**`

`t'=t/(L/a_o)`

`A'=A/(A^**)`

On simplifying using the terms we get,

Continuity equation:

`(delrho')/(delt')=-rho'(delV')/(delx')-rho'V'(del(lnA'))/(delx')-V(delrho')/(delx')` (1d)

Momentum Equation:

`(delV')/(delt')=-V'(delV')/(delx')-1/(gamma)((delT')/(delx')+(T')/(rho')*(delrho')/(delx'))` (1e)

Energy equation:

`(delT')/(delt')=-V'(delT')/(delx')-(gamma-1)T'[(delV')/(delx')+V'(del(lnA'))/(delx')]` (1f)

Applying the MacCormack's technique:

The MacCormack method is a variation of the two steps of Lax-Wendroff but is much simpler in application. The application of this technique proceeds in two steps technique namely the predicator and collector steps

In predicator steps, the provisional value for the selected variable is predicted using the forward differencing. In collector steps the selected value is corrected by the using the backward differencing.

Above the equation 1d,1e, and 1f are solved by the MacCormack's technique.

Before the explanation of the technique.

One can solve the problem in two ways. they are

• Conservative form
• Non-Conservative form

In this challenge, We are going to both forms of the method to solve the equation.

Note: Conservative and Non-Conservative forms are Mathematically the same, but the variables that are being differentiated does not exist as coefficients of any derivative in the equation. On the other hand, the non-conservative form has its differentiated variables as coefficients of some derivatives in the equation.

As part of the solution, we need to get the equation for both forms from the above derivation part

For Conservative form:

Continuity equation:

`(del(rho'A'))/(delt')+(del(rho'A'V'))/(delx')=0` (1g)

Momentum equation:

`(del(rho'A'V'))/(delt')+(del(rho'A'V'^2+(p'A')/gamma))/(delx')=1/(gamma)p'(delA')/(delx')` (1h)

Energy equation:

`(del[rho'{(T')/(gamma-1)+(gamma)/2V'^2}A'])/(delt')+(del[rho'{(T')/(gamma-1)+(gamma)/2V'^2}V'A'+p'A'V'])/(delx')=0` (1i)

For Non-Conservative form:

Continuity equation:

`(delrho')/(delt')=-rho(delV')/(delx')-rho'V'(del(lnA'))/(delx')-V'(delrho')/(delx')` (1j)

Momentum equation:

`(delV')/(delt')=-V'(delV')/(delx')-1/(gamma)((delT')/(delx')+(T')/(rho')(delrho')/(delx'))` (1k)

Energy equation:

`(delT')/(delt')=-V'(delT')/(delx')-(gamma-1)T'[(delV')/(delx')+V'(del(ln A))/(delx')]` (1L)

The above equation are discretized using the MacCormack Technique

NON-CONSERVATIVE FORM:

PREDICATOR STEP:

In this step, the equation are solved by the predicting the value for the density, velocity, and Temperature by the Forward difference Method.Once we find the value (Temporary derivative) is updates in the next time step.

Continuity equation:

`((delrho)/(delt))_(p(i))^t=-rho_i^t(V_(i+1)^t-V_(i)^t)/(Deltax)-rho_i^tV_i^t(ln A_(i+1)-ln A_(i))/(Deltax)-V_i^t(rho_(i+1)^t-rho_(i)^t)/(Deltax)`

Momentum equation:

`((delV)/(delt))_(p(i))^t=V_i(V_(i+1)^t-V_(i)^t)/(Deltax)-1/(gamma)[(T_(i+1)^t-T_i^t)/(Deltax)+(T_i^t)/(rho_i^t)(rho_(i+1)^t-rho_(i)^t)/(Deltax)]`

Energy equation:

`((delT)/(delt))_(p(i))^t=-V_i^t(T_(i+1)^t-T_(i)^t)/(Deltax)-(gamma-1)T_i^t[(V_(i+1)^t-V_(i)^t)/(Deltax)+V_i^t(lnA_(i+1)-lnA_i)/(Deltax)]`

After Finding the values using the above equation, the value are predicted by the following equations

Density given by,

`rho_i^(t+Deltat)=rho_i^t+((delrho)/(delt))_(p(i))^t(Deltat)`

Velocity is given by,

`V_i^(t+Deltat)=V_i^t+((delV)/(delt))_(p(i))^t(Deltat)`

Temperature is given by,

`T_i^(t+Deltat)=T_i^t+((delT)/(delt))_(p(i))^t(Deltat)`

CORRECTOR STEP:

In this step, the predicted values of density,Temperature, and velocity from the above step are collected by using backward differencing to find out the correct value and updating to the next time step.

Continuity equation:

`((delrho)/(delt))_(c(i))^(t+Deltat)=-rho_i^(t+Deltat)(V_(i)^(t+Deltat)-V_(i-1)^(t+Deltat))/(Deltax)-rho_i^(t+Deltat)V_i^(t+Deltat)(ln A_(i)-ln A_(i-1))/(Deltax)-V_i^(t+Deltat)(rho_(i)^(t+Deltat)-rho_(i-1)^(t+Deltat))/(Deltax)`

Momentum equation:

`((delV)/(delt))_(c(i))^(t+Deltat)=V_i(V_(i)^(t+Deltat)-V_(i-1)^(t+Deltat))/(Deltax)-1/(gamma)[(T_(i)^(t+Deltat)-T_(i-)^(t+Deltat))/(Deltax)+(T_i^t)/(rho_i^t)(rho_(i)^(t+Deltat)-rho_(i-1)^(t+Deltat))/(Deltax)]`

Energy Equation:

`((delT)/(delt))_(c(i))^(t+Deltat)=-V_i^(t+Deltat)(T_(i)^(t+Deltat)-T_(i-1)^(t+Deltat))/(Deltax)-(gamma-1)T_i^(t+Deltat)[(V_(i)^(t+Deltat)-V_(i-1)^(t+Deltat))/(Deltax)+V_i^(t+Deltat)(lnA_(i+1)-lnA_i)/(Deltax)]`

After finding the collector temporal values, we have to update the average of temporal terms calculated from predictor and collector steps

Density:

`((delrho)/(delt))_(avg)=0.5[((delrho)/(delt))_(p(i))^t+((delrho)/(delt))_(c(i))^(t+Deltat)]`

Velocity:

`((delV)/(delt))_(avg)=0.5[((delV)/(delt))_(p(i))^t+((delV)/(delt))_(c(i))^(t+Deltat)]`

Temperature:

`((delT)/(delt))_(avg)=0.5[((delT)/(delt))_(p(i))^t+((delT)/(delt))_(c(i))^(t+Deltat)]`

Finally, we update the properties to the next time step by the following equation

Density:

`rho_i^(t+Deltat)=rho_i^t+((delrho)/(delt))_(avg)^t(Deltat)`

Velocity:

`V_i^(t+Deltat)=V_i^t+((delV)/(delt))_(avg)^t(Deltat)`

Temperature:

`T_i^(t+Deltat)=T_i^t+((delT)/(delt))_(avg)^t(Deltat)`

CONSERVATIVE FORM:

Conservative form is solved by using some simple form

From J.A's book

U,F,J are some generic form is used to solve the equation and they are respresented as solution vector,flux vector, and source term

The terms are

`U_1=rho'A'`

`U_2=rho'A'V'`

`U_3=rho'((T')/(gamma-1)+(gamma)/2V'^2)A'`

FLUX TERMS are given by

`F_1=rho'A'V'`

`F_2=rho'A'V'+1/(gamma)p'A'`

`F_3=rho'((T')/(gamma-1)+(gamma)/2V'^2)V'A'+p'A'V'`

Source terms:

`J_2=1/(gamma)p'(delA')/(delx')`

With these equation, we can simplify the equation.,

Continuity equation

`(delU_1)/(delt')=-(delF_1)/(delx')`

Momentum equation

`(delU_2)/(delt')=-(delF_2)/(delx')+J_2`

Energy equation

`(delU_3)/(delt')=-(delF_3)/(delx')`

Here the above is solved by the Mac-Cormack method:

1.Predictor step:

This step is used to predict the temporal value for the terms using the forward differencing equation,

Continuity equation:

`((delU_1)/(delt))_(p(i))^t=-((F_1)_(i+1)^t-(F_1)_i^t)/(Deltax)`

Momentum equation:

`((delU_2)/(delt))_(p(i))^t=-((F_2)_(i+1)^t-(F_2)_i^t)/(Deltax)+(J_2)_i^t`

Energy Equation:

`((delU_3)/(delt))_(p(i))^t=-((F_3)_(i+1)^t-(F_3)_i^t)/(Deltax)`

We proceed to find the predicted values of `U_1,U_2, and U_3`,is calculated by the following reaction

`(U_1)_(p(i))^(t+Deltat)=(U_1)_i^t+((delU_1)/(delt))_(p(i))^t(Deltat)`

`(U_2)_(p(i))^(t+Deltat)=(U_2)_i^t+((delU_2)/(delt))_(p(i))^t(Deltat)`

`(U_3)_(p(i))^(t+Deltat)=(U_3)_i^t+((delU_3)/(delt))_(p(i))^t(Deltat)`

After this stage the predicted values are decoded from`U_1,U_2, and U_3` using the following equations

Density:

`(rho')_(p(i))^(t+Deltat)=((U_1)_(p(i))^(t+Deltat))/(A')`

Velocity:

`(V')_(p(i))^(t+Deltat)=((U_1)_(p(i))^(t+Deltat))/((U_2)_(p(i))^(t+Deltat))`

Temperature:

`(T')_(p(i))^(t+Deltat)=(gamma-1)[((U_1)_(p(i))^(t+Deltat))/(U_2)_(p(i))^(t+Deltat)-gamma/2((V')_(p(i))^(t+Deltat))]`

Pressure:

`p'=rho'T'`

Corrector step-Conservative

In this step, the updated value of the predicted step are used for the backward differencing spacial terms which is used to calculate the corrector temporal terms and further calculate the average terms and update in the next time step.

Discretization of equation

Continuity Equation:

`((delU_1)/(delt))_(c(i))^(t+Deltat)=-((F_1)_i^(t+Deltat)-(F_1)_(i-1)^(t+Deltat))/(Deltax)`

Momentum Equation:

`((delU_2)/(delt))_(c(i))^(t+Deltat)=-((F_2)_i^(t+Deltat)-(F_2)_(i-1)^(t+Deltat))/(Deltax)+(J_2)_i^(t+Deltat)`

Energy Equation:

`((delU_3)/(delt))_(c(i))^(t+Deltat)=-((F_3)_i^(t+Deltat)-(F_3)_(i-1)^(t+Deltat))/(Deltax)`

Average the terms for the predicted and collector time step

`((delU_1)/(delt))_(avg)=0.5[((delU_1)/(delt))_(p(i))^t+((delU_1)/(delt))_(c(i))^(t+Deltat)]`

`((delU_2)/(delt))_(avg)=0.5[((delU_2)/(delt))_(p(i))^t+((delU_2)/(delt))_(c(i))^(t+Deltat)]`

`((delU_3)/(delt))_(avg)=0.5[((delU_3)/(delt))_(p(i))^t+((delU_3)/(delt))_(c(i))^(t+Deltat)]`

After averaging, we need to update to the next time step

`(U_1)_(i)^(t+Deltat)=(U_1)_i^t+((delU_1)/(delt))_(avg)^t(Deltat)`

`(U_2)_(i)^(t+Deltat)=(U_2)_i^t+((delU_2)/(delt))_(avg)^t(Deltat)`

`(U_3)_(i)^(t+Deltat)=(U_3)_i^t+((delU_3)/(delt))_(avg)^t(Deltat)`

After this we need to decode the solution to obtain the required terms

Density:

`(rho')_(i)^(t+Deltat)=((U_1)_(i)^(t+Deltat))/(A')`

Velocity:

`(V')_(i)^(t+Deltat)=((U_1)_(i)^(t+Deltat))/((U_2)_(i)^(t+Deltat))`

Temperture:

`(T')_(i)^(t+Deltat)=(gamma-1)[((U_1)_(i)^(t+Deltat))/(U_2)_(i)^(t+Deltat)-gamma/2((V')_(i)^(t+Deltat))]`

Pressure:

`p'=rho'T'`

Boundary Condition and Initial Condition:

With referring to John D.Anderson book the boundary condition and initial condition are given

Non-Conservative form:

A.INITIAL CONDITION:

Nozzle shape for the present calculation is Non-dimensional distance along nozzle is 3m. The following condition describes the initial condition of the nozzle

The area is given by,

`A=1+2.2(x-1.5)^2`with the limits of `0<=x<=3`.

The choice of the initial condition is based on the Qualitatively behave in the same fashion. For simplicity, let us consider the flow-field variables, as a function of `x`. Let us assume the following values at the time `t=0`.

Density:

`rho=1-0.3146x`

Temperature:

`T=1-0.2314x`

Velocity:

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

B.BOUNDARY CONDITION:

INLET:

For the inlet values with respect to the reservoir are taken book is given by:

`rho(1)=1`

`T(1)=1`

`V(1)=2V(2)-V(3)`

OUTLET CONDITIONS:

The Non-dimensional variables are

`V_N=2V_(N-1)-V_(N-2)`

`rho_N=2rho_(N-1)-rho_(N-2)`

`T_N=2T_(N-1)-T_(N-2)`

2.Conservative Form:

A.INITIAL CONDITION:

The equations are solved by solution vectors, flux vectors, and source terms, therefore the initial condition must be given accordingly to the terms at the time `t=0`

`rho'=1, T'=1` for the limits of `0<=x'<=0.5`

`rho'=1-0.366(x'-0.5)` and `T'=0.1-0.167(x'-0.5)` for the limits of`0.5<=x'<=1.5`

`rho'=0.634-0.3879(x'-1.5)` and`T'=0.833-0.3507(x'-1.5)` for the limits of `1.5<=x'<=3.5`

Assuming the constant mass flow rate is given by

`V'=0.59/(rho'A')`

Using the equation above we get initial condition for `U_1, U_2, and U_3`

`U_1=rho'A'`

`U_2=rho'A'V'`

`U_3=rho'((e')/(gamma-1)+(gamma)/2V'^2)A'`

Where `e'=T'`. For the above initial condition `V'` is calculated such that `U_2=rho'A'V'=0.59`

B.BOUNDARY CONDITION:

The boundary condition for the CD-Nozzle are given by

A.INLET:

The inlet condition is

`(U_1)_(i=1)=(rho'A')_(i=1)`

The floating value of the `V'` in the inflow boundary is calculated at the end of each time step by linear extrapolating `U_2` from the known values at the internal grid points `i=2 and 3,`

`(U_2)_(i=1)=2(U_2)_(i-2)-(U_2)_(i=3)`

After Obtaining the `V'`

`U_3=U_1((T')/(gamma-1)+(gamma)/2V'^2)`

Time step calculation:

For both the conservative and non-conservative form in Quasi-dimensional hyperbolic equation are solved depends on the time step and CFL number is represented as

`Deltat=(lamdaDeltax)/(V+a)`

where,

`lamda=`CFL Number

`Deltax=` step size

`Deltat=`time step

`V=`Velocity

`a=` speed of sound

Mac-Cormack Method has stable criteria of the Courant number is should less than or equal to one for the Conservative and Non-conservative forms.

The Timestep for all nodes is calculated and the minimum time step is taken for the calculation(`dt=0.02304`).

Program:

1.Main program

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

%Actually plot
rho_act=0.634*ones(900,1400);
t_act = 0.833*ones(900,1400);
p_act = 0.528*ones(900,1400);
M_act =1*ones(900,1400);




tic
[rho,t_t,vel,rho_thor,vel_thor,t_t_thor,drhodt_Thor,dtdt_thor,pressure_t,mach_number1_t,dvdt_thor,non_massflow,throat,pressure,mach_number1]=Initial_program();
A=toc;
%%%%%%%%%%%%%%
R=rho;
T=t_t;
fprintf('Throat is %d',throat);

%%%%%%%%%%%%%%%
tic
[den,temp,V,error_den,error_t,error_v,v_throat,den_throat,T_throat,P_throat,M_throat,m_throat,nt,m,throat,p,M]=Conservative();
S=toc;



%plotting
figure(6)
subplot(3,1,1)
plot(x,R);
hold on 
plot(x,den,'*');
legend('NC','C');
ylabel('Density');
grid on
title('Comprasion between the properties NC and C');
subplot(3,1,2)
plot(x,T);
hold on
plot(x,temp,'o');
legend('NC','C');
ylabel('Temperature');
grid on
title('Comprasion between the properties NC and C');
subplot(3,1,3)
plot(x,vel);
hold on
plot(x,V,'+');
legend('NC','C');
ylabel('Velocity');
grid on
title('Comprasion between the properties NC and C');

%PlottingNEW
figure(30)
subplot(4,1,1)
plot(linspace(1,nt,nt),den_throat);
hold on
plot(linspace(1,nt,nt),rho_thor);
hold on
plot(linspace(1,nt,nt),rho_act);
grid on;
ylabel('Density');
legend('denC','denNC','actual');
%%%%
subplot(4,1,2)
plot(linspace(1,nt,nt),T_throat);
hold on
plot(linspace(1,nt,nt),t_t_thor);
hold on
plot(linspace(1,nt,nt),t_act);
grid on;
ylabel('Temperature');
legend('tempC','tempNC','actual');
%%%%
subplot(4,1,3)
plot(linspace(1,nt,nt),M_throat);
hold on
plot(linspace(1,nt,nt),mach_number1_t);
hold on
plot(linspace(1,nt,nt),M_act);
grid on;
ylabel('Mach number');
legend('machC','machNC','actual');
%%%
subplot(4,1,4)
plot(linspace(1,nt,nt),P_throat);
hold on
plot(linspace(1,nt,nt),pressure_t);
hold on
plot(linspace(1,nt,nt),p_act);
grid on;
ylabel('Pressure');
legend('pressureC','pressureNC','actual');
%%%%%%%
figure(39)
plot(x,non_massflow,'r');
hold on
plot(x,m,'b');
hold on
plot(x,linspace(0.579,0.579,n),'g');
legend('NC','C','actual');
title('Mass flow rate Comparsion Between the NC and C');
fprintf('Time taken for Non conservative %fs\n',A);
fprintf('Time taken for Conservative %fs\n',S);
fprintf('Density,Temperature,Velocity and Pressure is %fkg/m^3,%fK,%fm/s,and %fN/m^2 for NC =n is %d\n',rho(throat),t_t(throat),vel(throat),pressure(throat),n);
fprintf('Density,Temperature,Velocity and Pressure is %fkg/m^3,%fK,%fm/s,and %fN/m^2 for C =n is %d\n',den(throat),temp(throat),V(throat),p(throat),n);
fprintf('\nNC-Non Conservative,C-Conservative\n');
%PLotting
rhoact=[0.995 .994 .992 .990 .987 .983 .978 0.970 .958 .942 .920 .888 .844 .787 .716 .634 .547 .461 .382 .315 .258 .213 .176 .147 .124 .105 .090 .078 .068 .059 .052];
Mact=[0.098 0.110 .124 0.140 .160 .185 .214  .249 .293 .347 .413 .494 .592 .709 .845 1.0 1.169 1.348 1.531 1.715 1.896 2.071 2.240 2.402 2.557 2.706 2.848 2.983 3.114 3.239 3.359];   
figure(956)
yyaxis left
plot(x,rho);
hold on
plot(x,rhoact,'--','linewidth',1.25);
ylabel('Density');
yyaxis right
plot(x,mach_number1);
hold on
plot(x,Mact,'--','linewidth',1.25);
grid on;
legend('Numerical','Exact');
xlabel('Dimensionless distance');
ylabel('Mach number');
title('Numerical vs Exact for NC');
%%%%
figure(957)
yyaxis left
plot(x,den);
hold on
plot(x,rhoact,'--','linewidth',1.25);
ylabel('Density');
yyaxis right
plot(x,M);
hold on
plot(x,Mact,'--','linewidth',1.25);
grid on;
legend('Numerical','Exact');
xlabel('Dimensionless distance');
ylabel('Mach number');
title('Numerical vs Exact for NC');

2.Non-Conservative

function[rho,t_t,vel,rho_thor,vel_thor,t_t_thor,drhodt_Thor,dtdt_thor,pressure_t,mach_number1_t,dvdt_thor,non_massflow,throat,pressure,mach_number1]=Initial_program()

%Input 
n=31;
x= linspace(0,3,n);
dx=x(2)-x(1);
%tolerance=1e-8;


% Calculating the initial condition;

rho=1-(0.3146*x);
t_t=1-(0.2314*x);
vel=(0.1+1.09*x).*(t_t.^0.5);
%figure(12)
%plot(rho);
%figure(13)
%plot(t_t);
%figure(14)
%plot(vel);

c=0.5;
g=1.4;  %gamma

a=1+2.2*(x-1.5).^2;
throat=find(a==1);

%time step
nt=1400;
%courant's number
for i=1:n
del_t=(c*dx)./((sqrt(t_t)+vel));
end
dt=min(abs(del_t));


%time step

for k=1:nt
  rho_old=rho;
  vel_old=vel;
  t_t_old=t_t;
  
  
  
  
  %Preditor method
  for j=2:n-1
    %gnr equation
    dvdx=(vel(j+1)-vel(j))/dx;
    dadx=(log(a(j+1))-log(a(j)))/dx;
    drhodx=(rho(j+1)-rho(j))/dx;
    dtdx=(t_t(j+1)-t_t(j))/dx;
    
    %contiunity equation
    drhodt_P(j)=-(rho(j)*dvdx)-(rho(j)*vel(j)*(dadx))-(vel(j)*drhodx);
    %Momentum equation
    dvdt_p(j)=-(vel(j)*dvdx)-((1/g)*((dtdx)+((t_t(j)/rho(j))*(drhodx))));
    %Temperature equaton
    dtdt_p(j)=-(vel(j)*dtdx)-((g-1)*t_t(j)*((dvdx)+(vel(j)*dadx)));
    % solution update
    
    vel(j)=vel(j)+dvdt_p(j)*dt;
    rho(j)=rho(j)+drhodt_P(j)*dt;
    t_t(j)=t_t(j)+dtdt_p(j)*dt;
  end
  for j=2:n-1
    %gnr equation
    dvdx=(vel(j)-vel(j-1))/dx;
    dadx=(log(a(j))-log(a(j-1)))/dx;
    drhodx=(rho(j)-rho(j-1))/dx;
    dtdx=(t_t(j)-t_t(j-1))/dx;
    
    %contiunity equation
    drhodt_c(j)=-(rho(j)*dvdx)-(rho(j)*vel(j)*(dadx))-(vel(j)*drhodx);
    %Momentum equation
    dvdt_c(j)=-(vel(j)*dvdx)-((1/g)*((dtdx)+((t_t(j)/rho(j))*(drhodx))));
    %Temperature equaton
    dtdt_c(j)=-(vel(j)*dtdx)-((g-1)*t_t(j)*[(dvdx)+(vel(j)*dadx)]);

    
  end
  %average
  drhodt=0.5*(drhodt_P+drhodt_c);
  dvdt=0.5*(dvdt_p+dvdt_c);
  dtdt=0.5*(dtdt_p+dtdt_c);
  %Final update
 for m=2:n-1
   rho(m)=rho_old(m)+(drhodt(m)*dt);
   vel(m)=vel_old(m)+(dvdt(m)*dt);
   t_t(m)=t_t_old(m)+(dtdt(m)*dt);
  
 end


    %inlet
    vel(1) = (2*vel(2)) - vel(3);
    rho(1) = 1;
    t_t(1)=1;
    
    %BC at outlet
    vel(n) = 2*vel(n-1) - vel(n-2);
    rho(n) = 2*rho(n-1) - rho(n-2);
    t_t(n) = 2*t_t(n-1) - t_t(n-2);
    
    %Mass Flow Rate
    inflow = rho(1)*vel(1)*a(1);
    outflow = rho(n)*vel(n)*a(n);
    massflow(k) = (outflow-inflow);
    
    non_massflow=rho.*vel.*a;
hold on
figure(1)
if k == 50
        plot(x,non_massflow,'linewidth',1.75);
        %legend('50');
    else if k == 100
            plot(x,non_massflow,'g','linewidth',1.75);
            %legend('100');
    else if k == 200
         plot(x,non_massflow,'b','linewidth',1.75);
         %legend('200');
    else if k == 400
       plot(x,non_massflow,'r','linewidth',1.75);
       %legend('400');
    else if k == 600
     plot(x,non_massflow,'y','linewidth',1.75);
     %legend('600');
   else if k==1400
     plot(x,non_massflow,'k','linewidth',1.75);
     %legend('1400');
   end
 end
end
end
end
end
title('Dimensionless distance vs Mass flow rate');
xlabel('x/L :Non-dimensional distance through the nozzle');
ylabel('Non-Dimensional Mass flow rate');
%legend('50','100','200','400','600','1400');
grid on
    
    
    %Non-dimensional Pressure calculation pressure=density*temperature
    
    pressure=rho.*t_t;
    %Mach number calculation mach_number=velocity/sqrt(temperature)
    mach_number1=vel./sqrt(t_t);
    
    drhodt_Thor(k)=drhodt(throat);
    dvdt_thor(k)=dvdt(throat);
    dtdt_thor(k)=dtdt(throat);
    rho_thor(k)=rho(throat);
    vel_thor(k)=vel(throat);
    t_t_thor(k)=t_t(throat);
    mach_number1_t(k)=mach_number1(throat);
    pressure_t(k)=pressure(throat);
    
    %Error Tolerance
    
    error_v(k) = max(vel - vel_old);
    error_rho(k) = max(rho - rho_old);
    error_t(k) = max(t_t - t_t_old);
    figure(3)
    subplot(3,1,1)
    plot(error_v);
    axis([0 inf -0.005 0.04]);
    xlabel('Time step');
    ylabel('Error on velocity');
    title('NC:Time step VS Error on velocity');
    grid on;
    
    subplot(3,1,2)
    plot(error_rho);
    axis([0 inf -0.005 0.04]);
    xlabel('Time step');
    ylabel('Error on Density');
    grid on;
    title('NC:Time step VS Error on Density');
    subplot(3,1,3)
    plot(error_t);
    axis([0 inf -0.005 0.02]);
    xlabel('Time step');
    ylabel('Error on temperature');
    grid on;
    title('NC:Time step VS Error on temperature');   
    
    %if(error_v(k) < tolerance && error_rho(k) < tolerance && error_t(k) < tolerance&&massflow(k)<tolerance)
       %break;
    %end
    
end

figure(2)
subplot(4,1,1)
plot(x,mach_number1)
title('Mach Number vs Distance')
xlabel('Length of Nozzle')
ylabel('Mach Number')
grid on
subplot(4,1,2)
plot(x,pressure)
title('Pressure Ratio vs Distance')
xlabel('Length of Nozzle')
ylabel('Pressure Ratio')
grid on
subplot(4,1,3)
plot(x,rho)
title('Density ratio vs Distance')
xlabel('Length of Nozzle')
ylabel('Density ratio')
grid on
subplot(4,1,4)
plot(x,t_t)
title('Temperature vs Distance')
xlabel('Length of Nozzle')
ylabel('Temperature')
grid on
figure(0899)
    plot(drhodt_Thor);
    hold on;
    plot(dvdt_thor);
    hold on
    plot(dtdt_thor);
    grid on;
    legend('drhodt Thor','dvdt thor','dtdt thor');
    title('Residual Plot');
figure(1)
legend('50','100','200','400','600','1400');
end

3.Conservative

function[den,temp,V,error_den,error_t,error_v,v_throat,den_throat,T_throat,P_throat,M_throat,m_throat,nt,m,throat,p,M]=Conservative()


%Conservative for the better solution
n=31;
x=linspace(0,3,n);
dx=x(2)-x(1);

%tolerance=1e-4;

A=1+2.2*(x-1.5).^2;  %area

throat=find(A==1);

C=0.5;

gamma=1.4;

%Calculating the initial conditions
for i=1:n
    if(x(i)>=0&&x(i)<0.5)
      den(i)=1;
      temp(i)=1;
    else
      if(x(i)>=0.5&&x(i)<1.5)
      den(i)=1-(0.366*(x(i)-0.5));
      temp(i)=1-(0.1678*(x(i)-0.5));
        else
          if(x(i)>=1.5&&x(i)<3.5)
          den(i)=0.634-(0.3879*(x(i)-1.5));
          temp(i)=0.833-(0.3507*(x(i)-1.5));
        end
      end
    end
end
%Calculating Velocity

  V=0.59./(den.*A);
  U1=(den.*A);
  U2=den.*A.*V;
  %v=U2.*(den.A)^-1;
  U3=U1.*((temp/(gamma-1))+((gamma/2).*(V.^2)));
 

%dt=min(C.*(dx./(sqrt(temp)+v)));

nt=1400;

for k=1:nt

  %INitial
  U1_initial=U1;
  U2_initial=U2;
  U3_initial=U3;
  dt = min(0.5.*dx./(sqrt(temp) + V));
  %Arrange the values 
  den_old=U1_initial./A;
  V_old=U2_initial./U1_initial;
  temp_old=(gamma-1).*((U3_initial./U1_initial)-((gamma).*((V_old.^2)/2)));
  p_old=den_old.*temp_old;

  %force terms
  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));
  %predictor step
  for i=2:n-1
    dF1dx_p=(F1(i+1)-F1(i))/dx;
    dF2dx_p=(F2(i+1)-F2(i))/dx;
    dF3dx_p=(F3(i+1)-F3(i))/dx;
    dadx=(A(i+1)-A(i))/dx;
    J2(i)=(1/gamma).*den(i).*temp(i).*dadx;
     %contiunity equation
    dU1dt_p(i)=-dF1dx_p;
    
    %momentum equation
    dU2dt_p(i)=-dF2dx_p+J2(i);
    
    %Energy eqaution
    dU3dt_p(i)=-dF3dx_p;
    
    %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
  %Arrange the values 
  den=U1./A;
  V=U2./U1;
  temp=(gamma-1).*((U3./U1)-((gamma).*((V.^2)/2)));
  p=den.*temp;
  %F
  F1=U2;
  F2=((U2.^2)./(U1))+((gamma-1)/gamma).*(U3-(0.5*gamma.*U2.^2)./U1);
  F3=((gamma.*U2.*U3)./U1)-(gamma*0.5*(gamma-1).*(U2.^3)./((U1.^2)));
  %Collector step
  for j=2:n-1
    dF1dx_c=(F1(j)-F1(j-1))/dx;
    dF2dx_c=(F2(j)-F2(j-1))/dx;
    dF3dx_c=(F3(j)-F3(j-1))/dx;
    dadx_c=(A(j)-A(j-1))/dx;
    J3(j)=(1/gamma).*den(j).*temp(j).*dadx_c;
     %contiunity equation
    dU1dt_c(j)=-dF1dx_c;
    
    %momentum equation
    dU2dt_c(j)=-dF2dx_c+J3(j);
    
    %Energy eqaution
    dU3dt_c(j)=-dF3dx_c;
   end
   
   %Average
   dU1dt_avg=0.5*(dU1dt_p+dU1dt_c);
   dU2dt_avg=0.5*(dU2dt_p+dU2dt_c);
   dU3dt_avg=0.5*(dU3dt_p+dU3dt_c);
   
   % final solution update
      
      for m = 2:n-1
          U1(m) = U1_initial(m) +dU1dt_avg(m)*dt;
          U2(m) = U2_initial(m) +dU2dt_avg(m)*dt;
          U3(m) = U3_initial(m) +dU3dt_avg(m)*dt;
      end
   
   % apply boundary conditions
      %Inlet
      U2(1) = (2*U2(2)) - U2(3);
      
      %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);
      
      %update
      den = U1./A;
      V = U2./U1;
      temp= (gamma-1).*((U3./U1)-((gamma*(V.^2))/2));
      p = den.*temp;
      
      %mach number and mass flow rate
      M = V./((temp).^0.5);
      m = den.*V.*A;% value sat throat
      den_throat(k) = den(throat);
      v_throat(k) = V(throat);
      T_throat(k) = temp(throat);
      P_throat(k) = p(throat);
      M_throat(k) = M(throat);
      m_throat(k) = m(throat);
      dU1dt_avgthroat(k)=dU1dt_avg(throat);
      dU2dt_avgthroat(k)=dU2dt_avg(throat);
      dU3dt_avgthroat(k)=dU3dt_avg(throat);
      %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      figure(41)
      %plotting
      hold on
      if k==50
        plot(x,m,'r','linewidth',1.75);
        %legend('50');
      else
        if k==150
          plot(x,m,'b','linewidth',1.75);
          %legend('150');
        else
          if k==400
            plot(x,m,'g','linewidth',1.75);
            %legend('400');
          else 
            if k==800
              plot(x,m,'k','linewidth',1.75);
              %legend('800');
            else
              if k==1400
                plot(x,m,'y','linewidth',1.75);
                %legend('1400');
                %title('Dimensionless distance vs mass flow rate');
                %legend('50','150','400','800','1400');
                grid on;
               end
             end
           end
         end
       end
       title('Dimensionless distance vs mass flow rate');
       xlabel('x/L :Non-dimensional distance through the nozzle');
        ylabel('Non-Dimensional Mass flow rate');
       %legend('50','150','400','800','1400');
    error_v(k) = max(V - V_old);
    error_den(k) = max(den - den_old);
    error_t(k) = max(temp - temp_old);
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    figure(31)
    subplot(3,1,1);
    plot(error_v);
    axis([0 inf -0.005 0.04]);
    grid on;
    xlabel('Time step');
    ylabel('Error on velocity');
    title('Time step VS Error on velocity');
    subplot(3,1,2)
    %figure(15)
    plot(error_den);
    axis([0 inf -0.005 0.01]);
    grid on;
    xlabel('Time step');
    ylabel('Error on Density');
    title('Time step VS Error on Density');
    subplot(3,1,3)
    %figure(56)
    plot(error_t);
    axis([-0.5 inf -0.05 0.05]);
    grid on; 
    xlabel('Time step');
    ylabel('Error on temperature');
    title('Time step VS Error on Temperature');
    %if(error_v(k) < tolerance && error_den(k) < tolerance && error_t(k) < tolerance)
       %break;
    %end
end
figure(21)
subplot(4,1,1)
plot(x,M)
title('Mach Number vs Distance')
xlabel('Length of Nozzle')
ylabel('Mach Number')
grid on
subplot(4,1,2)
plot(x,p)
title('Pressure Ratio vs Distance')
xlabel('Length of Nozzle')
ylabel('Pressure Ratio')
grid on
subplot(4,1,3)
plot(x,den)
title('Density ratio vs Distance')
xlabel('Length of Nozzle')
ylabel('Density ratio')
grid on
subplot(4,1,4)
plot(x,temp)
title('Temperature vs Distance')
xlabel('Length of Nozzle')
ylabel('Temperature')
grid on
figure(41)
legend('50','150','400','800','1400');
figure(089)
    plot(dU1dt_avgthroat,'g');
    hold on;
    plot(dU2dt_avgthroat,'b');
    hold on
    plot(dU3dt_avgthroat,'r');
    grid on;
    legend('dU1dt avgthroat','dU2dt avgthroat','dU3dt avgthroat');
    title('Residual Plot');
end

RESULT AND DISCUSSION:

1.Non-Conservative:

Figure 1a: Error variation according to the time-step variation

Table 1a:

Description	Exact Solution at the throat	Numerical Solution at the throat
Density	0.634	0.638664
Temperature	0.833	0.836451
Pressure	0.528	0.534211
Mach number	1	0.9994

Figure 1a: This illustrates the error variance of the variables(Density, Temperature, Velocity). It is Observed as the initial time step has a large variation between the old step and new step and it is approaching the steady-state variation slowly decrease to the zero.Large Variation in first is due to the solution is transient and approaching the steady-state with sequence steps.

In the above plot, the error variation approaches to zero, and then the solution is converged. It seems around the 400th time step error approaches to zero, but to obtain more accuracy the solution is running towards the 1400 time steps.

In the table, the variables are very close to the exact variables. Therefore, the solution is said to be converged. Also, it has an error of less than 2%, which is acceptable.

b.Mass Flow variation through a non-dimensional distance of the nozzle for the NON-CONSERVATIVE approach

The above plot shows the variation in the mass flow for several time steps. It is observed the initial time steps are large in variation due to the product of the velocity and density. In the 50th time step, the drastic changes occur in the mass flow. In the 100th time step, the solution is approaching the time-marching.In the 200th time step, the solution approaches a constant state. It is observed, the solution is getting constant over the domain after the 400th time. The plot shows the solution is converged after the 400th time step.

C.Properties of NON-CONSERVATIVE:

In the nozzle, the inlet values mostly resemble from the reservoir. At the converging section, the flow is subsonic, the Mach number is less than one. The velocity tends to increase and the pressure, temperature, density, and area tends to decrease. In the throat, the Mach number equals the one and the flow becomes sonic. At the diverging section, the flow is supersonic, the Mach number tends to increase. The velocity and area tend to increase and the pressure, temperature, density, tend to further decrease.

D.Residual plot:

This plot shows the error (or) residual decrease with the increase in the time step. It seems like the initial plot has a large oscillation due to the solution is time marching. Once the solution approaches the steady-state the values of the `((drho)/dt)_(throat)`,`((dV)/dt)_(throat)`,`((dT)/dt)_(throat)`tends to zero. As the truncation error tends to zero, the Mac-Cormack technique is stable and gives the consistency for the given CFL Constant`lambda<=1`and the solution converges.

Comparison with the EXACT solution:

The exact solution is the condition, we assumed the velocity at the inlet tends to zero but in the actual case, there is some velocity in the flow. The exact solution is obtained from the JOHN ANDERSON Book

Above the table, the exact analytical results are taken for reference. The Numerical value obtained is given below

>> rho

rho =

  Columns 1 through 18

    1.0000    0.9978    0.9970    0.9943    0.9916    0.9873    0.9818    0.9738    0.9626    0.9467    0.9241    0.8925    0.8492    0.7921    0.7210    0.6387    0.5510    0.4650

  Columns 19 through 31

    0.3863    0.3184    0.2618    0.2159    0.1791    0.1496    0.1261    0.1072    0.0919    0.0795    0.0690    0.0610    0.0529

>> rhoact

rhoact =

  Columns 1 through 18

    0.9950    0.9940    0.9920    0.9900    0.9870    0.9830    0.9780    0.9700    0.9580    0.9420    0.9200    0.8880    0.8440    0.7870    0.7160    0.6340    0.5470    0.4610

  Columns 19 through 31

    0.3820    0.3150    0.2580    0.2130    0.1760    0.1470    0.1240    0.1050    0.0900    0.0780    0.0680    0.0590    0.0520

>> mach_number1

mach_number1 =

  Columns 1 through 10

    0.0988    0.1121    0.1253    0.1428    0.1625    0.1870    0.2162    0.2521    0.2958    0.3496

  Columns 11 through 20

    0.4155    0.4962    0.5938    0.7103    0.8459    0.9994    1.1674    1.3451    1.5276    1.7103

  Columns 21 through 30

    1.8901    2.0646    2.2330    2.3942    2.5490    2.6963    2.8386    2.9724    3.1053    3.2249

  Column 31

    3.3532

>> Mact

Mact =

  Columns 1 through 10

    0.0980    0.1100    0.1240    0.1400    0.1600    0.1850    0.2140    0.2490    0.2930    0.3470

  Columns 11 through 20

    0.4130    0.4940    0.5920    0.7090    0.8450    1.0000    1.1690    1.3480    1.5310    1.7150

  Columns 21 through 30

    1.8960    2.0710    2.2400    2.4020    2.5570    2.7060    2.8480    2.9830    3.1140    3.2390

  Column 31

    3.3590

The plot shows the small error between the exact and numerical solutions. The error, in this case, is acceptable. the error may be caused by assumption, Rounding-off values, and other factors.

2.CONSERVATIVE:

a.ERROR VARIATION:

Figure(2a):

This illustrates the error variance of the variables(Density, Temperature, Velocity). It is Observed as the initial time step has a large variation between the previous step and current step and it is approaching the steady-state variation slowly decrease to the zero. Large Variation in first is due to the solution is transient and approaching the steady-state with sequence steps.

In the above plot, the error variation approaches to zero, and then the solution is converged. It seems around the 300th time step error approaches to zero, but to obtain more accuracy the solution is running towards the 1400 time steps.

Table 2a

Description	Exact Solution at the throat	Numerical Solution at the throat
Density	0.634	0.650459
Temperature	0.833	0.840082
Pressure	0.528	0.546439
Mach number	1	0.9816

In the table, the variables are very close to the exact variables. Therefore, the solution is said to be converged. Also, it has an error of less than 3%, which is acceptable.

b.Mass Flow variation through a non-dimensional distance of the nozzle for the CONSERVATIVE approach

The above plot shows the variation in the mass flow for several time steps. It is observed the initial time steps are large in variation due to the product of the velocity and density. In the 50th time step, the drastic changes occur in the mass flow. In the 150th time step, the solution is approaching the time-marching.In the 400th time step, the solution approaches a steady state. It is observed, the solution is getting constant over the domain after the 400th time. The plot shows the solution is converged before the 400th time step.

C.Properties of CONSERVATIVE:

In the nozzle, the inlet values(`P_o,T_o,rho_o`) mostly resemble from the reservoir. At the converging section, the flow is subsonic, the Mach number is less than one`(M<1)`. The velocity tends to increase and the pressure, temperature, density, and area tends to decrease. In the throat, the Mach number equals the one`(M=1)` and the flow becomes sonic. At the diverging section, the flow is supersonic, the Mach number tends to increase`(M>1)`. The velocity and area tend to increase and the pressure, temperature, density, tend to further decrease.

D.Residual plot:

This plot shows the error (or) residual decrease with the increase in the time step. It seems like the initial plot has a large oscillation due to the solution is time marching. Once the solution approaches the steady-state the values of the `((dU1)/(dt))_(throat),((dU2)/(dt))_(throat),((dU3)/(dt))_(throat)`tends to zero. As the truncation error tends to zero, the Mac-Cormack technique is stable and gives the consistency for the given CFL Constant`lambda<=1`and the solution converges.

Comparison with the EXACT solution:

The exact solution is the condition, we assumed the velocity at the inlet tends to zero but in the actual case, there is some velocity in the flow. The exact solution is obtained from the JOHN ANDERSON Book.

>> den

den =

  Columns 1 through 18

    1.0000    0.9987    0.9971    0.9949    0.9918    0.9877    0.9819    0.9738    0.9623    0.9460    0.9229    0.8905    0.8464    0.7886    0.7173    0.6505    0.5473    0.4619

  Columns 19 through 31

    0.3837    0.3160    0.2597    0.2139    0.1772    0.1480    0.1246    0.1058    0.0906    0.0783    0.0681    0.0598    0.0527

>> rhoact

rhoact =

  Columns 1 through 18

    0.9950    0.9940    0.9920    0.9900    0.9870    0.9830    0.9780    0.9700    0.9580    0.9420    0.9200    0.8880    0.8440    0.7870    0.7160    0.6340    0.5470    0.4610

  Columns 19 through 31

    0.3820    0.3150    0.2580    0.2130    0.1760    0.1470    0.1240    0.1050    0.0900    0.0780    0.0680    0.0590    0.0520

>> M

M =

  Columns 1 through 18

    0.0983    0.1103    0.1245    0.1413    0.1614    0.1856    0.2150    0.2509    0.2950    0.3492    0.4159    0.4977    0.5971    0.7158    0.8540    0.9816    1.1821    1.3613

  Columns 19 through 31

    1.5448    1.7281    1.9081    2.0828    2.2511    2.4125    2.5673    2.7149    2.8571    2.9914    3.1241    3.2433    3.3760

>> Mact

Mact =

  Columns 1 through 18

    0.0980    0.1100    0.1240    0.1400    0.1600    0.1850    0.2140    0.2490    0.2930    0.3470    0.4130    0.4940    0.5920    0.7090    0.8450    1.0000    1.1690    1.3480

  Columns 19 through 31

    1.5310    1.7150    1.8960    2.0710    2.2400    2.4020    2.5570    2.7060    2.8480    2.9830    3.1140    3.2390    3.3590

The FOLLOWING plot shows the small error between the exact and numerical solutions. The error, in this case, is acceptable. the error may be caused by assumption, Rounding-off values, and other factors.

Comparison between the two methods of Mac-Cormack method:

A.MASS FLOW:

The above plot shows the comparison between the Mass flow rate and the dimensionalize distance of Nozzle. The comparison gives the advantage of the conservative form over the non-conservative form. The mass flow is the dependent variable in the conservative approach whereas the mass flow rate in the non-conservative form secondary dependent variable.

B.Properties obtained from the conservative and non-conservative form

This plot shows the properties of the flow after the 1400th time step in Mac-Cormack Method.

This plot shows the properties comparison between the Conservative, Non-Conservative, and Exact solution at the throat section.

C.Grid independent test:

Description	`(rho^**)/(rho_o)`	`(T^**)/(T_o)`	`(p^**)/(p_o)`	`(M^**)/(M_o)`
Exact Solution	0.634	0.833	0.528	1
Non-conservative form(n=31)	0.638664	0.836451	0.534211	0.9994
Conservative form(n=31)	0.650459	0.840082	0.546439	0.9816
Non-Conservative form(n=61)	0.637659	0.835439	0.532726	0.9998
Conservative form(n=61)	0.638355	0.835245	0.533182	0.9991

The table above gives the effect of an increase in the number of grids by comparing the steady-state (analytical solution) results with the conservative and non-conservative form at the throat. It seems like the number of grid increase the accuracy of the numerical solution in both forms. The Non-Conservative approach has better accuracy than the conservative form. Due to practical difficulties and processing time, we use the less number of iteration like 31 nodes.

CONCLUSION:

• When the Mass flow rate is compared between the two forms conservative and non-conservative.Conservative form is mainly based on the Principle of mass conservation. It is observed the conservative form conserve the mass flow rate more than the non-conservative one.
• From the conservative and non-conservative approach, it is observed that the Non-Conservative form has more accurate results than the conservative form. The Error percentage in the Non-Conservative form is less than 2%, whereas the conservative approach percentage is more likely 3%.
• The error or residual in the non-conservative approaches decreases. The convergence state of the solution is obtained in less time because the non-conservative is primarily uses the primitive variables for the solution.
• Whereas the conservative form uses the solution vector, flux vector to obtain the convergence state of the solution. Processing time increases in converting into the primitive variables. The error or residual increase because of this cause.