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

Objectives

• To write separate functions for conservative and non-conservative forms to and code to solve the 1D supersonic nozzle flow equations using the Macormack Method.
• To figure out the minimum number of cycles for which the simulation should be run in order for convergence
• To compare the normalized mass flow rate between the conservative forms and non-conservative forms - make suitable conclusions
• To find which method is faster?
• To perform grid dependence test.

In this case steady, isentropic flow through a convergent-divergent nozzle is studied. Inlet of nozzle is connected to the reservoir in which flow is subsonic, at throat flow is sonic and at the outlet flow is supersonic. The ultimate objective of the challenge is to write a MATLAB code which calculated mass flow rate through the nozzle using a conservative and non-conservative form of governing equations. At the later part of the challenge, the error in density, velocity, and temperature at each time step is discussed as well mass flow rate at different time steps for both forms of governing equations.The total length of the nozzle is taken as 3units in this case first the 1D nozzle is discretized into 31 grid points and gradually increased then grid dependency of mass flow rate is studied for different mesh parameters.

Nozzle -

A nozzle is a device designed to control the direction or characteristics of a fluid flow (especially to increase velocity) as it 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 the pressure of the stream that emerges from them. In a nozzle, the velocity of fluid increases at the expense of its pressure energy.

1D Supersonic nozzle flow:

The convergence(subsonic) section connected with reservoir, here mach number is less than one. At throat area, the thermodynamic properties will become critical property as Mach number equal to 1. The sonic flow at the throat means that the local velocity at this point is equal to the local speed of sound and at divergence section(supersonic) the flow expands isentropically to supersonic speed at the nozzle exit. The corresponding notations are clearly written in the figure itself.

The flow is locally subsonic in the convergence section of the nozzle, sonic at the throat and supersonic at the divergence section.

The governing continuity, momentum and energy equation for the quasi 1D steady isentropic flow can be expressed as

Continuity Equation - `rho_1A_1V_1=rho_2A_2V_2`

Momentum Equation - `p_1A_1+rho_1V_1^2+intpdA=p_2A_2+rho_2V_2^2`

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

The real physics of the problem lies on the integral forms of these equations. We will write the equations in their Non-conservative and conservative forms.

Rho - Density

V- Velocity

T - Temperature

t - Time

A - Area

MACORMACK METHOD:

This method is employed to solve the governing equations of quasi 1D supersonic nozzle flow rather than the differencing schemes as we have studied earlier. Maccormack method is an explicit finite difference technique which is second order accurate in both time and space.In computational fluid dynamics, the MacCormack method is a widely used discretization scheme for the numerical solution of hyperbolic partial differential equations. This method essentially uses a predictor step and a corrector step which will be depicted in the code below.

MACORMACK PROCEDURE:

According to the maccormack method, the value of variable at grid point i,j at time level t+dt is given as follows-

`v_(i,j)^(t+Deltat) = v_(i,j)^t + ((delv)/(delt))_(avg)*Deltat`

where the average time derivative is the mean value between time t and t+dt. The average time derivative is obtained from the predictor-corrector philosophy.

`((delv)/(delt))_(\"avg\") = 0.5[((delv)/(delt))_(\"predected\")+ ((delv)/(delt))_(\"corected\")]`

Predictor Step-

In the predictor step, the spatial terms on the right-hand side of the equation are replaced by forward differences. The RHS of the equation is known as all the values of variables are for time-step \'t\'. Next, the value of the variable on the LHS is calculated by employing the time marching solution and for time-step \'t+dt\' value of the variable at all grid points is obtained.

Consider the following equation for example-

`((delv)/(delt)) = ((delrho)/(delx))`

`((delv)/(delt)) _(predected)= ((rho_(i+1,j)^t-rho_(i,j)^t)/(Deltax) +(rho_(i+1,j)^t-rho_(i,j)^t)/(Deltay))`

`v_(i,j)^(t+Deltat) = v_(i,j)^t + ((delv)/(delt))_(predicted)*Deltat`

This is the predicted value of variable which will be used in the corrector step to evaluate the corrected value of variable.

Corrector Step-

In the corrector step, the spatial terms on the right-hand side are replaced by backward differences. The value of variables in RHS are taken which are obtained in the predicted step for time-step \'t+dt\'. Then the corrected time derivative is calculated.

`((delv)/(delt)) _(\"corrected\")= ((rho_(i,j)^(t+Deltat)-rho_(i-1,j)^(t+Deltat))/(Deltax) +(rho_(i,j)^(t+Deltat)-rho_(i-1,j)^(t+Deltat))/(Deltay))`

Now, the predicted and corrected value of time derivative is used to calculate the average time derivative which in turn gives the value of the variable at a given point for a particular time-step. The same steps are repeated until the values of variables are obtained for all grid points for a given number of time-steps

In this project the fluid properties are studied as it passes through the convergent-divergent nozzle. The fluid properties are studied for both Conservative and Non-Conservative types of flows. Though both conservative and non-conservative are mathematically equal, they are different in numerical descritisation techniques. One form can be converted to another exactly.

The governing continuity, momentum and energy equation for the quasi 1D steady isentropic flow can be expressed as

Continuity Equation - `rho_1A_1V_1=rho_2A_2V_2`

Momentum Equation - `p_1A_1+rho_1V_1^2+intpdA=p_2A_2+rho_2V_2^2`

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

The real physics of the problem lies on the integral forms of these equations. We will write the equations in their Non-conservative and conservative forms.

Rho - Density

V- Velocity

T - Temperature

t - Time

A - Area

The conservative and non-conservative form of governing equations are derived using the integral form of continuity, momentum and energy equation.The derivation is given in the book- Computational Fluid Dynamics by John Anderson Jr.

Non- Conservative Form

Contunity Equation: `(delrho)/(delt) = -rho(delu)/(delx) -rhou(del(lnA))/(delx)-u(delrho)/(delx)`

Momentum Equation: `(delu)/(delt) = -u(delu)/(delx)-1/gamma((delT)/(delx)+T/rho(delrho)/(delx))`

Energy Equation:`(delT)/(delt) = -u(delT)/(delx)-(gamma-1)T((delu)/(delx)+u(del(lnA))/(delx))`

All the variables are the non-dimensional values. These equations are solved using Maccormack method and the steady-state values of all variables i.e. Temperature, Velocity and density are obtained. From these variables, the Mach number and pressure is obtained.

`p = rhoT`

Mach number = `u/sqrtT`

Calculation of Time-step-

In order for the solution to be stable, the concept of Courant flow number was introduced in previous projects. For a linear Hyperbolic equation, the value of CFL<1 will give a stable solution.

For this project, CFL=0.5 and value of time-step is calculated as follows-

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

where, C is the Courant number, V is the local flow velocity at a point and a is the local speed of sound. Value of dt at all grid points is calculated and the minimum of them is chosen for the simulation.

Boundary Conditions-

At the inlet:

Velocity at the inlet is changing with time, while the density and temperature are constant with respect to time.

`u_1=2u_2+u_3`

Velocity at the inlet is calculated by the method of linear extrapolation. Values of velocity at node 2 and 3 are used for the same.

At the Outlet:

All the variables at the outlet are kept changing with respect to time and linear extrapolation is used to calculate the variables at the outlet.

`u_N = 2u_(N-1)+u_(N-2)`

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

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

Initial Conditions-

Area of the nozzle will remain constant throughout the time of simulation.

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

The flow-field variables vary with respect to time. At t=0, their values are given as-

`rho = 1+0.3146x`

`T = 1+0.2314x`

`u = (0.1+1.09x)T^0.5`

Conservative Form-

The governing equations are expressed in a generic form which consists of the solution vector U, the flux vector F and the source term J.

`(delU)/(delt)+(delF)/(delx) = J`

The equations for Quasi 1-D supersonic flow can be expressed in a similar generic form as follows-

`U_1 = rhoA` `F_1 = rho A u`

`U_2 = rho Au` `F_2 = rhoAu^2+1/gammapA` `J_2 = 1/gammap(delA)/(delx)`

`U_3 = rho(e/(gamma-1)+gamma/2u^2)` `F_3 = rho(e/(gamma-1)+gamma/2u^2)+1/gammapA`

The equations which are solved using Maccormack method are-

`(delU_1)/(delt) = (delF_1)/(delx)` `(delU_2)/(delt) = -(delF_2)/(delx)+J_2` `(delU_3)/(delt) = -(delF_3)/(delx)`

Further, the flux vectors and source terms are expressed in terms of solution vectors. Also, to find the field variables, they are written in terms of solution vectors.

`F_1 = U_1`

`F_2 = (U_2^2)/(U_1)+ (gamma-1)/gamma(U_3-gamma/2(U_2^2)/(U1))` `F_3 =gamma (U_2U_3)/(U_1) + (gamma(gamma-1))/2(U_2^3)/(U_1^2)` `J_2 = (gamma(gamma-1))/gamma(U_3-(gamma)/2U_2^2/U_1)del(lnA)/(delx)`

Also,

`rho = (U_1)/(A)` `u = (U_2)/(U_1)` `T = (gamma-1)((U_3)/U_1-gamma/2u^2)`

Calculation of time step is the same as that of the procedure followed for non-conservation form.

Boundary conditions-

At inlet:

U2 and U3 are changing with respect to time whereas U1 is kept as constant.

`U_(2(i=1))=U_(2(i=2))- U_(2(i=3)`

`U_3 = U_1(T/(gamma-1)+gamma/2u_2)`

At outlet:

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

Area of the nozzle is constant throughout the runtime of simulation and its value is same as that taken for non-conservative form.

Initial velocity is given as,

`u = (0.59)/(rhoA)`

At t=0, values of density and temperature throughout the nozzle are given as follows,

``rho\'=1.0` ``T\'=1.0` for `0lex\'le0.5`

`rho\' = 1.0-0.366(x\'-0.5)` `T\' = 1.0-0.167(x\'-0.5)` for `0.5 lex\'le1.5`

`rho\' = 0.634-0.387(x\'-1.5)` `T\' = 0.833-0.3507(x\'-1.5)` for `1.5 lex\'le3.5`

Main Code

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

Function Code for Non-Conservative Form

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

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

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

for k=1:nt
%copy of the variables
rho_old=rho;
T_old=T;
v_old=v;
    
%predictor step
    for j= 2:n-1
        dv_dx= (v(j+1)-v(j))/dx;
        drho_dx =(rho(j+1)-rho(j))/dx;
        dTdx = (T(j+1)-T(j))/dx;
        da_dx =(log(a(j+1))-log(a(j)))/dx;
        %continuity
        drhodt_p(j) = -rho(j)*dv_dx - rho(j)*v(j)*da_dx - v(j)*drho_dx;
        %momentum
        dvdt_p(j) = -v(j)*dv_dx - (1/gama)*(dTdx+(T(j)/rho(j))*drho_dx);
        %energy
        dTdt_p(j) = -v(j)*dTdx -(gama-1)*T(j)*(dv_dx+v(j)*da_dx);
       %solution updates
       v(j)= v(j)+dvdt_p(j)*dt;
       T(j)= T(j)+dTdt_p(j)*dt;
       rho(j)= rho(j)+drhodt_p(j)*dt;
    end
   
    % corrector Method
    for j= 2:n-1
        dv_dx= (v(j)-v(j-1))/dx;
        drho_dx =(rho(j)-rho(j-1))/dx;
        dTdx = (T(j)-T(j-1))/dx;
        da_dx =(log(a(j))-log(a(j-1)))/dx;
        %continuity
        drhodt_c(j) = -rho(j)*dv_dx-rho(j)*v(j)*da_dx-v(j)*drho_dx;
        %momentum
        dvdt_c(j) = -v(j)*dv_dx-(gama^-1)*(dTdx+(T(j)/rho(j))*drho_dx);
        %energy
        dTdt_c(j) = -v(j)*dTdx -(gama-1)*T(j)*(dv_dx+v(j)*da_dx);
    end
    
    %average
    dvdt_avg=0.5*(dvdt_c+dvdt_p);
    drhodt_avg = 0.5*(drhodt_c+drhodt_p);
    dTdt_avg = 0.5*(dTdt_p+dTdt_c);
    
    % Solution update
    for i =2:n-1
     v(i)= v_old(i)+dvdt_avg(i)*dt;
     T(i)= T_old(i)+dTdt_avg(i)*dt;
     rho(i)= rho_old(i)+drhodt_avg(i)*dt;
    end
    
    %Inlet condition by interpolation
    v(1)=2*v(2)-v(3);
    %outlet condition by interpolation
    v(n)=2*v(n-1)-v(n-2);
    T(n)=2*T(n-1)-T(n-2);
    rho(n)=2*rho(n-1)-rho(n-2);
     
   %properties
   mass = rho.*a.*v;
   M = v./((T).^0.5);
   p= rho.*T;
   
   %throat values
   rho_throat(k)=rho(th);
   v_throat(k)=v(th);
   T_throat(k)=T(th);
   M_throat(k)=M(th);
   p_throat(k)=p(th);
   mass_throat(k)=mass(th);
   
   %residues at throat.
   v_res(k)=dvdt_avg(th);
   rho_res(k)=drhodt_avg(th);
   T_res(k)=dTdt_avg(th);
   
%plotting the results
% plotting mass flow rate at different time steps
figure(1)
if k == 50
plot(x, mass, \'b\', \'linewidth\', 1.5)
hold on
elseif k == 100
plot(x, mass, \'r\', \'linewidth\', 1.5)
hold on
elseif k == 150
plot(x, mass, \'c\', \'linewidth\', 1.5)
hold on
elseif k == 200
plot(x, mass, \'g\', \'linewidth\', 1.5)
hold on
elseif k == 550
plot(x, mass, \'y\', \'linewidth\', 1.5)
hold on
elseif k == 700
plot(x, mass, \'k\', \'linewidth\', 1.5)
hold on
elseif k == 1400
plot(x, mass,\'m\', \'linewidth\', 1.5)
grid on
ylabel(\'mass flow rate\');
xlabel(\'nozzle length (x)\');
title(\'mass flow rate at different time steps Non conservative approach\');
legend(\'50 time steps\',\'100 time steps\',\'150 time steps\',\'200 time steps\',\'550 time steps\',\'700 time steps\',\'1400 time steps\');
end
end

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

Function Code for Conservative Form

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

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


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

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

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


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

rho=rho_initial;
T= T_initial;


%time loop

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


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

Results for Grid size n = 31

1. Variation of Flow Properties along the Length of Nozzle

Observations from graphs

It can be observed from the above graph that both conservative and non conservative approach are in close agreement at any point in a nozzle.

2. Variation of flow properties with time step at throat

Observations from graph

a. From the graph we can see that both the conservative and non-conservative tends to steady state after 600th time step.

b. Non-conservative solution is in more close agreement with analytical solution than the conservative solution.

c. The plot that the stability is achieved by the non-conservative form earlier than the conservative form but the amplitude of the initial variations is also higher in the non-conservation form. Also the convergence criteria has been kept equal for both the cases.

3. Mass flow rate variation at different time steps-Conservative approach

Observations from graphs

a.It can be seen that the curve for 550,700 and 1400 timesteps are close and lie in the same band

b.The fluctuations are less as we assumed to have constant mass flow rate intially.

4. Mass flow rate variation at different time steps non-conservative approach

Observations from graphs

a.Even here we can see that the curve for 550,700 and 1400 time step lie in the same band are close.

b.The fluctuations are high as the mass flow rate is derived from other parameter.

5. Comparision of mass flow rate for Conservative, non-conservative and analytical approaches

Observations from graphs

i) Conservative curve is smooth and lies close to the analytical value since mass flow rate is directly calculated and the assumption in the begining is to have to a constant mass flow rate.

ii) Non-conservative approaches close to the analytical value only at the throat.

iii) The profile of the conservative form is much closer to the exact solution and on the other hand the variation in the mass flow rate of the non-conservative form is larger and also far from the exact analytical solution.

6. Residuals graphs

Observations from graphs

all the resudals from graph for conservative and not conservative is is converging at time step 800 and 600 respectively

Grid Dependency:-

The solution is said to be grid indepenent if the solution has no effect even if the number of grid points are increased once a stable solution is oobtained

In our case let us double the no. of grid points to 61 and acoordingly the dt should be changed. This can be obtained from the timestep calculation from courant number for each grid and taking the least.

So for n=61 we should have 2800 time steps, therefore for grid-independency the solution for n=31 must repeat at exactly the double time step.

Results for Grid Size n=61

1. Variation of Flow Properties along the Length of Nozzle

2. Variation of flow properties with time step at throat

3. Mass flow rate variation at different time steps-Conservative approach

4. Mass flow rate variation at different time steps non-conservative approach

5. Comparision of mass flow rate for Conservative, non-conservative and analytical approaches

6. Residuals graphs

Grid Dependency Test

For n=31

For n=61