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

Objective:

• Steady-state distribution of primitive variables inside the nozzle.
• Time-wise variation of primitive variables.
• Variation of mass flow rate distribution inside the nozzle at the different time steps during the time-marching process.
• Comparision of normalized mass flow rate distribution of both forms.

Theory:

A nozzle is a device designed to control the direction or characteristics of a fluid flow (especially to increase velocity) as exits (or enters) an enclosed chamber or pipe. A nozzle is often a pipe or tube of varying cross-sectional area, and it can be used to direct or modify the flow of a fluid (liquid or gas). nozzles are frequently used to control the rate of flow, speed, direction, mass, shape,and / or pressue of the stream that emerges from them.

There are different types of nozzle but here in our project we uses a convergent Divergent nozzle. This nozzle is used to obtain a supersonic stream starting from low speed at the inlet. Then the mach number should increase from Ma=0 near the inlet to Ma.1 at the exit. it is clear that the nozzle must converge in subsonic portion and diverge in the supersonic portion. Such a nozzle is called a convergent-divergent nozzle.

Here we considering the steady, isentropic flow through a convergent-divergent nozzle as given below, The flow at the inlet to the nozzle comes from a reservoir where the pressure and temperature are denoted by `P_0 and T_0`. The cross-sectional area of the reservior is large `(A=infty)`, and hence the velocity is very small `(V=0)`. Thus `p_0 and T_0` are stagnation values, or total pressure and total temperature, respectively. The flow expands isentropically to supersonic speeds at the nozzle exit, where the exit pressure, teperature ,velocity and mach number are denoted by `p_e,T_e.V_e and M_e`. The flow is locally subsonic in convergent section of the nozzle, sonic at the throat and supersonic at the divergent section. The sonic flow (M=1) at the throat means that the local velocity at this location is equal to the local speed of sound.We denote the pressure, temperature, velocity and area as p*, t*,v* (which is equal to V*=a* (speed of sound)), and A*.

We assume that at a given section , where the cross-sectional area is A the flow properties are unifrom across the section. hence the area of the section changes as a function of distance along the nozzle x and therfore in reality the flow field is 2-dimensional, we make the assumption that thee flow properties vary only with x, this is equal to assuming unifrom flow properties across any given cross section. Such flow is defined as quasi-one dimensional flow.

Step-1-Governing equations

For Non Conservation:

• Continutity Equation

Continuity: `((delrho^')/(delt^'))=-rho^'*(∂V^')/(∂x^')-rhoV'*(∂(ln A6'))/(∂x)-V'*(∂rho^')/(∂x')`

• Momentum Equation

`Momentum: (∂V')/(∂t')=-V'(∂V')/(∂x')-1/gamma((∂T')/(∂x')+(T')/(rho')(∂rho')/(∂x'))`

• Energy Equation

Energy: `(∂T')/(∂t')=-V'(∂T')/(∂x')-(gamma-1)T'[(∂V')/(∂x')+V'((∂(lnA'))/(∂x'))`

For Conservative

• Continuity Equation

`(∂(rho'A'))/(∂t')+(∂(rho'A'V'))/(∂x')=0`

• Momentum Equation

`(∂(rho'A'V'))/(∂t')+(∂[rho'A'V'^2+(1/gamma)p'A'])/(∂x')=(1/gamma)p'(∂A')/(∂x')`

• Energy Equation

`(∂[rho'((e')/(gamma-1)+(gamma/2)V'^2)A'])/(∂t')+(∂[rho'((e')/(gamma-1)+(gamma/2)V'^2)V'A'+p'A'V'])/(∂x')=0`

Step-2-Finite difference expressions using MacCranmack's explicit technique for the numerical solution

To implement a finite-differene solution, we divide the x axis along the nozzle into a number of discrete frid points,

The first grid point, labeled as point 1 , is assumed to be in reservior. The points are evenly distributed along the x axis, with `Delta x` denoting the spacing between grid points. The last point, namely, that at the nozzle exit, is denoted by N; we have a total number of N grid points distributed along the axis. Point i is simply an arbitrary grid point, with points i-1 and i+1 as the adjacent points. MacCormack's technique is predictor- corrector method. In the time marching approach, remember that we know the flowfield variables at time t, and we use the differene equations to solve explicitly for the variables at time `t+Deltat`.

Non-Conservasion from

First-we consider the predictor step, we set up the spatial derivatives as forward differences.

Continuity Equation:

`((∂rho)/(∂t))_i^(t)=-rho_i^(t)(V_(i+1)^(t)_-V_i^(t))/(Deltax)-rho_i^(t)V_i^t(lnA_(i+1)-lnA_i)/(Deltax)-V_i^t(rho_(i+1)^t-rho_i^t)/(Deltax)`

Momentum Equation:

`((delV)/(delV))_i^t=-V_i^t(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))_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))`

Using the above values we will obtain the predicted values of `rho,V and T`.

`rho_i^(-t+Deltat)=rho_i^t+((delrho)/(delt))_i^tDeltat`

`V_i^(-t+Deltat)=V_i^t+((delV)/(delt))_i^tDeltat`

`T_i^(-t+Deltat)=T_i^t+((delT)/(delt))_i^tDeltat`

But these values are only First order accurate.So we do the corrector step to make it Second order accurate.

In the corrector step we first obtain a predicted values of the time derivative at time `t+Deltat, ((delrho)/(delt))_(i,j)^(t+Deltat)` by substituting the predicted values of `rho,u and v` into the right side of the continuity equation, replacing the spacial derivates with rearward differences.

`((delrho)/(delt))_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)(lnA_i-lnA_(i-1))/(Deltax)-V_i^(-t+Deltat)(rho_i^(-t+Deltat)-rho_(i-1)^(-t+Deltat))/(Deltax)`

likewise momentum and energy equations will be,

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

`((delT)/(delt))_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-lnA_(i-1))/(Deltax))`

Using above equations, the average time derivatives are

`((delrho)/(dely))_(avg)=0.5[((delrho)/(delt))_i^t+((delrho)/(delt))_i^(t+Deltat)]`

`((delV)/(delt))_(avg)=0.5[((delv)/(delt))_i^t+((delV)/(delt))_t^(t+Deltat)]`

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

Then we can substitute the values into the below formula to find the corrected values of the flow-field variables at time `t+Deltat`

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

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

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

Conservation from

The equations for quasi-one-dimensional flow can be expressed in generic from.Let us define the elements of the solutions vector U, the flux vector F,and the source term j as follows

`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)V'A'+p'A'V'`

`J_2=1/gammap'(delA')/(delx')`

First we will find the source vectors with the above equations with intial condition, then with them we can find the flux vectors.

`F_1=U_2`

`F_2=U_2^2/U_1+(gamma-1)/gamma(U3-gamma/2 U_2^2/U_1)`

`F_3=gamma(U_2U_3)/U_1-(gamma(gamma-1))/2U_2^3/U_1^2`

`J_2=(gamm-1)/gamma(U_3-gamma/2U_2^2/U_1)(del(lnA'))/(delx')`

Predictor step

We set the spacial derivatives as forward differences.

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

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

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

We find the pridicted values of U1,U2 and U3

`(U_1)_i^(t'+Deltat')=(U_1)_i^(t')+((delU_1)/(delt'))_i^(t')Deltat'`

`(U_2)_i^(t'+Deltat')=(U_2)_i^(t')+((delU_2)/(delt'))_i^(t')Deltat'`

`(U_3)_i^(t'+Deltat')=(U_3)_i^(t')+((delU_3)/(delt'))_i^(t')Deltat'`

Using these values we will find the predicted values of primitive variables

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

`(T')_i^(t'+Deltat')=(gamma-1){(U_3)_i^(t'+Delta')/(U1)_i^(t'+Deltat')-gamma/2[(U_2)_i^(t'+Deltat')/(U_1)_i^(t'+Deltat')]^2}`

Now update the flux vector F1,F2,F3 with the new density, temperature and velocity.

Corrector Step:

The Predicted time derivatives of U1,U2 and U3 are obtained above and correct time derivatives can be find using rearward differences for the x derivatives.

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

`((delU_2)/(delt'))_i^(t'+Deltat')=-((F_2)_i-(F_2)_(i-1))/(Deltax')+1/gammarho'T'(delA')/(delx')`

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

Then the average time derivatives are formed as follows:

`((delU_1)/(delt))avg=1/2[((delU1)/(delt'))_i^(t')+((delU_1)/(delt'))_i^(t'_+Deltat')]`

`((delU_2)/(delt))avg=1/2[((delU2)/(delt'))_i^(t')+((delU_2)/(delt'))_i^(t'_+Deltat')]`

`((delU_3)/(delt))avg=1/2[((delU3)/(delt'))_i^(t')+((delU_3)/(delt'))_i^(t'_+Deltat')]`

So, the final corrected values of U1,U2,u3 will become

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

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

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

Finally the corrected values of the primitive variables are obtained by decoding U1,U2 and u3 obtained above.

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

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

`(T')_i^(t'+Deltat')=(gamma-1)((U_3)/(U1)-gamma/2V'^2)_i^(t'+Deltat')`

Step-3-Setting time step calculation, intial and boundary conditions

Time step can be calculated using the following formula:

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

where

C=Courant number

a=speed of sound

V=Velocity

`Deltax`=Distance between to consective points.

Non Conservative from

Nozzle shape:

`A=1+2.2(x-1.5)^2` `0lexle3`

initial conditions:

Intial conditions at t=0

`rho=1-0.3146x`

`T=1-0.2314x`

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

Boundary Condition at inlet:

`V_1=2V_2-V_3`

Boundary condition at outlet

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

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

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

Conservative form

Nozzle shape:

`A=1+2.2(x-1.5)^2` `0lexle3`

Intial condition:

for`0lex'le0.5`

`rho'=1`

`T'=1`

For`0.5lex'le1.5`

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

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

For`1.5lex'le3.5`

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

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

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

`U_1=rho'A'`

`U_2=rho'A'V'`

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

Boundary condition at inflow:

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

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

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

Boundary Condition at outflow:

`(U_1)_N=2(U_1)_(N-1)-(U_1)_(N-2)`

`(U_2)_N=2(U_2)_(N-1)-(U_2)_(N-2)`

`(U_3)_N=2(U_3)_(N-1)-(U_3)_(N-2)`

Code and Explanation:

Main Code:

clear all
close all
clc

n=61;
gamma=1.4;
x=linspace(0,3,n);
dx=x(2)-x(1);
nt=1400;

%Non Conservative form
[rho,v,t,p,M,M_t,rho_t,v_t,p_t,t_t,x,a]=non_conservative(n,gamma,x,dx,nt);

%Conservative form
[rho_c,v_c,t_c,p_c,M_c,rho_t_c,t_t_c,v_t_c,p_t_c,M_t_c,x,a]=conservative(n,gamma,x,dx,nt);

figure(3)
plot(x,M)
title_text=sprintf('Non conservative \n mach number vs Distance (No. of Timestep=%d',nt)
title(title_text)
xlabel('Distance')
ylabel('Mach number')

figure(4)
plot(x,p)
title_text=sprintf('Non conservative \n pressure vs Distance (No. of Timestep=%d',nt)
title(title_text)
xlabel('Distance')
ylabel('pressure')

figure(5)
plot(x,rho)
title_text=sprintf('Non conservative \n density vs Distance (No. of Timestep=%d',nt)
title(title_text)
xlabel('Distance')
ylabel('density')

figure()
plot(x,t)
title_text=sprintf('Non conservative \n Temperature vs Distance (No. of Timestep=%d',nt)
title(title_text)
xlabel('Distance')
ylabel('Temperature')

figure()
plot(x,M_c)
title_text=sprintf(' conservative \n mach number vs Distance (No. of Timestep=%d',nt)
title(title_text)
xlabel('Distance')
ylabel('Mach number')

figure()
plot(x,p_c)
title_text=sprintf('conservative \n pressure vs Distance (No. of Timestep=%d',nt)
title(title_text)
xlabel('Distance')
ylabel('pressure')

figure()
plot(x,rho_c)
title_text=sprintf('conservative \n density vs Distance (No. of Timestep=%d',nt)
title(title_text)
xlabel('Distance')
ylabel('density')

figure()
plot(x,t_c)
title_text=sprintf(' conservative \n Temperature vs Distance (No. of Timestep=%d',nt)
title(title_text)
xlabel('Distance')
ylabel('Temperature')


nt=(1:1400);

figure()
plot(nt,M_t)
title_text=sprintf('Non conservative \n mach number vs Distance at throat')
title(title_text)
xlabel('Distance')
ylabel('Mach number')

figure()
plot(nt,p_t)
title_text=sprintf('Non conservative \n pressure vs Distance at throat')
title(title_text)
xlabel('Distance')
ylabel('pressure')

figure()
plot(nt,t_t)
title_text=sprintf('Non conservative \n Temperature vs Distance at throat')
title(title_text)
xlabel('Distance')
ylabel('Temperature')

figure()
plot(nt,rho_t)
title_text=sprintf('Non conservative \n density vs Distance at throat')
title(title_text)
xlabel('Distance')
ylabel('density')

figure()
plot(nt,M_t_c)
title_text=sprintf('conservative \n mach number vs Distance at throat')
title(title_text)
xlabel('Distance')
ylabel('Mach number')


figure()
plot(nt,p_t_c)
title_text=sprintf('conservative \n pressure vs Distance at throat')
title(title_text)
xlabel('Distance')
ylabel('pressure')


figure()
plot(nt,t_t_c)
title_text=sprintf('conservative \n temperature vs Distance at throat')
title(title_text)
xlabel('Distance')
ylabel('Temperature')


figure()
plot(nt,rho_t_c)
title_text=sprintf('conservative \n density vs Distance at throat')
title(title_text)
xlabel('Distance')
ylabel('density')

%mass flow rate distribution inside the nozzle

figure()
m=rho.*a.*v;
m_c=rho_c.*a.*v_c;
plot(x,m,'LineWidth',2)
hold on
plot(x,m_c,'LineWidth',2)
title_text=sprintf('Mass Flow rate vs Distance');
title(title_text)
xlabel('Distance')
ylabel('Mass flow rate')
legend('non conservative from','Conservative from')

Explanation:

• n=31 : Sets the number of points to 31.
• gamma=1.4; Defines the adiabatic index (ratio of specific heats)
• x=linspace(0,3,n): Creates a linearly spaces vector 'x' from 0 to 3 with 31 points.
• Comptes the spacial grid spacing.
• Set the number of time step to 1400.
• Calls a function 'non_conservative' with the specified parameters and receives multiple output variables.
• Calls a function 'conservative' with the specified parameters and receives multiple output variables.
• Creates the figure and plotting for non-conservative and conservative forms of mach number, pressure, temperature, density against distance ('x').
• Set the title and labels for the plot.

Function Code For Non Conservative:

function [rho,v,t,p,M,M_t,rho_t,v_t,p_t,t_t,x,a]=non_conservative(n,gamma,x,dx,nt);

%calculate intial profiles
rho=1-0.3146*x;
t=1-0.2314*x;%t=temperature
v=(0.1+1.09*x).*t.^0.5;
%area
a=1+2.2*(x-1.5).^2;
%Time steps
delta_t=0.5*(dx./(sqrt(t)+v));
dt=min(delta_t);

%outer time loop
for k=1:nt

    %pridector method
    rho_old=rho;
    v_old=v;
    t_old=t;
    for i=2:n-1
       
        dvdx=(v(i+1)-v(i))/dx;
        drhodx=(rho(i+1)-rho(i))/dx;
        dtdx=(t(i+1)-t(i))/dx;
        dlogadx=(log(a(i+1))-log(a(i)))/dx;

        %continuity equation
        drhodt_p(i)=-rho(i)*dvdx-rho(i)*v(i)*dlogadx-v(i)*drhodx;
        %momentum equation
        dvdt_p(i)=-v(i)*dvdx-(1/gamma)*(dtdx+(t(i)/rho(i))*drhodx);
        %energy
        dtdt_p(i)=-v(i)*dtdx-(gamma-1)*t(i)*(dvdx+v(i)*dlogadx);
        %solution update
        rho(i)=rho(i)+drhodt_p(i)*dt;
        v(i)=v(i)+dvdt_p(i)*dt;
        t(i)=t(i)+dtdt_p(i)*dt;

    end
    for j=2:n-1
       
        dvdx=(v(j)-v(j-1))/dx;
        drhodx=(rho(j)-rho(j-1))/dx;
        dtdx=(t(j)-t(j-1))/dx;
        dlogadx=(log(a(j))-log(a(j-1)))/dx;

        %continuity equation
        drhodt_c(j)=-rho(j)*dvdx-rho(j)*v(j)*dlogadx-v(j)*drhodx;
        %momentum equation
        dvdt_c(j)=-v(j)*dvdx-(1/gamma)*(dtdx+(t(j)/rho(j))*drhodx);
        %energy
        dtdt_c(j)=-v(j)*dtdx-(gamma-1)*t(j)*(dvdx+v(j)*dlogadx);
        
    end
    %compute average time derivation
    dvdt=0.5*(dvdt_p+dvdt_c);
    drhodt=0.5*(drhodt_p+drhodt_c);
    dtdt=0.5*(dtdt_p+dtdt_c);
    %final solution update
    for i=2:n-1
        rho(i)=rho_old(i)+drhodt(i)*dt;
        v(i)=v_old(i)+dvdt(i)*dt;
        t(i)=t_old(i)+dtdt(i)*dt;
    end
    delta(k)=0.5*(dx/(a(i+1)+v(i+1)));
    %apply the boundary conditions
    %inlet
    v(1)=2*v(2)-v(3);
    %oulet
    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);
    p=rho.*t;
    M=(v./sqrt(t));
    %throat intial data
    throat=find(a==1);
    rho_t(k)=rho(throat);
    v_t(k)=v(throat);
    t_t(k)=t(throat);
    M_t(k)=M(throat);
    p_t(k)=p(throat);
    %mass
    mass=rho.*v.*a;
    figure(1)
    if k==50
        plot(x,mass,'LineWidth',2)
        hold on
    elseif k==100
        plot(x,mass,'LineWidth',2)
    elseif k==150
        plot(x,mass,'LineWidth',2)
    elseif k==200
        plot(x,mass,'LineWidth',2)  
    elseif k==700
        plot(x,mass,'LineWidth',2)
    elseif k==1400
        plot(x,mass,'LineWidth',2)    
    end
    title('Plot of mass flow rate at different time step Non conservative form')
    xlabel('Distance')
    ylabel('Non-dimensional mass flow rate')
    legend('50','100','150','200','700','1400')
end

• This line defines a function named 'non_conservative' that takes parameters 'n' ,'gamma' ,'x',dx', and 'nt'. It returns several variables: 'rho', 'v', 't', 'p', 'M', 'M_t', 'rho_t', v_t','p_t','t_t', 'x' amd 'a'.
• Intialize the initial profiles for density('rho') , temperature('t'), and velocity('v).
• Comptes the area 'a' based on the given formula.
• Computes the time step 'dt' based on the Courant-Friedrichs-Lewy(CFL) condition.
• Initiates the outer time loop that iterates from 1 to 'nt' (number of time steps).
• Stores the current values of density ('rho'), velocity('v'), and temperature ('t') in temporary variables for the predictor method.
• Intiates a loop over spatial points excluding the boundaries.
• Within the spatial loop, it calculates the preditor step for continuity, momentum, and energy equations, and updates the solution.
• Intiates another loop for the corrector method, calculating similar quantities for continuity, momentum, and energy equations.
• Computes the average timestep derivaties using both predictor and corrector methods.
• Updates the solution with the final values.
• After that we applying the boundary conditions, computing propertirs at the throat plotting mass flow rate at specfic time steps, and updating the variables for visualization.

Function Code For Conservative:

function [rho_c,v_c,t_c,p_c,M_c,rho_t_c,t_t_c,v_t_c,p_t_c,M_t_c,x,a]=conservative(n,gamma,x,dx,nt)

for i=1:length(x)
    if(x(i)<=0.5);
        rho(i)=1;
        t(i)=1;
    elseif(x(i) >0.5 && x(i)<=1.5)
            rho(i)=1-0.366*(x(i)-0.5);
            t(i)=1-0.167*(x(i)-0.5);
    else
        rho(i)=0.634-0.366*(x(i)-1.5);
        t(i)=0.833-0.3507*(x(i)-1.5);
    end
end
a=1+(2.2*(x-1.5).^2);  %area
v=0.59./(rho.*a);
delta_t=0.5*(dx./(sqrt(t)+v));
dt=min(delta_t);
u1=rho.*a
u2=rho.*a.*v
u3=rho.*((t/(gamma-1))+((gamma/2).*(v).^2)).*a;


for k=1:nt
    u1_old=u1;
    u2_old=u2;
    u3_old=u3;

    F1=u2;
    F2=((u2.^2)./u1)+((gamma-1)/gamma)*(u3-(gamma/2)*(((u2).^2)./u1));
    F3=(((gamma)*(u2.*u3))./u1)-(((gamma*(gamma-1))/2)*((u2.^3)./u1.^2));
    for j=2:n-1
        J2(j)=(1/gamma)*rho(j)*t(j)*((a(j+1)-a(j))/dx);

        du1dt_p(j)=-(F1(j+1)-F1(j))/dx;
        du2dt_p(j)=-((F2(j+1)-F2(j))/dx)+J2(j);
        du3dt_p(j)=-((F3(j+1)-F3(j))/dx);

        u1(j)=u1_old(j)+du1dt_p(j)*dt
        u2(j)=u2_old(j)+du2dt_p(j)*dt;
        u3(j)=u3_old(j)+du3dt_p(j)*dt;
    end

    rho=(u1./a);
    %v1=u2./u1;
    t=(gamma-1)*((u3./u1)-((gamma/2)*((u2./u1).^2)));

    F1=u2;
    F2=((u2.^2)./u1)+((gamma-1)/gamma)*(u3-(gamma/2)*(((u2).^2)./u1));
    F3=gamma*(u2.*u3./u1)-(((gamma*(gamma-1))/2)*((u2.^3)./u1.^2));

    for m=2:n-1
        J2(m)=(1/gamma)*rho(m)*t(m)*((a(m)-a(m-1))/dx);
        du1dt_c(m)=-(F1(m)-F1(m-1))/dx;
        du2dt_c(m)=-((F2(m)-F2(m-1))/dx)+J2(m);
        du3dt_c(m)=-(F3(m)-F3(m-1))/dx;
    end

    %average time derivative
    du1dt=0.5*(du1dt_p+du1dt_c);
    du2dt=0.5*(du2dt_p+du2dt_c);
    du3dt=0.5*(du3dt_p+du3dt_c);
    for i=2:n-1
        u1(i)=u1_old(i)+(du1dt(i)*dt);
        u2(i)=u2_old(i)+(du2dt(i)*dt);
        u3(i)=u3_old(i)+(du3dt(i)*dt);
    end
    %solution update
    rho_c=u1./a;
    v_c=(u2./u1);
    t_c=(gamma-1)*((u3./u1)-(gamma/2)*(u2./u1).^2);

    %defining inletflow Boundary conditions
    u1(1)=rho_c(1)*a(1);
    u2(1)=2*u2(2)-u2(3);
    v(1)=u2(1)/u1(1);
    u3(1)=u1(1)*(t_c(1)/(gamma-1)+(gamma/2)*v(1).^2);

    %defining outflow boundary conditions
    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);

    p_c=rho.*t;
    M_c=(v./sqrt(t));
    %throat intial data
    throat=find(a==1);
    rho_t_c(k)=rho_c(throat);
    v_t_c(k)=v_c(throat);
    t_t_c(k)=t_c(throat);
    M_t_c(k)=M_c(throat);
    p_t_c(k)=p_c(throat);
    %mass
    mass=rho_c.*v_c.*a;
    
    figure(2)
    if k==50
        plot(x,mass,'LineWidth',2)
         hold on
    elseif k==100
        plot(x,mass,'LineWidth',2)
    elseif k==150
        plot(x,mass,'LineWidth',2)
    elseif k==200
        plot(x,mass,'LineWidth',2)  
    elseif k==700
        plot(x,mass,'LineWidth',2)
    elseif k==1400
        plot(x,mass,'LineWidth',2)    
    end
    title('plot of mass flow rate at different time steps conservative from')
    xlabel('Distance')
    ylabel('Non-dimensional mass flow rate')
    legend('50','100','150','200','700','1400')
end

Explanation:

• This line defines a function named 'non_conservative' that takes parameters 'n' ,'gamma' ,'x',dx', and 'nt'. It returns several variables: 'rho_c', 'v_c', 't_c', 'p_c', 'M_c', 'M_t_c', 'rho_t_c', v_t_c','p_t_c','t_t_c', 'x' amd 'a'.
• Initialize density ('rho') and temperature('t') based on different condition for different regions of 'x'.
• Computes the srea ('a'), velocity('v') and initialize conservative u1, u2 and u3.
• Initiates the time loop for the conservative form, similar to the non-conservative case.
• Calculates the conservative variables (rho_c, v_c, t_c) from the conservative variables ('u1', 'u2', 'u3').
• Calculates pressure (p_c) and mach umber(M_c) from the conservative variables.
• Concludes the time loop and the function.

Results:

Conservative:

Grid independence test:

Flow variables at throat	anlytical solution	31grid points (n=16) non-conservative	31grid points (n=16) conservative	61grid points (n=16) non-conservative	61grid points (n=16) conservative
Density	0.634	0.6387	0.6473	0.6377	0.6376
Pressure	0.528	0.5342	0.5427	0.5327	0.522
Temperature	0.833	0.8365	0.8384	0.8354	0.8348
Mach number	1.000	0.994	0.9876	0.998	1.002

Observations:

• When we go through the graphs of distribution of primitive variables inside the nozzle, we can see that as we moves from the inlet to outlet the mach number goes on increaseing. At the inlet Mach number will be less than 1, at the throat the mach number will be equa to 1 and the at outlet the mach number will be greater than 1.
• Lokking at the pressure, density, temperature and mach number at the nozzle throat plotted vaerse the number of time step we can see that te largest changes take place at early times after which the final, steady state values is approched. At early times, the time derivatives are large and they oscillate in value. These oscillations are associated with various unstedy compression and expansion waves which propagate through the nozzle during the transient process. Later the time drivatives rapidly grow small over a finite number of time steps.
• While comparing the steady state variation of mass flow obtained with the non-conservation and conservation from of governing equation we can see that conservation from of governing equations is much more satisfactory than that obtained from non conservation from. The conservation from gives a distribution that is uch closer to being a constant. But the non conservation results have a sizeable variation, with some false oscillations at the both the inflow and outflow boundaries. Also the consevation from is much more closser to exact analytical solution.

Conclusion:

• Better mass flow distribution is done in conservation from.
• Non-conservation from is better for programming convenience because they are solved in term of primitive variables and their values are obtained directly. But in consrvation from equation are expressed in terms of flux variables and first solve solution vector then the primitive variables which makes things complex.
• For accuracy also non-conservative is better due to the direct calculation of primative variables. While in conservation from primitive variables are calculated indirectly from flux vaariables. Increased number of calculations will lead to more errors. This is the case whn we take primitive variables. But if we are taking Flux variables conservaton from is better.
• In con conservation from behavior of the residuals are good as compare to conservation from. The residuals decay faster in non conservation from when compared to conservation from. The amount by which the residual decay is often used as an index 'quality' of the numerical algoritm. So the non conservation from is better.