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

OBJECTIVE: This project's main aim is to perform numerical-solution analysis of a Quasi 1D Supersonic Nozzle Flow using the Maccormack Method. Although analytically identical, the numerical solution of conservation and non-conservation forms of the governing equations is usually different. That is the reason why both ways are solved and compared in this project. For that, the number of iterations and computation time will be examined, as long as the final profile for the normalized mass flow rate. During this project, both a grid dependence test and a CFL-based adaptative time loop control was programmed to assure the fidelity of the numerical results. Besides, this report also includes the variation of essential fluid properties at the throat section of our nozzle and mass flow rates at different time steps.

PROBLEM STATEMENT:

Consider a steady isentropic flow through a convergent-divergent nozzle, as shown in Figure 1. The flow at the inlet to the nozzle comes from a reservoir where the pressure and temperature are constants. The cross-sectional area of the reservoir is so large that the velocity at the inlet is almost Zero (Vin​=∼0). The flow expands isentropically to supersonic speeds at the nozzle exit. The flow is locally subsonic in the convergent section of the nozzle (M < 1).sonic at the throat section (M = 1) and supersonic at the divergent section (M > 1).Though in reality, the flow field varies in the two-dimensional (XY plane), an assumption is made that the flow properties varies only along the x-axis, this flow is defined as quasi-one-dimensional flow.

THEORY

Quasi 1D Supersonic Nozzle flow

Quasi one dimensional implies that we still have variations of flow quantities in one direction only. Still, we allow the cross-section area of stream tubes to vary along the same path as well. This means that we will be able to analyze axisymmetric channels with changing cross-section areas such as example, convergent-divergent nozzles. The convergent-divergent nozzle or CD nozzle is a tube that is pinched in the middle, making a carefully balanced, asymmetric hourglass shape. It is used to accelerate a hot, pressurized gas passing through it to a higher supersonic speed in the axial (thrust) direction, by converting the heat energy of the flow into kinetic energy. Because of this, the nozzle is widely used in some types of steam turbines and rocket engine nozzles. It also sees use in supersonic jet engines.

Since, area(A), Mach number(M ) and velocity(v) changes along the x-axis and thus there exist a relation which relates all the three variables, written as,

`{dA}/{A}=(M^2 −1)*(dv)/v`

This relation is an important contribution to the compressible flow engineering field since it describes gives us information on how to accelerate a flow to supersonic speeds. This can be further illustrated by observing that the nozzle first has a convergent part (where A increases along positive x ) in which we observed sub-sonic flow (M<1) thus from our previous relation, it can be concluded that v increases along positive x while for the second part of the nozzle having a divergent part (where A decreases along positive x ), we have supersonic flow (M>1) and thus from above relation, v will further increase along the positive x-axis.

In our problem statement, we consider the isentropic flow(whenever there is a reversible and adiabatic flow, the constant value of entropy is maintained) flowing throw CD nozzle. This is solved by numerical solution techniques using different forms of governing equations to reach a steady-state solution after some time interval.

Non-Conservative form of governing equation

continuity equation

`(drho')/(dt')=-rho'*(dV')/(dx')-rho'*V'*(d(logA'))/(dx')-V'*(drho')/(dx')`

Momentum equation

`(dV')/(dt')=-V'(dV')/(dx')-(1/gamma)((dT')/(dx')+(T')/(rho')(drho')/(dx'))`

energy equation

`(dT')/(dt')=-V'(dT')/(dx')-(gamma-1)T'((dV')/(dx')+V'(d(logA'))/(dx'))`

All the variables are in non dimensional form.

Initial conditions

`rho=1-0.3146x`

`t=1-0.2314x`

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

boundary conditions

Inlet
v(1)=2*v(2)-v(3)
Outlet
v(n)=2*v(n-1)-v(n-2)
rho(n)=2*rho(n-1)-rho(n-2)
t(n)=2*t(n-1)-t(n-2)

Conservative form of governing equation

continuity equation

`(d(rho'A'))/(dt')=-(d(rho'A'V'))/(dx')`

Momentum equation

`(d(rho'A'V'))/(dt')+(d(rho'A'V'^2+(1/gamma)rho'A'))/(dx')=(1/gamma)rho'(dA')/(dx')`

energy equation

`(d(rho'(T'1/(gamma-1)+(gamma/2)V'^2)A'))/(dt')=-(d(rho'(T'1/(gamma-1)+(gamma/2)V'^2)V'A'+rho'A'V'))/(dx')`

The above equation can be simpified into the following equations.

`U_1=rho'A'`

`U_2=rho'A'V'`

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

`F_1=rho'A'V'`

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

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

`J_2=1/gammarho'(dA')/(dx')`

where `U_1`, `U_2`, `U_3`is called the solution vector, `F_1`, `F_2` ,`F_3` are flux vectors and `J_2` is called the source term.

substituting the above term in governing equation we get;

`(dU_1)/(dt')=-(dF_1)/(dx')`

`(dU_2)/(dt')=J_2-(dF_2)/(dx')`

`(dU_3)/(dt')=-(dF_3)/(dx')`

Initial conditions

`rho=1-0.3146x`

`t=1-0.2314x`

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

boundary conditions

Inlet
v(1)=2*v(2)-v(3)
Outlet
v(n)=2*v(n-1)-v(n-2)
rho(n)=2*rho(n-1)-rho(n-2)
t(n)=2*t(n-1)-t(n-2)

Discretization using Macormack Method

Steps in Macormack Method

Predictor step

`(df)/(dt_p)`

Solution Update

`f^'=f_(old)+((df)/(dt_p))*Deltat`

Corrector step

`(df^')/(dt_c)`

calculating average

`((df)/dt)_(avg)=0.5((df)/(dt_p)+(df')/(dt_c))`

final solution update

`f=f_(old)+((df)/dt)_(avg)*Deltat`

MATLAB CODE:

Main Program

close all
clear all
clc
n=31;
x=linspace(0,3,n);
tol=1e-4;
m_analytical=0.579*ones(1,n);
[m_nc]=supersonic_nozzle_non_conservative(n,tol);
[m_c]=supersonic_nozzl_conservative(n,tol);
    figure(5)
    plot(x,m_nc,'linewidth',2)
    hold on
    plot(x,m_c,'linewidth',2)
    hold on
     plot(x,m_analytical,'linewidth',2)
    title('mass flow rate vs x')
    xlabel('x')
    ylabel('Mass flow rate')
    legend('Non-conservative','conservative','analytical')
    

function for Non-conservative form

function [m_nc]=supersonic_nozzle_non_conservative(n,tol)
tic
%grid
x=linspace(0,3,n);

dx=x(2)-x(1);
gamma=1.4;
%specify throat
throat=(n+1)/2;
 m_analytical=0.579*ones(1,n);

C=0.5;
%initial profiles
rho=1-0.3146*x;
t=1-0.2314*x;   %Temperature
v=(0.1+1.09*x).*t.^0.5;
rho_old=rho;
v_old=v;
t_old=t;
%area
a=1+2.2*(x-1.5).^2;
%time steps
nt=1400;
nt_step=linspace(1,nt,nt);
for p=1:n
    dt_a(p)=C*dx/(t(p)^0.5+v(p));
end
dt=min(dt_a);
%outer time loop
iter=1;
error_rho=9e9;
error_v=9e9;
error_t=9e9;
for k=1:nt
     rho_old=rho;
     v_old=v;
     t_old=t;
   
  
    %predictor method
    for j=2:n-1
        dvdx=(v(j+1)-v(j))/dx;
        drhodx=(rho(j+1)-rho(j))/dx;
        dtdx=(t(j+1)-t(j))/dx;
        dlogadx=(log(a(j+1))-log(a(j)))/dx;
        
        %continuity equation
        drhodt_p(j)=-rho(j)*dvdx-rho(j)*v(j)*dlogadx-v(j)*drhodx;
        
        %momentum equation
        dvdt_p(j)=-v(j)*dvdx-(1/gamma)*(dtdx+(t(j)/rho(j))*drhodx);
        
        %energy equation
        dtdt_p(j)=-v(j)*dtdx-(gamma-1)*t(j)*(dvdx+v(j)*dlogadx);
        
        %solution update
        rho(j)=rho(j)+drhodt_p(j)*dt;
        v(j)=v(j)+dvdt_p(j)*dt;
        t(j)=t(j)+dtdt_p(j)*dt;
    end
    %corrector method
    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 equation
        dtdt_c(j)=-v(j)*dtdx-(gamma-1)*t(j)*(dvdx+v(j)*dlogadx);
        
    end
 %compute avg time derivative
 drhodt=0.5*(drhodt_p+drhodt_c);
 dvdt=0.5*(dvdt_p+dvdt_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
 M=v./t.^0.5;
 m_nc=rho.*a.*v;
  m_nct(k,:)=rho.*a.*v;
   error_rho=max(abs(rho-rho_old));
  error_v=max(abs(v-v_old));
  error_t=max(abs(t-t_old));
  if (error_rho>tol && error_v>tol && error_t>tol)
      iter=iter+1;

  end
 %apply boundary conditions
 %inlet
 v(1)=2*v(2)-v(3);
 %outlet
  v(n)=2*v(n-1)-v(n-2);
   rho(n)=2*rho(n-1)-rho(n-2);
    t(n)=2*t(n-1)-t(n-2);
    
   
end
rho_throat=rho(throat);
  v_throat=v(throat);
   t_throat=t(throat);
    M_throat=M(throat);
 figure(1)
    subplot(4,1,2)
    plot(x,rho,'linewidth',2)
    
    xlabel('nozzle length')
    ylabel('density ratio')
    subplot(4,1,3)
    plot(x,v,'r','linewidth',2)
   
    xlabel('nozzle length')
    ylabel('velocity ratio')
    subplot(4,1,4)
    plot(x,t,'g','linewidth',2)
   
    xlabel('nozzle length')
    ylabel('Temperature ratio')
    subplot(4,1,1)
    plot(x,M,'m','linewidth',2)
    title(sprintf('No of timestep=%d, no of iteration to steady state=%d',k,iter))
    xlabel('nozzle length')
    ylabel('Mach number')
    suptitle('Non-conservative form')
    
  
     figure(3)
     plot(x,m_nct(50,:),x,m_nct(100,:),x,m_nct(200,:),x,m_nct(300,:),x,m_nct(500,:),'o',x,m_analytical,'linewidth',2)
     xlabel('non-dimensional distance x')
     ylabel('non-dimentional mass flow rate')
     title('variation of mass flow rate(non-conservative)')
     legend('50dt','100dt','200dt','300dt','700dt','analytical')
    time_nc=toc;
    fprintf('time of simlulation for nonconservative=%gsn',time_nc);
end

function for conservative form

function [m_c]=supersonic_nozzl_conservative(n,tol)
tic
%grid 
x=linspace(0,3,n);
dx=x(2)-x(1);
gamma=1.4;
throat=(n+1)/2;
m_analytical=0.579*ones(1,n);
%initial profiles
for q=1:n
    if ((x(q)>=0) && (x(q)<=0.5))
    rho(q)=1;
    t(q)=1;
    elseif ((x(q)<1.5) && (x(q)>0.5))
        rho(q)=1-0.366*(x(q)-0.5);
        t(q)=1-0.167*(x(q)-0.5);
    elseif ((x(q)>=1.5) && (x(q)<=3.5))
        rho(q)=0.634-0.3879*(x(q)-1.5);
        t(q)=0.833-0.03507*(x(q)-1.5);
    end
end
%area
a=1+2.2*(x-1.5).^2;
 v=0.59./(rho.*a);
U1=rho.*a;
U2=rho.*a.*v;
U3=rho.*(t./(gamma-1)+(gamma/2)*v.^2).*a;
%time steps
nt=1400;
nt_step=linspace(1,nt,nt);
C=0.5;
for p=1:n
    dt_a(p)=C*dx/(t(p)^0.5+v(p));
end
    dt=min(dt_a);
%outer time loop
iter=1;
 error_rho=9e9;
 error_v=9e9;
 error_t=9e9;
for k=1:nt
     U1_old=U1;
     U2_old=U2;
     U3_old=U3;
     rho_old=U1_old./a;
    v_old=U2_old./U1_old;
    t_old=(gamma-1)*((U3_old./U1_old)-(gamma/2)*v_old.^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);
    
  
    %predictor method
    for j=2:n-1
       J2(j)=(1/gamma)*rho(j)*t(j)*(a(j+1)-a(j))/dx;
                
        %continuity equation
        dU1dt_p(j)=-(F1(j+1)-F1(j))/dx;
        
        %momentum equation
        dU2dt_p(j)=J2(j)-(F2(j+1)-F2(j))/dx;
        
        %energy equation
        dU3dt_p(j)=-(F3(j+1)-F3(j))/dx;
        
        
        
        %solution update
        U1(j)=U1(j)+dU1dt_p(j)*dt;
        U2(j)=U2(j)+dU2dt_p(j)*dt;
        U3(j)=U3(j)+dU3dt_p(j)*dt;
    end
    rho=U1./a;
    v=U2./U1;
    t=(gamma-1)*((U3./U1)-(gamma/2)*v.^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 method
    for j=2:n-1
         J2(j)=(1/gamma)*rho(j)*t(j)*(a(j)-a(j-1))/dx;
         %continuity equation
        dU1dt_c(j)=-(F1(j)-F1(j-1))/dx;
        
        %momentum equation
        dU2dt_c(j)=J2(j)-(F2(j)-F2(j-1))/dx;
        
        %energy equation
        dU3dt_c(j)=-(F3(j)-F3(j-1))/dx;
       
        
    end
 %compute avg time derivative
 dU1dt=0.5*(dU1dt_p+dU1dt_c);
 dU2dt=0.5*(dU2dt_p+dU2dt_c);
 dU3dt=0.5*(dU3dt_p+dU3dt_c);
 
 %final solution update
 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

 %apply boundary conditions
 %inlet

U1(1)=rho(1).*a(1);
  U2(1)=2*U2(2)-U2(3);
   v(1)=U2(1)/(rho(1)*a(1));
U3(1)=U1(1)*(t(1)/(gamma-1)+(gamma/2)*v(1)^2);
 %outlet
  U2(n)=2*U2(n-1)-U2(n-2);
   U1(n)=2*U1(n-1)-U1(n-2);
    U3(n)=2*U3(n-1)-U3(n-2);
    %updating vlues
    rho=U1./a;
    v=U2./U1;
    t=(gamma-1)*((U3./U1)-(gamma/2)*v.^2);
    M=v./t.^0.5;
 m_c=rho.*a.*v;
 m_nc(k,:)=rho.*a.*v;
 
    error_rho=max(abs(rho-rho_old));
    error_v=max(abs(v-v_old));
 error_t=max(abs(t-t_old));
  if (error_rho>tol && error_v>tol && error_t>tol)
      iter=iter+1;
  end
   
end
rho_throat=rho(throat);
  v_throat=v(throat);
   t_throat=t(throat);
    M_throat=M(throat);
 figure(2)
   subplot(4,1,2)
    plot(x,rho,'linewidth',2)
    
    xlabel('nozzle length')
    ylabel('density ratio')
    subplot(4,1,3)
    plot(x,v,'r','linewidth',2)
   
    xlabel('nozzle length')
    ylabel('velocity ratio')
    subplot(4,1,4)
    plot(x,t,'g','linewidth',2)
  
    xlabel('nozzle length')
    ylabel('Temperature ratio')
    subplot(4,1,1)
    plot(x,M,'m','linewidth',2)
    title(sprintf('No of timestep=%d ,no of iteration to steady state=%d',k,iter))
    xlabel('nozzle length')
    ylabel('Mach number')
    suptitle('Conservative form')
    
    
     figure(4)
     plot(x,m_nc(50,:),x,m_nc(100,:),x,m_nc(200,:),x,m_nc(300,:),x,m_nc(700,:),'o',x,m_analytical,'linewidth',2)
     xlabel('non-dimensional distance x')
     ylabel('non-dimentional mass flow rate')
     title('variation of mass flow rate(conservative)')
     legend('50dt','100dt','200dt','300dt','700dt','analytical')
      time_c=toc; 
      fprintf('time of simlulation for conservative=%gsn',time_c);
end

Results:

The profile of flow variables for Non-conservative form and Conservative form.

Convergence-

Through the above two graphs, the convergence time step can be analyzed. We can see that at a 700-time step the mass flow rate corresponds to the analytical solution which means a steady-state is achieved.

The convergence for non-conservative form and conservative form was found to be 404 and 471 iterations when the convergence criteria were used 1e-4. So we can conclude that the non-conservative form has a faster convergence than the conservative form.

Comparison of mass flow rate between non-conservative form and conservative form.

From the above graph, it can be concluded that the mass flow rate in conservative form showed better results and corresponded to the analytical solution.

Simulation time for conservative and non-conservative form

Grid Dependence test

Non-conservative form(31 grids)

Conservative form (31 grids)

Non-conservative form (61 grids)

Conservative form(61 grids)

From the grid dependence test, it can be seen that for 31 and 61 grid points respectively there was hardly any change in density, velocity, temperature, and Mach number at the throat. This is noticed for both non-conservative as well as conservative. So the solution has no change for any number of grids.