1D Super-sonic nozzle flow simulation using Macormack Method in Octave / MATLAB

Aim:-> To perform numerical flow simulation of quasi 1D subsonic-supersonic nozzle using MacCormack technique with conservative and non conservative methods.

Theory:-

MacCormack technique:- MacCormack method is an explicit finite difference technique which is second order accurate in both time and space. The method is applied for discretizing hyperbolic partial differential equation. As in MacCormack method there are two steps i.e. predictor step with forward differencing and corrector step with backward differencing, the solution obtained by taking average of predictor and corrector step gives second order accurate. The procedure for this technique is as follows.

1. calculate initial values for equation
2. To calculate predictor step values by using forward differencing for space terms.
3. Update values on RH side of equation
4. To calculate corrector step values by using backward differencing for space term.
5. Take average of corrector step values and predictor step values. `("Corrector value"+"Predictor value")/2`
6. Update RH side values by taking previous timestep values
7. Repeat above all steps for successive timestep until solution gets converged.

There are two methods to solve almost every CFD problem numerically i.e. by using conservative equations and non conservative equations. From theoretical and mathematical view both forms are equivalent and can be converted into one another. For the engineering problems, user can choose any of the above flow method, each method possess its own advantages and disadvantages.

~ Conservative form

Consider an infinitesimal small fluid element inside nozzle, and the fluid moving through that element, such type of condition is conservative form and the equation obtained in such cases are called as conservative form of governing equations.

~ Non conservative form

Consider an infinitesimal small fluid element inside nozzle, the fluid same fluid element moving along streamline, such type of condition is non conservative form and the equation obtained in such cases are called as non conservative form of governing equations.

Assumptions for problem:-

1. Consider steady isentropic flow through a convergent divergent nozzle.
2. Flow is compressible.
3. Ignore friction and heat transfer at surface.
4. The pressure ratio is maintained to avoid shockwaves.
5. The flow is quasi 1D i.e. velocity varies along X-axis only (flow properties are function of x only).
6. Length of nozzle = 3.
7. The area of nozzle varies according to equation, `A = 1+2.2(x-1.5)^2 " "0<=x<=3`

The above figure shows schematic diagram of convergent divergent nozzle.

The flow at inlet of nozzle come from reservoir. For the convergent section flow is subsonic i.e. Mach number < 1, divergent section flow is supersonic i.e. Mach number > 1 and at throat flow is sonic i.e. Mach number = 1.

Where, `P_0 =` Pressure at inlet of nozzle.

`T_0 =` Temperature at inlet of nozzle.

`rho_0 =` Density at inlet of nozzle.

A = Area, V= Velocity, M = Mach number.

`P^star =` Pressure at throat.

`T^star =` Temperature at throat.

`V^star =` `a^star =` Sound velocity.

`P_e =` Pressure at exit of nozzle.

`T_e =` Temperature at exit of nozzle.

`V_e =` Velocity at exit of nozzle.

Governing equations:-

As we considering the flow through is only function of X and not Y hence we must apply integral form of Navier-Stokes equation to a control volume consistent with the quasi 1D assumption. Hence by referring below figure we obtain PDE form of governing equations from integral forms.

The shaded region in below figure shows control volume i.e. slice of nozzle flow with infinitesimal thickness of dx.

For the Quasi 1D flow governing equations are as follows-

~Continuity equation

`rho_1 A_1 V_1 = rho_2 A_2 V_2`

~Momentum equation

`P_1 A_1 + rho_1 A_1V_1^2 + int_(A_1)^(A2) P.da = P_2A_2 + rho_2 A_2 V_2^2`

~Energy equation

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

Where, `h=c_pT`

1. Analytical values:-

The mach number can be calculated from below equation.

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

Hence from above Mach number we can able to calculate pressure, density and temperature values easily.

-> `(P/P_0)=(1+(gamma-1)/2 M^2)^((-gamma)/(gamma-1))`

-> `(rho/rho_0)=(1+(gamma-1)/2 M^2)^((-1)/(gamma-1))`

->`(T/T_0)=(1+(gamma-1)/2 M^2)^-1`

By above equations we get exact analytical values for nondimensional pressure ratio, nondimensional density ratio and nondimensional temperature ratio.

The following figures consisting graphs which are plotted by using analytical values v/s nozzle length.

2. Non conservation form:-

~Continuity equation

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

~Momentum equation

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

~Energy equation

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

Where, dash superscript denotes non dimensional variables.

`rho'=rho/(rho_0)`, `A'=A/(A^star)`, `T'=T/(T_0)`, `V'=V/(a_0)`

`a_0=(gamma. R. T_0)^0.5`, `t'=t/(L/(a_0))`, `x'=x/L`

For computation optimization purpose we have converted Navier-Stokes equation in non dimensional parameters.

3. Conservation form:-

~Continuity equation

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

~Momentum equation

`(del(rho' A' V'))/(delt') + (del[rho'A'V'^2 + (1/gamma)P'A'])/(delx') =(1/gamma)P'(delA')/(delx')`

~Energy equation

`(del[rho'((e')/(gamma-1)+gamma/2V'^2)A'])/(dt') + (del[rho'((e')/(gamma-1)+gamma/2V'^2)A'V' + P'A'V'])/(dx')=0`

Let us convert primitive variables in to solution vector 'U', flux vector 'F' and source term 'J'.

Hence we get governing equation as-

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

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

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

Setup of nozzle flow problem:-

Mesh preperation:-

The 1D nozzle domain i.e. X-axis along nozzle can be divided in number of discrete grid points. The last grid which is at nozzle exit can be denote by 'N'. Hence the domain is divided into 'N-1' elements. Here `i` is arbitrary grid point in domain the grid points `i+1` and `i+1` are adjacent grid points. The below figure gives clear understanding of mesh preparation. Where `Deltax` denotes mesh element size.

Discretization:-

It can be done by setting a finite difference expressions using MacCormak's explicit technique.

1. Non conservative form:-

Predictor step:[Forward differencing]

~Continuity equation

~Momentum equation

~Energy equation

Now to obtain predicted values of `'rho'," "'V'" & " 'T'`

Corrector values:[Backward differencing]

~Continuity equation

~Momentum equation

~Energy equation

The average time derivative values can be calculated by taking average of corrector time derivative and predictor time derivative values as follows,

The final value of `'rho'," "'V'" & " 'T'` can be obtained from following equations.

Initial condition:-

`rho=1-0.3146x`

`T=1-0.2314x`

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

Boundary conditions:-

->Subsonic inflow condition:-

`V_1 = 2V_2 -V_3`

->Supersonic outflow condition:-

`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. Conservation form:-

Predictor step:

`(dU_1)/(dt) = (F_(1(i+1))-F_(1(i)))/dx`

`(dU_2)/(dt) = (F_(2(i+1))-F_(2(i)))/dx + (J_(2(i+1))-J_(2(i)))/dx`

`(dU_3)/(dt) = (F_(3(i+1))-F_(3(i)))/dx`

For corrector method updating values of U1, U2 & U3.

Corrector step:

`(dU_1)/(dt) = (F_(1(i))-F_(1(i-1)))/dx`

`(dU_2)/(dt) = (F_(2(i))-F_(2(i-1)))/dx + (J_(2(i))-J_(2(i-1)))/dx`

`(dU_3)/(dt) = (F_(3(i))-F_(3(i-1)))/dx`

The average time derivative values can be calculated by taking average of corrector time derivative and predictor time derivative values as follows,

Hence the final values of U1, U2 & U3 are obtain from following equations.

Initial condition:-

`"For "0<=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)),(T',=,0.833-0.3507(x'-1.5)):}`

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

Boundary conditions:-

->Subsonic inflow condition:-

`U_(1(i=1)) = A'_((i=1))`

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

`U_(3(i=1)) = U_(1(i=1))((T'_(i=1))/(gamma -1) + (gamma/2)V'_(i=1)^2)`

->Supersonic outflow condition:-

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

Timestep calculation `(Deltat)`:

By following relation we can able to calculate timestep value.

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

Where, C is courant number and by going through some literature review the courant number should be `<=` 1.

`Deltax` is a grid spacing or mesh element size.

`a_0` is a speed of sound.

So by using some relations the formula for timestep becomes,

`Deltat'_(i)= C(dx)/(T'_(i)^(0.5) + V'_(i))`

As there is no certainty for accuracy of local time stepping so for selecting timestep we can use minimum value of timestep in grids.

`Deltat="minimum"(Deltat_1^t,Deltat_2^t,Deltat_3^t,.....,Deltat_i^t,......,Deltat_(N-1)^t,Deltat_N^t)`

Although this consistent time marching method may take more time than local time marching method to give converged solution but consistent time marching method gives us physically meaningful transient variation which frequently are of intrinsic value by themselves.

~Non conservation function code:

%-----Non conservative function-----%

function [dt,rho,vel,temp,m_nc,Pressure_nc,mach_nc,nt,rho_t,Pressure_t,m_t,temp_t,vel_t]=N_C(n,x,dx,area,th,c,gamma)

  %Initial conditions
  rho = 1-0.3146*x;
  temp = 1-0.2314*x;
  vel = (0.1+1.09*x).*temp.^(0.5);

  %Time step calculation
  for i=1:n
    dt(i)= c*dx/((temp(i))^(0.5)+vel(i));
  endfor
  dt=min(dt);
  nt = 0:1:1400;
  
  %Time loop
  for i=1:length(nt)
    
    vel_old=vel;
    rho_old=rho;
    temp_old=temp;
    
    %Space (grid) loop
    %predictor method
    for j=2:n-1
      dv_dx = (vel(j+1)-vel(j))/dx;
      dlnA_dx = (log(area(j+1))-log(area(j)))/dx;
      dtemp_dx = (temp(j+1)-temp(j))/dx;
      drho_dx = (rho(j+1)-rho(j))/dx;
      
      %Continuity equation
      drho_dt_p(j) = -(rho(j)*dv_dx) - (rho(j)*vel(j)*dlnA_dx) - (vel(j)*drho_dx);
      
      %Momentum equation
      dv_dt_p(j) = -(vel(j)*dv_dx) - (1/gamma)*(dtemp_dx+(temp(j)/rho(j))*drho_dx);
      
      %Energy equation
      dtemp_dt_p(j) = -(vel(j)*dtemp_dx)-(gamma-1)*temp(j)*(dv_dx+vel(j)*dlnA_dx);
    
      %Solution update
      rho(j) = rho(j) + drho_dt_p(j)*dt;
      vel(j) = vel(j) + dv_dt_p(j)*dt;
      temp(j)= temp(j)+ dtemp_dt_p(j)*dt;
      
    endfor
    
    %corrector method
    for j=2:n-1
      dv_dx = (vel(j)-vel(j-1))/dx;
      dlnA_dx = (log(area(j))-log(area(j-1)))/dx;
      dtemp_dx = (temp(j)-temp(j-1))/dx;
      drho_dx = (rho(j)-rho(j-1))/dx;

      %Continuity equation
      drho_dt_c(j) = -(rho(j)*dv_dx) - (rho(j)*vel(j)*dlnA_dx) - (vel(j)*drho_dx);
      
      %Momentum equation
      dv_dt_c(j) = -(vel(j)*dv_dx) - (1/gamma)*(dtemp_dx+(temp(j)/rho(j))*drho_dx);
      
      %Energy equation
      dtemp_dt_c(j) = -(vel(j)*dtemp_dx)-(gamma-1)*temp(j)*(dv_dx+vel(j)*dlnA_dx);
    
    endfor

    %Finding average values
    drho_dt_nc = 0.5*(drho_dt_p + drho_dt_c);
    dv_dt_nc = 0.5*(dv_dt_p + dv_dt_c);
    dtemp_dt_nc = 0.5*(dtemp_dt_p + dtemp_dt_c);

    %Updating old values
    for k=2:n-1
      rho(k) = rho_old(k) + drho_dt_nc(k)*dt;  
      vel(k) = vel_old(k) + dv_dt_nc(k)*dt;
      temp(k) = temp_old(k) + dtemp_dt_nc(k)*dt;
    endfor
    
    %Boundary condition
    %Inflow boundary condition
    vel(1) = 2*vel(2)-vel(3);

    %Outflow boundary condition
    rho(n) = 2*rho(n-1)-rho(n-2);  
    vel(n) = 2*vel(n-1)-vel(n-2);
    temp(n) = 2*temp(n-1)-temp(n-2);
     
    %mass flow rate
    m_nc = rho.*area.*vel;
    
    %Pressure
    Pressure_nc = rho.*(temp);
    
    %Mach number
    mach_nc = vel./(temp).^0.5;
    
    figure(1);
    if i==50
      plot(x,m_nc,'linewidth',1.75);
      hold on;
      xlabel('Non dimensional length of nozzle');
      ylabel('Non dimensional mass flow rate');
    elseif i==100 
      plot(x,m_nc,'linewidth',1.75);
    elseif i==150
      plot(x,m_nc,'linewidth',1.75);
    elseif i==300
     plot(x,m_nc,'linewidth',1.75);
    elseif i==700
      plot(x,m_nc,'linewidth',1.75);
    elseif i==1400
      plot(x,m_nc,'linewidth',1.75);
    title('Mass flow rate for non conservation form at diffrent timesteps','fontsize',14);
    legend('50 timestep','100 timestep','150 timestep','300 timestep','700 timestep','1400 timestep');
    end
    

    %Values at throat
    rho_t(i)=rho(th);
    temp_t(i)=temp(th);
    vel_t(i)=vel(th);
    Pressure_t(i)=Pressure_nc(th);
    m_t(i)=m_nc(th);

  endfor

    figure(2)
    subplot(5,1,1);
    plot(nt,rho_t,'linewidth',1.75,'r');
    hold on;
    xlabel('No of time step');
    ylabel('rho / rho_o'); 
    title('Value at throat for variable timesteps (Non conservation form)','fontsize',14);
    subplot(5,1,2);
    plot(nt,temp_t,'linewidth',1.75,'b');
    hold on;
    xlabel('No of time step');
    ylabel('T / T_o');
    subplot(5,1,3);
    plot(nt,vel_t,'linewidth',1.75,'g');
    hold on;
    xlabel('No of time step');
    ylabel('Velocity');
    subplot(5,1,4);
    plot(nt,Pressure_t,'linewidth',1.75,'c');
    hold on;
    xlabel('No of time step');
    ylabel('P / P_o');
    subplot(5,1,5);
    plot(nt,m_t,'linewidth',1.75,'m');
    hold on;
    xlabel('No of time step');
    ylabel('Massflow rate');
    
endfunction

All of the above plots shows result for non conservative form.

~Conservation function code:

%-----Conservative function-----%

function [dt,rho_c,vel_c,temp_c,m_c,p_c,mach_c,nt,rho_tc,p_tc,m_tc,temp_tc,vel_tc]=C_f(x,n,dx,area,th,c,gamma)

  %Initial conditions
  for i=1:length(x)
    if x(i)>=0 && x(i)<0.5
      rho_c(i) = 1;
      temp_c(i) = 1;
    elseif x(i)>=0.5 && x(i)<1.5
      rho_c(i) = 1-0.366*(x(i)-0.5);
      temp_c(i) = 1-0.167*(x(i)-0.5);
    elseif x(i)>=1.5 && x(i)<3.5
      rho_c(i) = 0.634-0.3879*(x(i)-1.5);
      temp_c(i) = 0.833-0.3507*(x(i)-1.5);
    endif
  endfor

  vel_c = (0.59)./(rho_c.*area);
  u1 = rho_c.*area;
  u2 = rho_c.*area.*vel_c;
  u3 = u1.*((temp_c./(gamma-1))+(gamma/2).*vel_c.^2);
  
  %Time step calculation
  for i=1:n
    dt(i)= (c*dx)/((temp_c(i)^(0.5)) + vel_c(i));
  endfor
  dt=min(dt);
  nt = 0:1:1400;
    
  for i=1:length(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)));
    
    %predictor step
    for j=2:n-1
      
      %Continuity equation
      dF1_dx_p(j) = (F1(j+1)-F1(j))/dx;
      du1_dt_p(j) = -dF1_dx_p(j);

      J2_p(j) = ((gamma-1)/gamma)*(u3(j)-((gamma/2)*(((u2(j))^2)/u1(j))))*((log(area(j+1))-log(area(j)))/dx);
      
      %Momentum equation
      dF2_dx_p(j) = (F2(j+1)-F2(j))/dx;
      du2_dt_p(j) = -dF2_dx_p(j) + J2_p(j);

      %Energy equation
      dF3_dx_p(j) = (F3(j+1)-F3(j))/dx;
      du3_dt_p(j) = -dF3_dx_p(j);
      
      %updating old values
      u1(j) = u1(j) + du1_dt_p(j)*dt;
      u2(j) = u2(j) + du2_dt_p(j)*dt;
      u3(j) = u3(j) + du3_dt_p(j)*dt;
    endfor 
    
    rho_c = u1./area;
    vel_c = u2./u1;
    temp_c = (gamma-1)*((u3./u1)-(gamma/2)*vel_c.^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)));

    %corrector step
    for j=2:n-1
  
      %Continuity equation
      dF1_dx_c(j) = (F1(j)-F1(j-1))/dx;
      du1_dt_c(j) = -dF1_dx_c(j);

      J2_c(j) = ((gamma-1)/gamma)*(u3(j)-(gamma/2)*((u2(j))^2/u1(j)))*((log(area(j))-log(area(j-1)))/dx);
      
      %Momentum equation
      dF2_dx_c(j) = (F2(j)-F2(j-1))/dx;
      du2_dt_c(j) = -dF2_dx_c(j) + J2_c(j);

      %Energy equation
      dF3_dx_c(j) = (F3(j)-F3(j-1))/dx;
      du3_dt_c(j) = -dF3_dx_c(j);
    endfor
    
    %Average values
    du1_dt = 0.5*(du1_dt_c + du1_dt_p);
    du2_dt = 0.5*(du2_dt_c + du2_dt_p);
    du3_dt = 0.5*(du3_dt_c + du3_dt_p);
    
    %Updating old values
    for k=2:n-1
      u1(k) = u1_old(k) + du1_dt(k)*dt;
      u2(k) = u2_old(k) + du2_dt(k)*dt;
      u3(k) = u3_old(k) + du3_dt(k)*dt;
    endfor
    
    %Boundary conditions for problem 
    %Inflow boundary condition
    u1(1)=rho_c(1)*area(1);
    u2(1)=2*u2(2) - u2(3);
    u3(1)=u1(1)*[((temp_c(1))/(gamma-1))+((gamma/2)*(u2(1)/u1(1))^2)];
    
    %Outflow boundary condition
    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);

    rho_c = u1./area;
    vel_c = u2./u1;
    temp_c = (gamma-1).*((u3./u1)-((gamma/2).*vel_c.^2));
    p_c = rho_c.*temp_c;
    mach_c = vel_c./((temp_c).^0.5);
    m_c = rho_c.*area.*vel_c;
    
    figure(4);
    if i==50
      plot(x,m_c,'linewidth',1.75);
      hold on;
      xlabel('Non dimensional length of nozzle');
      ylabel('Non dimensional mass flow rate');
    elseif i==100 
      plot(x,m_c,'linewidth',1.75);
    elseif i==150
      plot(x,m_c,'linewidth',1.75);
    elseif i==300
     plot(x,m_c,'linewidth',1.75);
    elseif i==700
      plot(x,m_c,'linewidth',1.75);
    elseif i==1400
      plot(x,m_c,'linewidth',1.75);
      title('Mass flow rate for conservation form at diffrent timesteps','fontsize',14);
      legend('50 timestep','100 timestep','150 timestep','300 timestep','700 timestep','1400 timestep');
    end
    
    %Values at throat
    rho_tc(i) = rho_c(th);
    temp_tc(i) = temp_c(th);
    vel_tc(i) = vel_c(th);
    p_tc(i) = p_c(th);
    m_tc(i) = m_c(th);

  endfor

    figure(5);
    subplot(5,1,1);
    plot(nt,rho_tc,'linewidth',1.75,'r');
    hold on;
    xlabel('No of time step');
    ylabel('rho/rho_o');
    title('Value at throat for variable timesteps (Conservation form)','fontsize',14)
    subplot(5,1,2);
    plot(nt,temp_tc,'linewidth',1.75,'b');
    hold on;
    xlabel('No of time step');
    ylabel('T / T_o');
    subplot(5,1,3);
    plot(nt,vel_tc,'linewidth',1.75,'g');
    hold on;
    xlabel('No of time step');
    ylabel('Velocity');
    subplot(5,1,4);
    plot(nt,p_tc,'linewidth',1.75,'c');
    hold on;
    xlabel('No of time step');
    ylabel('P / P_o');
    subplot(5,1,5);
    plot(nt,m_tc,'linewidth',1.75,'m');
    hold on;
    xlabel('No of time step');
    ylabel('Massflow rate');
          
endfunction

All of the above plots shows result for conservative form.

Main octave code:

clear all
close all
clc

%Number of grids
n=31;
%Mesh creation
x=linspace(0,3,n);
dx=x(2)-x(1);

%Cross section area (function of x) 
area = 1+2.2*(x-1.5).^2;
%Throat grid calculation
th = find(area==1);
%Courant number
c=0.5;
gamma = 1.4;
gp = gamma+1;
gm = gamma-1;

%radius of nozzle at varius cross section
r=sqrt(area/pi);

%Analytical mach number calculation
for i=1:n
  eqn = @(M) (1/M^2)*((2/gp)+((gm*M^2)/gp))^(gp/gm) - area(i)^2;
  if ith
    x0 = [1+1e-3 5];
    M(i) = fzero(eqn,x0);
  end
end

%Analytical mach number
mach_ana = M;

%Analytical pressure
Pressure_ana=[1+((gamma-1)/2).*mach_ana.^2].^(-gamma/(gamma-1));
%Analytical density
rho_ana=[1+((gamma-1)/2).*mach_ana.^2].^(-1/(gamma-1));
%Analytical temperature
temp_ana=[1+((gamma-1)/2).*mach_ana.^2].^(-1);

%-----Non conservative function-----%
tic;
  [dt,rho,vel,temp,m_nc,Pressure_nc,mach_nc,nt,rho_t,Pressure_t,m_t,temp_t,vel_t]=N_C(n,x,dx,area,th,c,gamma);
ellapsed_nc=toc;

%error for NC form
error_p_nc = abs(Pressure_ana-Pressure_nc);
error_rho_nc = abs(rho_ana-rho);
error_temp_nc = abs(temp_ana-temp);
error_mach_nc = abs(mach_ana-mach_nc);

figure(3);
subplot(4,1,1);
plot(x,Pressure_nc,'linewidth',1.75);
xlabel('Nozzel length');
ylabel('P / P_o');
hold on;
plot(x,Pressure_ana,'linewidth',1.75,'*');
legend('Numerical','Analytical');
title('Non conservation form','fontsize',14)
subplot(4,1,2);
plot(x,rho,'linewidth',1.75);
xlabel('Nozzel length');
ylabel('rho / rho_o');
hold on;
plot(x,rho_ana,'linewidth',1.75,'*');
legend('Numerical','Analytical');
subplot(4,1,3);
plot(x,temp,'linewidth',1.75);
xlabel('Nozzel length');
ylabel('T / T_o');
hold on;
plot(x,temp_ana,'linewidth',1.75,'*');
legend('Numerical','Analytical');
subplot(4,1,4);
plot(x,mach_nc,'linewidth',1.75);
xlabel('Nozzel length');
ylabel('Mach number');
hold on;
plot(x,mach_ana,'linewidth',1.75,'*');
legend('Numerical','Analytical');

%-----Conservative function-----%
tic;
  [dt,rho_c,vel_c,temp_c,m_c,p_c,mach_c,nt,rho_tc,p_tc,m_tc,temp_tc,vel_tc]=C_f(x,n,dx,area,th,c,gamma);
ellapsed_c=toc;

%error for Conservation form
error_p_c = abs(Pressure_ana-p_c);
error_rho_c = abs(rho_ana-rho_c);
error_temp_c = abs(temp_ana-temp_c);
error_mach_c = abs(mach_ana-mach_c);

figure(6);
subplot(4,1,1);
plot(x,p_c,'linewidth',1.75);
xlabel('Nozzel length');
ylabel('P / P_o');
hold on;
plot(x,Pressure_ana,'linewidth',1.75,'*');
legend('Numerical','Analytical');
title('Conservation form','fontsize',14)
subplot(4,1,2);
plot(x,rho_c,'linewidth',1.75);
xlabel('Nozzel length');
ylabel('rho / rho_o');
hold on;
plot(x,rho_ana,'linewidth',1.75,'*');
legend('Numerical','Analytical');
subplot(4,1,3);
plot(x,temp_c,'linewidth',1.75);
xlabel('Nozzel length');
ylabel('T / T_o');
hold on;
plot(x,temp_ana,'linewidth',1.75,'*');
legend('Numerical','Analytical');
subplot(4,1,4);
plot(x,mach_c,'linewidth',1.75);
xlabel('Nozzel length');
ylabel('Mach number');
hold on;
plot(x,mach_ana,'linewidth',1.75,'*');
legend('Numerical','Analytical');

%Plot for comparing errors
figure(7);
subplot(4,1,1);
plot(x,error_p_c,'linewidth',1.75);
hold on;
plot(x,error_p_nc,'linewidth',1.75);
xlabel('Nozzel length');
ylabel('pressure error');
legend('Conservation form','Non conservation form');
title('Comparison of error values','fontsize',14);
subplot(4,1,2);
plot(x,error_rho_c,'linewidth',1.75);
hold on;
plot(x,error_rho_nc,'linewidth',1.75);
xlabel('Nozzel length');
ylabel('rho error');
legend('Conservation form','Non conservation form');
subplot(4,1,3);
plot(x,error_temp_c,'linewidth',1.75);
hold on;
plot(x,error_temp_nc,'linewidth',1.75);
xlabel('Nozzel length');
ylabel('Temp error');
legend('Conservation form','Non conservation form');
subplot(4,1,4);
plot(x,error_mach_c,'linewidth',1.75);
hold on;
plot(x,error_mach_nc,'linewidth',1.75);
xlabel('Nozzel length');
ylabel('Mach error');
legend('Conservation form','Non conservation form');

%Plot for comparing non dimensional flow rate
figure(8);
plot(x,m_c,'linewidth',1.75);
hold on;
plot(x,m_nc,'linewidth',1.75);
xlabel('Non dimensional nozzle length');
ylabel('Masss flow rate');
legend('Conservation form','Non-conservation form');
title('Comaprison of mass flow rate','fontsize',14);

%Plot for comparing ellapsed time
figure(9);
bar(1,ellapsed_nc,'r');
hold on;
bar(3, ellapsed_c,'b');
ylabel('Ellapsed time');
title('Comparison for ellapsed time','fontsize',14)
legend('Non-conservative form','Conservative form');

R = max(r);
y =linspace(-R, R, n);
Y = meshgrid(y);
Y = -1*Y';
Mach = zeros(n,n);
for i = 1:n
  for j = 1:n
    if (Y(j,i) >= -r(i) && Y(j,i) <= r(i))
      Mach_nc(j,i) = mach_nc(i);
      Mach_c(j,i) = mach_c(i);
    else
      Mach_nc(j,i) = 0;
      Mach_c(j,i) = 0;
    end
  endfor
endfor

figure(10);
subplot(2,1,1);
colormap(viridis);
contourf(Mach_nc);
title('Non conservative form mach number in nozzle');
subplot(2,1,2);
colormap(viridis);
contourf(Mach_c);
title('Conservative form mach number in nozzle');

For above plot the error for property is calculated by taking difference between analytical value and numerical value. And then that error is plotted on Y-axis against non dimensional nozzle length.

Above plot shows that non dimensional mass flow rate for conservative form shows lots of irregularities than non conservation scheme.

In conservation part we convert our primitive terms into flux terms and then solve it by numerical discretization and again convert those flux values into primitive variables. But in non conservation form we directly solve primitive variables by numerical discretization, hence due to conversion of flux term, source term from primitive variable and after final answer again converting to primitive variables conservation part requires more computational memory and time as compared to the non conservative part for calculations. The above time comparison shows the same i.e. Non conservation form required less time than conservation form.

There are many parameters affect the exact solution and creates error in numerical schemes. Such parameters are listed below-

1. Inflow boundary condition error:- At x=0 we have assumed that density, pressure and temperature are equal to reservoir properties, but its strictly possible if Mach number = 0. In above problem as the mach number is function of area ratio so mach number has finite value.
2. Courant number:- For linear hyperbolic equation CFL should be `<=0`. But as seen in given problem as C increases the error is also increasing. Hence, we have taken C=0.5 for this problem.
3. Truncation error/ Grid independence:- Doubling the number of grid points from 31 to 61 changes solution marginally. By comparing the below table, it can be seen that non conservation part provides more close solution to analytic solution than conservation part.

	Flow parameters at nozzle throat
	`(rho^star)/(rho_0)`	`(T^star)/(T_0)`	`(P^star)/(P_0)`	`M`
Non conservation part n=31	0.63866	0.83645	0.53421	0.99939
Conservation part n=31	0.69286	0.85881	0.59504	0.91377
Non conservation part n=61	0.63766	0.83544	0.53272	0.99977
Conservation part n=61	0.63998	0.83535	0.53461	0.99868
Exact analytical solution	0.63394	0.83333	0.52828	1

Conclusion:-

From referring to all above plots and table, we can conclude that for above Quasi 1D nozzle flow problem Non conservation scheme is more accurate, and faster than conservation scheme.

Reference:- A book of 'Computational fluid dynamics' by "John D. Anderson".