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

# The aim of the project is to simulate the flow through a supersonic (1D) nozzle using MacCormack's method:

* Supersonic nozzle: It is a convergent and divergent nozzle with a pinched section in the middle. It is mainly used in steam turbines and jet engines to accelarate a high pressurised gas to supersonic speed.

Below is the diagramatic representation of the supersonic nozzle:

As shown in the above diagram it has a convergent section, a throat section and a divergent section, the high pressurised gas enters from convergent section and comes out of diergent section at supersonic speed.

Mach number is a dimensionless number used to estimate the speed at different cross sections of th nozzle.

Mach number, `M = u/c`, where u - Local velocity of flow (m/s)

c - speed of sound in the medium (m/s).

at convergent section the mach number is less than 1, M < 1. i.e the local velocity of the flow is less than the speed of sound in that particular medium

at throat section the mach number , M = 1 i.e, the flow velocity is equal to speed of sound in the medium.

at divergent section mach number M > 1, i.e, the flow velocity is greater than the speed of sound in the medium which means the output velocity of the flow has attained the supersonic speed.

Here both the non conservative form as well as the conservative form of the navier stokes equation are solved and the results are compared.

# Non conservative form:

The continuity, momentum and energy equation for 1d nozzle in its non conservative form is as given below:

* Continuity equation:

`(delrho)/(delt) = -rho*(delv)/(delx)-(rho*v*(dellog(a))/(del(x)))-v*(delrho)/(delx)`

* Momentum equation:

`(delv)/(delt) = -v*(delv)/(delt)-(1/gamma*((delT)/(delx)+(T/rho*(delrho)/(delx))))`

* Energy equation:

`(delT)/(delt) = -v*(delT)/(delx)-((gamma-1)*T*((delv)/(delx)+(v*(del(log(a)))/(delx))))`

where rho - density

T - Temperature

v - Velocity

a - area

all the above variables are dimensionless variables.

Above set of equations are solved numerically using MacCormack's method,

Below is the representation of the Method of solving using MacCormack's method:

It basically involves 6 step:

1. Setting the initial conditions:

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)

2. Predictor method:

continuity equation:

`(delrho)/(delt)_p = -rho(i)*((v(i+1)-v(i))/dx)-(rho(i)*v(i)*((log(a(i+1))-log(a(i)))/dx))-v(i)*((rho(i+1)-rho(i))/dx)`

momentum equation:

`(delv)/(delt)_p = -v(i)*((v(i+1)-v(i))/dx)-(1/gamma*(((T(i+1)-T(i))/dx))+((T(i))/(rho(i)))*((rho(i+1)-rho(i))/dx)`

Energy equation:

`(delT)/(delt)_p = -v(i)*((T(i+1)-T(i))/dx)-((gamma-1)*T(i)*(((v(i+1)-v(i))/dx)+(v(i)*((log(a(i+1))-log(a(i)))/dx))))`

In predictor mehod as shown above the partial differentials on the right hand side is forward differenced and hence the partial diferential terms which are on the left side of the equation are computed for all the grid points

3. Updating the values of variables:

After the predictor step is completed the next step is to update the value of all the variables at all the grid points

`rho(i) = rhoold(i)+((delrho)/(delt)_p*dt)`

`v(i) = vold(i) + ((delv)/(delt)_p*dt)`

`T(i) = Told(i) +((delT)/(delt)*dt)`

4. Corrector step:

After updating the values of all the variables, then the updated values of the variables are used to compute the partial differentials in the left hand side using the corrector method

continuity equation:

`(delrho)/(delt)_c = -rho(i)*((v(i)-v(i-1))/dx)-(rho(i)*v(i)*((log(a(i))-log(a(i-1)))/dx))-v(i)*((rho(i)-rho(i-1))/dx)`

momentum equation:

`(delv)/(delt)_c = -v(i)*((v(i)-v(i-1))/dx)-(1/gamma*(((T(i)-T(i-1))/dx))+((T(i))/(rho(i)))*((rho(i)-rho(i-1))/dx)`

Energy equation:

`(delT)/(delt)_c = -v(i)*((T(i)-T(i-1))/dx)-((gamma-1)*T(i)*(((v(i)-v(i-1))/dx)+(v(i)*((log(a(i))-log(a(i-1)))/dx))))`

In this method as shown above the partial differentials on the right hand side of the equation is backward differenced, hence the partial differential of variable which is present in the left side of the equation is computed for all the gird points

5. Computing the average:

Using the partial differential of variable that is obtained from predictor method as well as corrector method is used to calculate the average of the same:

`((delrho)/(delt))_(avg) = 0.5* (((delrho)/(delt))_(p) + ((delrho)/(delt))_(c))`

`((delv)/(delt))_(avg) = 0.5*(((delv)/(delt))_(p)+((delv)/(delt))_(c))`

`((delT)/(delt))_(avg) = 0.5* (((delT)/(delt))_(p)+((delT)/(delt))_(c))`

6. Final update:

The average value that is computed above is used to update the values of all the variables at all the grid points:

`rho(i) = rhoold(i)+(((delrho)/(delt))_(avg)*dt)`

`v(i) = vold(i) + (((delv)/(delt))_(avg)*dt)`

`T(i) = Told(i) +(((delT)/(delt))_(avg)*dt)`

The above mentioned steps are repeated for all each and every time step. In this simulation we will be simulating the above variables for 1400 time steps, for a total number of grid points n = 31.

After each time step the boundary conditions are updated as follow:

Inlet:

v(1) = 2*v(2)-v(3)

outlet:

v(n) = (2*v(n-1))-v(n-2)

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

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

# Conservative form:

In conservative form the equations are not expressed directly in terms of primitive variables, but in terms of solution vector's, flux vector's and source term.

Navier stokes equation in conservative form is given as: (In terms of nozzle flow)

* Continuity equation:

`(del(rho*a))/(delt) +(del(rho*a*v))/(delx) = 0`

* Momentum equation:

`(del(rho*a*v))/(delt)+((del(rho*v^2*a+(1/gamma)*p*a))/(delt)) = p*(1/gamma)*(dela)/(delx)`

* Energy equation:

`(del((rho*a*(e/(gamma-1))) + ((gamma/2)*v^2)))/(delt) + (del((rho*a*v*(e/(gamma-1)+(gamma/2)*v^2))+(p*a*v)))/(delx) = 0`

Above equations are expressed in terms of solution vectors, flux vectors and source term

% Solution vectors are given as follow:

`U1 = rho*a`

`U2 = rho*a*v`

`U3 = U1*((T/(gamma-1))+((gamma/2)*v^2))`

% Flux terms are given as follow:

`F1 = U2`

`F2 = ((U2^2)/(U1))+(((gamma-1)/gamma)*(U3-((gamma*0.5*(U2^2))/(U1)))`

`F3 = ((gamma*U3*U2)/(U1))-((gamma*0.5*(gamma-1)*(U2^3))/(U1^2))`

% Source term is given as follow:

`J2 = (1/gamma)*rho*T*(dela)/(delx)`

so the above equations of conservative form which are in terms of primitive variables can be expressed in terms of solution vector, flux vector and source term.

* Coninuity equation:

`(delU1)/(delt) = -(delF1)/(delx)`

* Momentum equation:

`(delU2)/(delt) = -(delF2)/(delx) + J2`

* Energy equation:

`(delU3)/(delt) = -(delF3)/(delx)`

Now these above equations are used in MacCormack's method,

1. Setting the initial conditions:

The initial condition set for conservative form is bit different from the one which is set for the non conserative form except the area, Initial conditions used are as below:

For rho and T:

for 0 <= x <= 0.5

rho = 1.0

T = 1.0

for 0.5 <= x <= 1.5

rho = 1.0-0.366*(x-0.5)

T = 1.0-0.167*(x-0.5)

for 1.5 <= x <= 3.5

rho 0.634-0.3879*(x-1.5)

T = 0.833-0.3507*(x-1.5)

2. Predictor method:

continuity equation:

`((delU1)/(delt))_(p) = -1*((F1(i+1)-F1(i))/dx)`

momentum equation:

`((delU2)/(delt))_(p) = -1*((F2(i+1)-F2(i))/dx)+J2`

Energy equation:

`((delU3)/(delt))_(p) = -1*((F3(i+1)-F3(i))/dx)`

In predictor mehod as shown above the partial differentials on the right hand side is forward differenced and hence the partial diferential terms which are on the left side of the equation are computed for all the grid points

3. Updating the values of variables:

After the predictor step is completed the next step is to update the value of all the solution vectors at all the grid points

`U1(i) = U1_(old)(i)+(((delU1)/(delt))_(p)(i)*dt)`

`U2(i) = U2_(old)(i) + (((delU2)/(delt))_(p)(i)*dt)`

`U3(i) = U3_(old)(i) +(((delU3)/(delt))_(P)(i)*dt)`

4. Corrector step:

After updating the values of all the solution vectors, these updated values are in turn used to update the values of flux terms and these updated flux terms are used in corrector method.

continuity equation:

`((delU1)/(delt))_(c) = -1*((F1(i)-F1(i-1))/dx)`

momentum equation:

`((delU2)/(delt))_(c) = -1*((F2(i)-F2(i-1))/dx)+J2`

Energy equation:

`((delU3)/(delt))_(c) = -1*((F3(i)-F3(i-1))/dx)`

In this the partial differenials on the right hand side is backward differenced and hence the partial differentials of all the solution vectors on the right hand side is calculated at all grid points.

5 . Computing the average:

Using the partial differentials of solution vector obtained from predictor method and corrector method the average of all the solution vector is calculated:

`((delU1)/(delt))_(avg) = 0.5*(((delU1)/(delt))_(p)+((delU1)/(delt))_(c))`

`((delU2)/(delt))_(avg) = 0.5*(((delU2)/(delt))_(p)+((delU2)/(delt))_(c))`

`((delU3)/(delt))_(avg) = 0.5*(((delU3)/(delt))_(p)+((delU3)/(delt))_(c))`

6. Final update:

The average value that is computed above is used to update the values of all the solution vectors at all the grid points:

`U1(i) = U1_(old)(i)+(((delU1)/(delt))_(avg)(i)*dt)`

`U2(i) = U2_(old)(i) + (((delU2)/(delt))_(avg)(i)*dt)`

`U3(i) = U3_(old)(i) +(((delU3)/(delt))_(avg)(i)*dt)`

These steps are repeated for all the time steps, after every time step the boundary conditions are updated as follow,

U1(1) = 2*(U1(n-1))-U1(n-2)

U2(1) = (2*U2(n-1))-U2(n-2)

U3(1) = (2*U2(n-1))-U3(n-2)

and also the values of primitive variables are extracted from solution vector as follow:

`rho = (U1)/a`
`v = (U2)/(U1)`
`T = (gamma-1)*((U3)/(U1) - ((gamma/2)*v^2))`

Time step calculation:

dt = min((C*dx)/(sqrt(T)+v)), time step s calculated using this equation, where in the courant number C = 0.5 is taken

* Below are the codes:

# Main function:

clear all
close all
clc

% Defining the No. of grid points;
n=31;

% Defining the nozzle domain;
x = linspace(0,3,n);

dx = x(2)-x(1);

% Defining the length of the nozzle domain;
L = 3;

% Defining the No. of time steps;
nt = 1400;

% Solving using Non conservative form;
[rho_nc, T_nc, v_nc, rho_tr_nc, T_tr_nc, pr_tr_nc, mach_tr_nc, mass_flow_nc, time_nc] = nonconserv_n(x,n,nt,dx);

% Solving using conservative form;
[rho_c, T_c, v_c, rho_tr_c, T_tr_c, pr_tr_c, mach_tr_c, mass_flow_c, time_c] = conserv_n(x,n,L,nt,dx);

% Defining the exact value of variables;
rho_exact = 0.634*ones(1,1400);

T_exact = 0.833*ones(1,1400);

p_exact = 0.528*ones(1,1400);

Mach_exact = 1*ones(1,1400);

% Ploting for the non conservative form;

% Ploting for variables at throat;
figure(1)
subplot(4,1,1)
plot(T_tr_nc,'linewidth',2);
hold on;
plot(T_exact,'linewidth',2,'g');
xlabel('time');
ylabel('T/T_0');
title({"Plot of variables at throat for Non conservative form";"Plot of temperature at throat"});
legend("Time wise variation","exact solution");

subplot(4,1,2);
plot(rho_tr_nc,'linewidth',2);
hold on;
plot(rho_exact,'linewidth',2,'g');
xlabel('time');
ylabel('rho/rho_0');
title('Plot of density values at throat');
legend("Time wise variation","exact solution");

subplot(4,1,3);
plot(pr_tr_nc,'linewidth',2);
hold on;
plot(p_exact,'linewidth',2,'g');
xlabel('time');
ylabel('p/p_0');
title('Plot of pressure values at throat');
legend("Time wise variation", "exact solution");

subplot(4,1,4);
plot(mach_tr_nc,'linewidth',2);
hold on;
plot(Mach_exact,'linewidth',2,'g');
xlabel('time');
ylabel('M/M_0');
title('Plot of mach number at throat');
legend("Time wise variation","exact solution");

figure(2);
% Ploting for variables at nozzle doamin;
subplot(3,1,1);
plot(x,T_nc,'linewidth',1);
xlabel('Nozzle domain');
ylabel('Temperatuer');
title({"Plot of variables at nozzle domain for non conservative form";"Temperature variation at nozzle domain"});
axis([0.5 3 0 max(T_nc)]);

subplot(3,1,2);
plot(x,rho_nc,'linewidth',1);
xlabel('Nozzle domain');
ylabel('Density');
title('Density variation at nozzle domain');
axis([0.5 3 0 max(rho_nc)]);

subplot(3,1,3);
plot(x,v_nc,'linewidth',1);
xlabel('Nozzle domain');
ylabel('Velocity');
title('Velocity variation at nozzle domain');
axis([0.5 3 0 max(v_nc)]);

% Plotting the mass flow rate at different time steps;
figure(3);
plot(x,mass_flow_nc(50,:),'linewidth',2,'r');
hold on;
plot(x,mass_flow_nc(100,:),'linewidth',2,'g');
hold on;
plot(x,mass_flow_nc(150,:),'linewidth',2,'b');
hold on;
plot(x,mass_flow_nc(600,:),'linewidth',2,'y');
hold on;
plot(x,mass_flow_nc(700,:),'linewidth',2,'c');
xlabel('nozzle domain');
ylabel('Mass flow rate');
title('Plot of mass flow rate through nozzle domain - Non conservative form');
legend("50dt","100dt","150dt","600dt","700dt");
axis([0 3 0 1.8]);

% Ploting for conservative form;

% Ploting of variables at throat;
figure(4);
subplot(4,1,1);
plot(rho_tr_c,'linewidth',2);
hold on
plot(rho_exact,'linewidth',2,'g');
xlabel('Time');
ylabel('rho/rho_0');
title({"Plot of variables at throat for conservative form";"Plot of temperature at throat"});
legend("Numerical values","Exact values");

subplot(4,1,2);
plot(pr_tr_c,'linewidth',2);
hold on;
plot(p_exact,'linewidth',2,'g');
xlabel('Time');
ylabel('p/p_0');
title('Pressure values at throat');

subplot(4,1,3);
plot(T_tr_c,'linewidth',2);
hold on;
plot(T_exact,'linewidth',2,'g');
xlabel('Time');
ylabel('T/T_0');
title('Temperature values at throat');

subplot(4,1,4);
plot(mach_tr_c,'linewidth',2);
hold on;
plot(Mach_exact,'linewidth',2,'g');
xlabel('Time');
ylabel('M/M_0');
title('Mach number at throat');

% Plotting of variables at nozzle domain;
figure(5);
subplot(3,1,1);
plot(x,T_c,'linewidth',1);
xlabel('Nozzle domain');
ylabel('Temperature');
title({"Plot of variables at nozzle domain for conservative form";"Temperature variation at nozzle domain"});
axis([0.5 3 0 max(T_c)]);

subplot(3,1,2);
plot(x,rho_c,'linewidth',1);
xlabel('Nozzle domain');
ylabel('Density');
title('Density variation at nozzle domain');
axis([0.5 3 0 max(rho_c)]);

subplot(3,1,3);
plot(x,v_c,'linewidth',1);
xlabel('Nozzle domain');
ylabel('Velocity');
title('Velocity variation at nozzle domain');
axis([0.5 3 0 max(v_c)]);

% Plotting of mass flow rates;
figure(6);
plot(x,mass_flow_c(50,:),'linewidth',2,'r');
hold on;
plot(x,mass_flow_c(150,:),'linewidth',2,'b');
hold on;
plot(x,mass_flow_c(200,:),'linewidth',2,'y');
hold on;
plot(x,mass_flow_c(500,:),'linewidth',2,'g');
hold on;
plot(x,mass_flow_c(700,:),'linewidth',2,'c');
xlabel('nozzle domain');
ylabel('Mass flow rate');
title('Plot of mass flow rate through nozzle domain - conservative form');
legend("50dt","150dt","200dt","500dt","700dt");
axis([0 3 0 1.8]);

% Printing out the computational time for conservative and non conservative form;
fprintf ('The computational time for conservative form of equation is %f sec', time_c);
fprintf ('The computational time for non conservative form of equation is %f sec', time_nc);

% Plot for comparing the mass flow rate of conservative form and non conservative form;
figure(7);
plot(x,mass_flow_c(1400,:),'linewidth',2,'r');
hold on;
plot(x,mass_flow_nc(1400,:),'linewidth',2,'b');
title('Comparing the mass flow rate of non conservative and conservative form');
xlabel('Nozzle domain');
ylabel('Mass flow');
legend("Conservative form","Non conservative form");

# Non conservative form: (Function for Non conservative form)

function [rho, T, v, rho_tr, T_tr, pr_tr, mach_tr, mass_flow, time_nc] = nonconserv_n(x,n,nt,dx)

% Defining the initial condition;
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);

thrt = find(a==1);

% Defining the value for gamma;
gamma = 1.4;

% Defining the dt;
% Assuming courant number, c = 0.5;
% dt = min((c*dx)/(sqrt(T)+v));

dt = min((0.5*dx)./(sqrt(T)+v));

% Creating a copy of all the above initial conditions;
rho_old= rho;
T_old = T;
v_old = v;
a_old = a;

% Solving 1d nozzle flow equation using macormack method;

% Non conservative form of fluid flow equations are given as;

% drho/dt = -rho*(dv/dx)-rho*v*d(loga)/dx-v*drho/dx, continuity equation;
% dv/dt =-v*(dv/dt)-1/gamma*(dT/dx+(T/rho*drho/dx)), Momentum equation;
% dT/dt = -v*dT/dx-(gamma-1)*T*(dv/dx+v*d(loga)/dx); Energy equation;

time_nc = cputime;

for k = 1:nt
rho_old = rho;
T_old = T;
v_old = v;

for i = 2:n-1

% Predictor step;

dvdx = (v(i+1)-v(i))/dx;
dadx = (log(a(i+1))-log(a(i)))/dx;
drhodx = (rho(i+1)-rho(i))/dx;
dTdx = (T(i+1)-T(i))/dx;

% Continuity equation;
drhodt_p(i) = (-1*rho(i)*dvdx)-(rho(i)*v(i)*dadx)-(v(i)*drhodx);

% Momentum equation;
dvdt_p(i) = (-1*v(i)*dvdx)-((1/gamma)*((dTdx)+((T(i)/rho(i))*drhodx)));

% Energy equation;
dTdt_p(i) = (-1*v(i)*dTdx)-((gamma-1)*T(i)*(dvdx+(v(i)*dadx)));

% Solution update;
rho(i) = rho_old(i)+(drhodt_p(i)*dt);
v(i) = v_old(i)+(dvdt_p(i)*dt);
T(i) = T_old(i)+(dTdt_p(i)*dt);

end

% Corrector step;

for i = 2:n-1;

dvdx = (v(i)-v(i-1))/dx;
dadx = (log(a(i))-log(a(i-1)))/dx;
drhodx = (rho(i)-rho(i-1))/dx;
dTdx = (T(i)-T(i-1))/dx;

% Continuity equation;
drhodt_c(i) = (-1*rho(i)*dvdx)-(rho(i)*v(i)*dadx)-(v(i)*drhodx);

% Momentum equation;
dvdt_c(i) = (-1*v(i)*dvdx)-((1/gamma)*((dTdx)+((T(i)/rho(i))*drhodx)));

% Energy equation;
dTdt_c(i) = (-1*v(i)*dTdx)-((gamma-1)*T(i)*(dvdx+(v(i)*dadx)));

end

% Finding out the average derivative;
for m = 2:n-1
drhodt(m) = 0.5*(drhodt_p(m)+drhodt_c(m));
dvdt(m) = 0.5*(dvdt_p(m)+dvdt_c(m));
dTdt(m) = 0.5*(dTdt_p(m)+dTdt_c(m));
end
% Final update;

for j = 2:n-1
rho(j) = rho_old(j)+(drhodt(j)*dt);
v(j) = v_old(j)+(dvdt(j)*dt);
T(j) = T_old(j)+(dTdt(j)*dt);
end

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

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

% calculating the non dimesional variables;
mass_flow(k,:) = rho.*a.*v;

pr = rho.*T;

mach_num = v./sqrt(T);

T_tr(k) = T(thrt);

rho_tr(k) = rho(thrt);

pr_tr(k) = pr(thrt);

mach_tr(k) = mach_num(thrt);

end

end

# Conservative form: (Function for conservative form)

function [rho, T, v, rho_tr, T_tr, pr_tr, mach_tr, mass_flow, time_c] = conserv_n(x,n,L,nt,dx)

% Defining initial conditions;
% Initial condition for density and temperature;
for i = 1:n
if x(i)>= 0 && x(i)<=L/6
rho(i) = 1;
T(i) = 1;
end
if x(i)>=L/6 && x(i)<=L/2
rho(i) = 1 -0.366*(x(i)-0.5);
T(i) = 1 -0.167*(x(i)-0.5);
end
if x(i)>=L/2
rho(i) = 0.634-0.3879*(x(i)-1.5);
T(i) = 0.833-0.3507*(x(i)-1.5);
end
end

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

thrt = find(a==1);

% Defining the initial condition for velocity;
v = 0.59./(rho.*a);

% Defining the value of gamma;
gamma = 1.4;

% Defining the solution vectors;
U1 = rho.*a;

U2 = rho.*a.*v;

U3 = U1.*((T./(gamma-1))+((gamma/2).*v.^2));

time_c = cputime;

% Solving using macormack's method;
for k = 1:nt

U1_old = U1;
U2_old = U2;
U3_old = U3;

dt = min((0.5*dx)./(sqrt(T)+v));
% Defining flux vector's;
F1 = U2;
F2 = ((U2.^2)./U1) + (((gamma-1)/gamma)*(U3 - ((gamma*0.5*(U2.^2))./U1)));
F3 = ((gamma*U2.*U3)./U1) - ((gamma*0.5*(gamma-1)*(U2.^3))./(U1.^2));

% predictor method
for i = 2:n-1

J2 = (1/gamma)*rho(i)*T(i)*((a(i+1)-a(i))/dx);

% dU1/dt = -dF1/dx - continuity equation;
% dU2/dt = -dF2/dx+j2 - Momentum equation;
% dU3/dt = -dF3/dx - Energy equation;

dU1dt_p(i) = -1*((F1(i+1)-F1(i))/dx);

dU2dt_p(i) = (-1*(F2(i+1)-F2(i))/dx) + J2;

dU3dt_p(i) = -1*((F3(i+1)-F3(i))/dx);

% Solution update;
U1(i) = U1_old(i)+(dU1dt_p(i)*dt);
U2(i) = U2_old(i)+(dU2dt_p(i)*dt);
U3(i) = U3_old(i)+(dU3dt_p(i)*dt);
end

% Updating the flux terms for corrector step;
F1 = U2;
F2 = ((U2.^2)./U1) + (((gamma-1)/gamma)*(U3 - ((gamma*0.5*(U2.^2))./U1)));
F3 = ((gamma*U2.*U3)./U1) - ((gamma*0.5*(gamma-1)*(U2.^3))./(U1.^2));

% Updating primitive variables for J2;
rho = U1./a;

v = U2./U1;

T = (gamma-1)*((U3./U1)-((gamma/2)*(v.^2)));

% Corrector step;
for j = 2:n-1

J2 = (1/gamma)*rho(j)*T(j)*((a(j)-a(j-1))/dx);

dU1dt_c(j) = -1*((F1(j)-F1(j-1))/dx);

dU2dt_c(j) = (-1*((F2(j)-F2(j-1))/dx)) + J2;

dU3dt_c(j) = -1*((F3(j)-F3(j-1))/dx);

end

% Calculating the average;

for m = 2:n-1
dU1dt_avg(m) = 0.5*(dU1dt_p(m)+dU1dt_c(m));
dU2dt_avg(m) = 0.5*(dU2dt_p(m)+dU2dt_c(m));
dU3dt_avg(m) = 0.5*(dU3dt_p(m)+dU3dt_c(m));
end

% Final update;
for i = 2:n-1
U1(i) = U1_old(i)+(dU1dt_avg(i)*dt);
U2(i) = U2_old(i)+(dU2dt_avg(i)*dt);
U3(i) = U3_old(i)+(dU3dt_avg(i)*dt);
end

% Applying the boundary conditions;
% inlet;

U1(1) = rho(1)*a(1);
U2(1) = 2*U2(2) - U2(3);
U3(1) = U1(1)*((T(1)/(gamma-1))+(0.5*gamma*(v(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);

% Extracting primitive variables from solution vector;
rho = U1./a;
v = U2./U1;
T = (gamma-1)*((U3./U1) - ((gamma/2)*(v.^2)));

mass_flow(k,:) = rho.*a.*v;

pr = rho.*T;

mach_num = v./sqrt(T);

rho_tr(k) = rho(thrt);

pr_tr(k) = pr(thrt);

T_tr(k) = T(thrt);

mach_tr(k) = mach_num(thrt);

end

end

% Plots:

# Plots of Non conservaive form:

1. Steady state distribution of primitive variables inside the nozzle:

2. Time wise variation of primitve variables ( In here the time wise variation of primitive variables at throat is plotted)

3. Variation of mass flow rates inside the nozzle at different ime step:

# Plots of conservative form:

1. Steady state distribution of primitive variables inside the nozzle:

2. Time wise variation of primitve variables ( In here the time wise variation of primitive variables at throat is plotted)

3 . Variation of mass flow rates inside the nozzle at different ime step:

Explaination of the plots:

* From the plots of Variation of variables in the nozzle domain (Both conservative and non cnservative form) it is clear that temperature is droping and the velocity is increasing, this is what exactly a supersonic nozzle does it converts the heat energy into kinetic energy.

* Comparing the plots of variables taken at the throat, in non conservative form the variables have settled to steady state (Exact values) after 500 time step. In conservative form there is a difference in between the exact value and numerical values, this is because instead of primitive variables solution vectors and flux vectors are used in the macCormack method. But still even the numerical values are getting converged after 400 time step.

* Comparing the mass flow rates of conservative form and non conservative form, it is clear from the graph that the mass flow is becoming steady at about 500 time step.

# Grid dependence test: (Values taken at throat at 1400 time step)

	rho	Pressure	Temp	Mach No.
Exact values	0.634	0.528	0.833	1
Non conservative form (n = 31)	0.6385	0.53454	0.83718	0.99880
Conservative form (n = 31)	0.65039	0.54635	0.84005	0.98173
Conservative form (n = 61)	0.63821	0.53302	0.83517	0.99935

It is clear from the above observation that non conservative form at n = 31 gives promising results i.e, results are close to exact values. But at the same time conservative form at n = 31 is not giving results close to the exact values there is a considerable amount of difference between the exact values and analytical values (refer the above Time wise variation of primitve variables plot for conservative form).

For conservative form at n = 61 the values are very much close to the exact values. A plot of time wise variation of primitive variables at throat for conservative form at n = 61 is provided below from which it is very much clear that the numerical values are converging into exact values at n = 61 and also this is much better compared to n = 31.

Below is the plot of primitive variables at throat for conservative form at n = 61, from which it is clear that at n = 61 the primitve variables are very much equal to exact values.

4. Comparison of normalised mass flow rates of non conservative form and non conservative form:

above plot is taken at 1400 time step for both conservative and non conservative form ( n = 31)

and it is clear that for conservative form the mass flow rate is very much stable whereas for non conservative form it is not stable.

# Comparing the computational time:

above are the screenshots of the output from the command window, it is well understood from the above observation that time taken by conservative form to do the computation is more beacuse it takes time to retrive or extract the primitive variables from the solution vectors. Moreover above observation of compuational time is done for n = 31. For n = 61 where we are getting accurate details it will ake much more time.

# Conclusion:

1. With less number of grid points if we want very much accurate result then, non conservative form is the best shot.

2. comparing the mass flow rates, if we want much stable mass flow rate then conservative form will be a good to go option.

3. Computational time wise non conservative form is way ahead than conservative form.

so overall we can say that non conservative form is much better compared to conservative form.