Quasi 1D transient super-sonic nozzle flow simulation using Macormack Method

OBJECTIVES -

1. Developing MATLAB based solvers to simulate quasi 1D supersonic nozzle flows using both the non conservative and conservative forms of the governing equations.
2. Determine the steady state distributions of the flow field variables (density, temperature, velocity etc.) for both the forms.
3. Determine the timewise variation of the flow field variables at a particular position in the nozzle for both the forms.
4. Determine the mass flow rate variation with time at a particular location in the nozzle.
5. Compare the results obtained by solving both the forms.
6. Perform grid independence test.

TOOL AND TECHNIQUE USED -

MATLAB R2016a was used to build the solvers for solving both the forms of the equations.

MacCormack's method of Finite difference discretization was used due to the hyperbolic nature of the governing equations.

ASSUMPTIONS -

1. The flow properties are uniform at a given cross section. Although the area varies with nozzle length x (which makes it a 2D flow in reality), the flow properties are asumed to vary only with x. Thus the flow is defined as Quasi 1D flow.
2. `gamma` is the specific heat ratio. Air has been considered as the working fluid(`gamma = 1.4`)

GOVERNING EQUATIONS AND METHODOLOGY -

A. The physical problem - Subsonic Supersonic Isentropic Flow

A steady, isentropic flow through a convergent - divergent nozzle (CD Nozzle) has been considered. A sketch of the nozzle has been shown below (Ref. Computational Fluid Dynamics by JD Anderson).

At inlet, the flow comes from a reservoir (pressure `p_0`, temperature `T_0`). `p_0 and T_0` are total pressure and total temperature, respectively.

The flow expands isentropically to supersonic speeds at the nozzle exit. The pressure, temperature, velocity and temperature at the exit are `p_e,T_e,V_e and M_e` respectively.

Thus, the flow is locally subsonic in the convergent section, sonic at the throat and supersoninc in the divergent section.

Reason for using the non dimensional forms of the governing equations - We can also use the dimensional form of the governing equations, if we want to determine the magnitudes of the flow variables. However, our aim is to do a qualitative study of the variation of flow field variables. Thus, we use non dimensional variables `rho/rho_0, T/T_0 and p/p_0` which vary between 0 and 1. Thus, we use the following non dimensional variables(with primes as superscript).

`T' = T/T_0;rho' = rho/rho_0; x' = x/L; V' = V/a_0 and t' = t/(L/a_0) and A' = A/A*`

Here `T, rho, x,V,t`are the dimensional temperature, density, nozzle length, velocity and time respectively.

`T_0, rho_0, a_0 (=sqrt(gammaRT_0))` are the reservoir temperature, density and speed of sound.

A* is the sonic throat area.

Non Convervative form of the quasi 1D flow -

Continuity equation

`(delrho')/(delt') = -rho'(delV')/(delx') - rho'V'(del(ln A'))/(delx') - V'(delrho')/(delx')`...(1)

Momentum equation

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

Energy equation

`(delT')/(delt') = -V'(delT')/(delx') - (gamma-1)T'((delV')/(delx')+V'(del(ln A'))/(delx'))` ......(3)

Convervative form of the quasi 1D flow -

Contiuity equation

`(del(rho'A'))/(delt') + (del(rho'A'V'))/(delx') = 0` .....(4)

Momentum equation

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

Energy equation

`(del(rho'(e/(gamma-1)+gamma/2V'^2)A'))/(delt') + (del(rho'(e/(gamma-1)+gamma/2V'^2)A'V' + p'A'V'))/(delx') = 0` .....(6`

In a simple generic form, the equations 4,5 and 6 are expressed as

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

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

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

respectively.

`U_1,U_2 and U_3` are the elements of the solutions vector U.

`F_1, F_2 and F_3` are the elements of the flux vector F.

`J_2`is the element of the source vector J.

`U_1 = rho'A'`

`U_2 = rho'A'V'`

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

`F_1 = rho'A'V'`

`F_2 = rho'A'V'^2 + 1/gammap'A'`

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

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

Numerical discretization technique using Maccormack method -

Due to the hyperbolic nature of the problem, the numerical discretization has been done using the Finite difference base Maccormack's explicit technique.

As with any numerical discretization method, we divide the domain into finite number of grod points (n).

The MacCormack's technique is a predictor -corrector method. The approach has been demonstrated below using the non conservation equation. The same approach ahs been utilized in the conservation form solution but instead of solving for dlow variables directly, we solve for the solutions vector and then detemine the flow variables using the results.

All the variables are non dimensional. The use of prime has been dropped to reduce the complexity of notation.

1. Predictor StepThe spatial derivatives are discretized using the Forward differencing method.

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

2. Intermediate solution update

`rho_(i,p)^(t+Deltat) = rho_(i)^t +((delrho)/(delt))_(i,p)^tDeltat`

3. Corrector Step

`((delrho)/(delt))_(i,c)^t = -rho_(i,p)^(t+Deltat)(V_(i+1)^t - V_(i)^t)/(Deltax) -rho_(i,p)^(t+Deltat)V_i^t(lnA_(i+1) - lnA_i)/(Deltax)-V_i^t(rho_(i+1,p)^(t+Deltat) - rho_(i,p)^(t+Deltat))/(Deltax)`

4. Average time derivatives

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

5. Final solution update

`rho_(i)^(t+Deltat) = rho_(i)^t +((delrho)/(delt))_(i,av)^tDeltat`

Here the first term on the right side is the value of density ratio at the start of that particular time step.

Note - The above 5 steps are repeated for all the flow field variables (solutions vecotor in case of conservation form).

Calculation of time step -

For the hyperbolic equations, as per the stability criterion a constraint exists on the time step size. It is given as -

`Deltat = CFL((Deltax)/(a+V))`

Here CFL is the Courant - Friedrichs–Lewy number ( = 0.5, in both forms).

`Deltax` is the grid spacing

`a`is the local speed of sound

V is the local velocity.

In non dimensional form, it becomes,

`Deltat = CFL((Deltax)/(sqrtT+V))`

At the start of each time step, the valued of `Deltat`is calculated at all the grid points, and the minimum value is chosen for use in all the equations i.e.

`Deltat = min(Deltat_1^t,Deltat_2^t,Deltat_3^t,Deltat_4^t .... Deltat_n^t)`

This gives a time accurate solution of the actual transient flow field, consistent with the unsteady, continuity, momentum and energy equations.

Boundary Conditions -

All the flow variables are non dimensional.

1. Non Conservative Form

Flow Variable	Value
`V_1`	`2V_2 - V_3`
`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

Solution vector element (or flow variable)	Value
`rho_1`	1.0
`T_1`	1.0
`U_1(1)`	`rho.A = 1`
`U_2(1)`	`2U_2(2) - U_2(3)`
`v_1`	`U_2/U_1`
`U_3(1)`	updated using U1 and U2
`U_1(n)`	`2U_1(n-1)-U_1(n-2)`
`U_2(n)`	`2U_2(n-1))-U_2(n-2)`
`U_3(n)`	`2U_3(n-1)-U_3(n-2)`

Initial Conditions-

All the flow variables are non dimensional.

1. Non Conservative form

Nozzle shape -

`A = 1+2.2(x-1.5)^2` where `x in [0,3]`

`rho = 1-0.3146x`

`T = 1-0.2314x`

`V =(0.1+1.09x)sqrtT`

2. Conservative Form

For x<=0.5

`rho = 1`

`T = 1`

For 0.5<x<=1.5

`rho = 1-0.366(x-0.5)`

`T = 1-0.167(x-0.5)`

For 1.5<x<=3.5

`rho = 0.634-0.3879(x-1.5)`

`rho = 0.833-0.3507(x-1.5)`

Thus, for initial conditions only,

`U_2 = 0.59`

Thus,

`V = 0.59/(rho*A)`

RESULTS AND DISCUSSIONS -

In all the plots shown below, the solid lines represent solution by non servative form and the dashed lines (black) represent the solution by conservative form.

1. Steady state distribution of the flow variables

Steady state solution can be either decided by a convergence criterion or by stopping the solution after a large numer of time steps. In this project, we have used the latter method and the solution achieves steady state after about 500 time steps. However, for a check, the number of time steps chosen are 1000.

From the above graph, it is clear that for this problem there is negligible difference betweeen the solution of the flow variables by any of the two forms.

2. Timewise variation of the flow variables along the nozzle

This plot compares the timewise variation of density ratio, temperature ratio and velocity ratio at two locations in the nozzle. The left column of the graphs is at the throat while the right column is at the nozzle exit.

After observing the graphs at both the locations, it can be conculuded that the solution by conservation form of equations show little deviation from the steady state value. These deviations are quite high and oscillatory in nature when solving the non conservation form of equations. The quantitaive values have been mentioned when discussing the grid independence test results.

3. Normalized mass flow rate variation at different time steps

This depicts the most important difference in the two approaches.

The initial condition of constant mass flow gives quite interesting results towards the steady state as evident by the conservation form plot shared above. The normalized mass flow rate in the conservation form is very close to the exact value of 0.579 at the steady state.

When compared with the non conservation form, we observe that the conservaiton form mass flow rate has very less variations. However, this can be attributed to the fact that a constant initial mass flow has been assumed in the conservation form.

4. Comparison of mass flow rate at steady state

We observe much satisfactory values of normalized mass flow rate by the conservation form of equation because of the following reasons -

1. It gives a distribution which is much closer to being constant. On contrast, the non conservation form has spurious variations at the inflow and outflow boundaries.

2. The steady state mass flow results are much closer to the exact values (=0.579).

This is because the mass flow is itself a dependent variable(U2) in the conservation form while it is a secondary result in the non conservation form.

However, this can't be used as a criterion to establish that conservation form is numerically superior to the non conservation form. When comparing the flow variables (values in the grid independence test results), we observe that the flow field variables from non conservation form are much closer to the exact values than those obtained from the conservation form of equations.

5. Grid Independence Test

All the values are at the throat of the nozzle.

Form of Equation	Grid points	`rho/rho_0`	`T/T_0`	`p/p_0`	M
Non Conservation	31	0.6421	0.8380	0.5381	0.9995
Conservative	31	0.6675	0.8441	0.5634	0.9733
Non Conservative	61	0.6412	0.8372	0.5368	0.9997
Conservative	61	0.6480	0.8387	0.5434	0.9877
Exact	-	0.634	0.833	0.528	1.0000

It can clearly be concluded from the above table that the values of flow variables are closer to the exact values in the non conservative form solution.

Upon comparing the results of non conservative and conservative forms at different number of grid points, it can be concluded that grid convergence is achieved at 31 grid points.

MATLAB Codes

1. Main Algorithm

%Quasi 1D flow in a CD Nozzle
%Macormack method has been used
%NON DIMENSIONAL Non conservative form of the GDE used
%Chapter 7, CFD by JD ANderson

%Flow is isentropic
%Time marchnig
% All the primal variables are non dimensional

clear all
clc

%Input parameters
n = 31;     %Number of nodes as used in JD Anderson
x = linspace(0,3,n); %One dimensional mesh
dx = x(2)-x(1); %Grid spacing
gamma = 1.4;    %specific heat ratio

%Initial conditions. As given in JD Anderson
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; %Nozzle cross section
p = rho.*T;     %Pressure ratio


%Initial conditions for conservative forms.Conservative form
rho_con = zeros(1,n);
T_con = zeros(1,n);
v_con = zeros(1,n);
p_con = zeros(1,n);     %Pressure ratio
for z = 1:n
    if(x(z)>=0 && x(z)<=0.5)
        rho_con(z) = 1;
        T_con(z) = 1;
    elseif(x(z)>0.5 && x(z)<=1.5)
        rho_con(z) = 1-0.366*(x(z)-0.5);
        T_con(z) = 1-0.167*(x(z)-0.5);
    elseif(x(z)>1.5 && x(z)<=3.5)
        rho_con(z) = 0.634-0.3879*(x(z)-1.5);
        T_con(z) = 0.833-0.3507*(x(z)-1.5);
    end
end
v_con = 0.59./(rho_con.*a);
p_con = rho_con.*T_con;
%Solution variables of Conservative form
U1 = rho_con.*a;
U2 = rho_con.*v_con.*a;
U3 = rho_con.*((T_con/(gamma-1))+(0.5*gamma*(v_con.^2))).*a;
%Time step and step size
nt = 1400;
CFL = 0.5;      %Used for calculating the time step size

%Location of minimum area (index)
k = find(a==1);

%Variables to store time distribution data at thorat and nozzle exit.
%First two columns for non conservative variables. Last two for
%conservative  variables.
rho_t = zeros(nt,4);
T_t = zeros(nt,4);
v_t = zeros(nt,4);
p_t = zeros(nt,4);
mf = ones(nt,1)*(rho.*v.*a);     %to store mass flow timewise distribution
mf_con = ones(nt,1)*(rho_con.*v_con.*a);     %to store mass flow timewise distribution

for time = 1:nt
    [rho,T,v,M,p] = non_cons(n,dx,rho,T,v,a,gamma,CFL);
    %Capturing time wise variation at throat (given by k) and nozzle exit.
    rho_t(time,1) = rho(k); rho_t(time,2) = rho(end);
    T_t(time,1) = T(k); T_t(time,2) = T(end);
    v_t(time,1) = v(k); v_t(time,2) = v(end);
    p_t(time,1) = p(k); p_t(time,2) = p(end);
    mf(time,:) = rho.*v.*a;
    [rho_con,T_con,v_con,M_con,p_con] = cons(n,dx,rho_con,T_con,v_con,a,gamma,CFL);
    %Capturing time wise variation at throat (given by k) and nozzle exit.
    rho_t(time,3) = rho_con(k); rho_t(time,4) = rho_con(end);
    T_t(time,3) = T_con(k); T_t(time,4) = T_con(end);
    v_t(time,3) = v_con(k); v_t(time,4) = v_con(end);
    p_t(time,3) = p_con(k); p_t(time,4) = p_con(end);
    mf_con(time,:) = rho_con.*v_con.*a;
    
end

%% Plotting the steady state results
figure(1)
subplot(3,1,1)
plot(x,rho,'LineWidth',1.5);
hold on
plot(x,rho_con,'--k','LineWidth',1.5);
legend('Non Conservative Form','Conservative Form');
xlabel('Non dimensional nozzle length');
ylabel('Density ratio');
title(sprintf('Steady state distribution - Density ratio (Time steps = %0d, CFL = %0.2f, Grid points = %0d) ',time,CFL,n));
grid on
subplot(3,1,2)
plot(x,T,'LineWidth',1.5);
hold on
plot(x,T_con,'--k','LineWidth',1.5);
legend('Non Conservative Form','Conservative Form');
xlabel('Non dimensional nozzle length');
ylabel('Temperature ratio');
title(sprintf('Steady state distribution - Temperature ratio (Time steps = %0d, CFL = %0.2f, Grid points = %0d) ',time,CFL,n));
grid on;
subplot(3,1,3)
plot(x,v,'LineWidth',1.5);
hold on
plot(x,v_con,'--k','LineWidth',1.5);
legend('Non Conservative Form','Conservative Form');
xlabel('Non dimensional nozzle length');
ylabel('Velocity ratio');
title(sprintf('Steady state distribution - Velocity ratioTime steps = %0d, CFL = %0.2f, Grid points = %0d) ',time,CFL,n));
grid on;

%% Plotting the time distribution data
figure(2)
subplot(3,2,1)
plot(1:nt,rho_t(:,1),'LineWidth',1.5);     %At throat
hold on
plot(1:nt,rho_t(:,3),'--k','LineWidth',1.5);
legend('Non Conservative Form','Conservative Form');
xlabel('Number of time steps');
ylabel('Density ratio');
title(sprintf('Timewise variation of Density ratio at nozzle throat(CFL = %0.2f) ',CFL));
grid on
subplot(3,2,2)
plot(1:nt,rho_t(:,2),'LineWidth',1.5);     %At exit
hold on
plot(1:nt,rho_t(:,4),'--k','LineWidth',1.5);
legend('Non Conservative Form','Conservative Form');
xlabel('Number of time steps');
ylabel('Density ratio');
title(sprintf('Timewise variation of Density ratio at nozzle exit(CFL = %0.2f) ',CFL));
grid on

subplot(3,2,3)
plot(1:nt,T_t(:,1),'LineWidth',1.5);     %At throat
hold on
plot(1:nt,T_t(:,3),'--k','LineWidth',1.5);
legend('Non Conservative Form','Conservative Form');
xlabel('Number of time steps');
ylabel('Temperature ratio');
title(sprintf('Timewise variation of Temperature ratio at nozzle throat(CFL = %0.2f) ',CFL));
grid on
subplot(3,2,4)
plot(1:nt,T_t(:,2),'LineWidth',1.5);     %At exit
hold on
plot(1:nt,T_t(:,4),'--k','LineWidth',1.5);
legend('Non Conservative Form','Conservative Form');
xlabel('Number of time steps');
ylabel('Temperature ratio');
title(sprintf('Timewise variation of Temperature ratio at nozzle exit(CFL = %0.2f) ',CFL));
grid on

subplot(3,2,5)
plot(1:nt,v_t(:,1),'LineWidth',1.5);     %At throat
hold on
plot(1:nt,v_t(:,3),'--k','LineWidth',1.5);
legend('Non Conservative Form','Conservative Form');
xlabel('Number of time steps');
ylabel('Velocity ratio');
title(sprintf('Timewise variation of Velocity ratio at nozzle throat(CFL = %0.2f) ',CFL));
grid on
subplot(3,2,6)
plot(1:nt,v_t(:,2),'LineWidth',1.5);     %At exit
hold on
plot(1:nt,v_t(:,4),'--k','LineWidth',1.5);
legend('Non Conservative Form','Conservative Form');
xlabel('Number of time steps');
ylabel('Velocity ratio');
title(sprintf('Timewise variation of Velocity ratio at nozzle exit(CFL = %0.2f) ',CFL));
grid on

%% Plotting timewise variation of mass flow rate along the nozzle
figure(3)
subplot(1,2,1)
plot(x,mf(1,:),'--');     %At start of the simulation i.e. 0*dt
leg1 = strcat('Mass flow rate at time step = ',num2str(1));
hold on
plot(x,mf(50,:));    %At time step = 50
leg2 = strcat('Mass flow rate at time step = ',num2str(50));
hold on
plot(x,mf(100,:));   %At time step = 100
leg3 = strcat('Mass flow rate at time step = ',num2str(100));
hold on
plot(x,mf(150,:));
leg4 = strcat('Mass flow rate at time step = ',num2str(150));
hold on
plot(x,mf(200,:));
leg5 = strcat('Mass flow rate at time step = ',num2str(200));
hold on
plot(x,mf(end,:));
leg6 = strcat('Mass flow rate at time step = ',num2str(nt));
legend(leg1,leg2,leg3,leg4,leg5,leg6);
title({'Variation of Mass flow rate distribution along the nozzle at different time steps',...
    'Solution of the Non Conservation Forms of Governing Equations'});
grid on

subplot(1,2,2)
plot(x,mf_con(1,:),'--');     %At start of the simulation i.e. 0*dt
leg1 = strcat('Mass flow rate at time step = ',num2str(1));
hold on
plot(x,mf_con(50,:));    %At time step = 50
leg2 = strcat('Mass flow rate at time step = ',num2str(50));
hold on
plot(x,mf_con(100,:));   %At time step = 100
leg3 = strcat('Mass flow rate at time step = ',num2str(100));
hold on
plot(x,mf_con(150,:));
leg4 = strcat('Mass flow rate at time step = ',num2str(150));
hold on
plot(x,mf_con(200,:));
leg5 = strcat('Mass flow rate at time step = ',num2str(200));
hold on
plot(x,mf_con(end,:));
leg6 = strcat('Mass flow rate at time step = ',num2str(nt));
legend(leg1,leg2,leg3,leg4,leg5,leg6);
title({'Variation of Mass flow rate distribution along the nozzle at different time steps',...
    'Solution of the Conservative Forms of Governing Equations'});
grid on

%% Comparison of normalized mass flow of both the solutions at steady state
figure(4)
plot(x,mf(end,:));
hold on
plot(x,mf_con(end,:));
title({'Comparison of normalized mass flow rate distributions obtained with',...
    'non conservation and conservation forms of governing equations at Steady State'});
xlabel('Non dimensional nozzle length');
ylabel('Non dimensional mass flow rate');
legend('Non Conservation form Mass flow distribution','Conservative form Mass flow distribution');
grid on

2. Non Conservative form solver

function [rho,T,v,M,p]= non_cons(n,dx,rho,T,v,a,gamma,CFL)

rho_old = rho;
T_old = T;
v_old = v;
%Calculation of time step as per CFL number. Derivation in JD anderson
dt = CFL*dx./((T.^0.5)+v);
dt = min(dt);
for i = 2:n-1
    %Predictor Step
    drho_dx = (rho(i+1)-rho(i))/dx;
    dv_dx = (v(i+1)-v(i))/dx;
    dlna_dx = (log(a(i+1))-log(a(i)))/dx;
    dT_dx = (T(i+1)-T(i))/dx;
    
    %Continuity eq
    drho_dt_p = -(rho(i)*dv_dx)-(rho(i)*v(i)*dlna_dx)-(v(i)*drho_dx);
    %Momentum eq
    dv_dt_p = -(v(i)*dv_dx)-((dT_dx + (T(i)*drho_dx/rho(i)))/gamma);
    %Energy eq
    dT_dt_p = -(v(i)*dT_dx)-((gamma-1)*T(i)*(dv_dx+(v(i)*dlna_dx)));
    
    %Solution update after predictor step
    rho(i) = rho(i) + drho_dt_p*dt;
    T(i) = T(i) + dT_dt_p*dt;
    v(i) = v(i) + dv_dt_p*dt;
    
    %Corrector step
    drho_dx = (rho(i)-rho(i-1))/dx;
    dv_dx = (v(i)-v(i-1))/dx;
    dlna_dx = (log(a(i))-log(a(i-1)))/dx;
    dT_dx = (T(i)-T(i-1))/dx;
    
    %Continuity eq
    drho_dt_c = -(rho(i)*dv_dx)-(rho(i)*v(i)*dlna_dx)-(v(i)*drho_dx);
    %Momentum eq
    dv_dt_c = -(v(i)*dv_dx)-((dT_dx + (T(i)*drho_dx/rho(i)))/gamma);
    %Energy eq
    dT_dt_c = -(v(i)*dT_dx)-((gamma-1)*T(i)*(dv_dx+(v(i)*dlna_dx)));
    
    %Average values
    drho_dt_avg = 0.5*(drho_dt_p+drho_dt_c);
    dv_dt_avg = 0.5*(dv_dt_p+dv_dt_c);
    dT_dt_avg = 0.5*(dT_dt_p+dT_dt_c);
    
    %Final solution update
    rho(i) = rho_old(i) + drho_dt_avg*dt;
    T(i) = T_old(i) + dT_dt_avg*dt;
    v(i) = v_old(i) + dv_dt_avg*dt;
end
%BCs. From JD Anderson
v(1) = 2*v(2)-v(3);   %Inlet
v(n) = 2*v(n-1)-v(n-2); %Outlet
rho(n) = 2*rho(end-1)-rho(n-2); %Outlet
T(n) = 2*T(n-1)-T(n-2); %Outlet

M = v./sqrt(T);
p = rho.*T;

3. Conservative Form solver

%To solve conservative form of equations for quasi 1D supersoninc nozzle.
function [rho_con,T_con,v_con,M_con,p_con] = cons(n,dx,rho_con,T_con,v_con,a,gamma,CFL)
%Refer JD Anderson chapter 7 for conservative form
%Continuity eq - delU1/delt = -delF1/delx
%Momnetum eq - delU2/delt = -delF2/delx + J2
%Energy eq - delU3/delt = -delF3/delx

U1 = rho_con.*a;
U2 = rho_con.*v_con.*a;
U3 = rho_con.*((T_con/(gamma-1))+(0.5*gamma*(v_con.^2))).*a;

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

dt = CFL*dx./((T_con.^0.5)+v_con);
dt = min(dt);

k = (gamma-1)/gamma;

%Inlet BCs. Outlet BCs at the end of the iterations on all nodes.
U2(1) = 2*U2(2) - U2(3); %Inlet
v_con(1) = U2(1)/U1(1);
U3(1) = U1(1)*((T_con(1)/(gamma-1))+(0.5*gamma*v_con(1)^2));

for j = 2:n-1
    %Predictor step
    F1 = U2;
    F2 = ((U2.^2)./U1) + (k*(U3-(0.5*gamma*(U2.^2)./U1)));
    F3 = gamma*((U2.*U3)./U1) - (0.5*gamma*(gamma-1)*((U2.^3)./(U1.^2)));
    
    del_lna_dx = (log(a(j+1))-log(a(j)))/dx;
    %     J2 = (rho_con.*T_con.*dela_dx)/gamma;
    J2 = k*(U3-(0.5*gamma*(U2.^2)./U1))*del_lna_dx;
    dU1_dt_p = -(F1(j+1)-F1(j))/dx;
    dU2_dt_p = -(F2(j+1)-F2(j))/dx + J2(j);
    dU3_dt_p = -(F3(j+1)-F3(j))/dx;
    
    %Solution update
    U1(j) = U1(j) + dU1_dt_p*dt;
    U2(j) = U2(j) + dU2_dt_p*dt;
    U3(j) = U3(j) + dU3_dt_p*dt;
    
    %Corrector step
    F1 = U2;
    F2 = (U2.^2./U1) + (((gamma-1)*(U3-(0.5*gamma*(U2.^2)./U1)))/gamma);
    F3 = gamma*((U2.*U3)./U1) - (0.5*gamma*(gamma-1)*((U2.^3)./(U1.^2)));
    
    del_lna_dx = (log(a(j))-log(a(j-1)))/dx;
    J2 = k*(U3-(0.5*gamma*(U2.^2)./U1))*del_lna_dx;
    dU1_dt_c = -(F1(j)-F1(j-1))/dx;
    dU2_dt_c = -(F2(j)-F2(j-1))/dx + J2(j);
    dU3_dt_c = -(F3(j)-F3(j-1))/dx;
    
    %Averaging the derivatives
    dU1_dt_avg = 0.5*(dU1_dt_p+dU1_dt_c);
    dU2_dt_avg = 0.5*(dU2_dt_p+dU2_dt_c);
    dU3_dt_avg = 0.5*(dU3_dt_p+dU3_dt_c);
    
    %Final Solution update
    U1(j) = U1_old(j) + dU1_dt_avg*dt;
    U2(j) = U2_old(j) + dU2_dt_avg*dt;
    U3(j) = U3_old(j) + dU3_dt_avg*dt;
end

%Outlet BCs
U1(n) = 2*U1(n-1)-U1(n-2); %Outlet
U2(n) = 2*U2(n-1)-U2(n-2); %Outlet
U3(n) = 2*U3(n-1)-U3(n-2); %Outlet

%Primal Variables
rho_con = U1./a;
v_con = U2./U1;
T_con = (gamma-1)*((U3./U1)-(0.5*gamma*v_con.^2));

M_con = v_con./sqrt(T_con);
p_con = rho_con.*T_con;