Quasi 1-D Supersonic nozzle flow simulation

To convert subsonic flow into supersonic flow, we use C-D nozzle. Here, we simulate 1-D supersonic flow using two approches, using conservative and non-conservative forms of flow equations.

We use basic principles of Continuity, Momentum and Energy equations that are in 3-D adapted for this One Dimensional simulation. As named, in 1-D we only consider the change of properties such as Pressure, Temperature and Mass flow rate as a function of length.

Above figure shows how a C-D nozzle looks like. C-D nozzles mostly find their applications in rockets, supersonic aircrafts, missiles, military aircraft and supersonic wind tunnels. If we want to achieve the supersonic flow through the nozzle, the throat have to be choked. Throat is a part of C-D nozzle which is having minimum area of the nozzle.

In this project we take the nozzle in which area varies with respect to this equation `A=1+2.2*(x-1.5)^2` where x is the length of nozzle varies from 0 to 3.

The equations of Continuity, Momentum and Energy that are adapted for this 1-D flow are

1.Continuity

2.Momentum

3. Energy

• The prime on the variables indicate the non-dimensionalised form. The above shown equations are in non-conservative form.We perform this simulation by using MacCormack\'s explicit technique. In this method we have two step Predictor step and Corrector step
• In predictor step we use forward difference for space terms and after then we update the obtained solution. After updating the step we proceed with corrector step.
• In corrector step we use backward difference for space terms and after then we update the obtained solution. In this corrector step we used the updated values of previous predictor step.
• After obtaining the corrector step we take the average of both corrector step and predictor step. After obtaining the average value one full time step will be completed. For next time step the above proces will be repeated again till we get the convergence.

FDE\'s:- Difference equations for above non-conservative form are given below.

In the same way the remaining energy and momentum equations are discretized and will follow predictor and corrector step as above.

Boundary Conditions:-

The boundary conditions of any PDE\'s are dependent upon the type of the equation. For Hyperbolic PDE\'s there\'s something called method of characteristics. These are the lines which are responsible. The one dimensional inviscid unsteady flows comes into the category of hyperbolic equations. For hyperbolic equations characteristic llines will be shown in the below figure.For one dimensional inviscid unsteady flow we always march in the direction of time whether it\'s a subsonic or supersonic flow. For hyperbolic pde\'s there exist two characteristic lines at any point in xt plane.

Here these characteristic lines are physically called as Mach waves. Method of characteristics tells us at a boundary if one characteristic is propagates into the boundary then one dependent variable must be specified at that boundary. If a characteristic is propagates out of the flow then another variable must be allowed to float i.e., it must be calculated in steps of time. At subsonic entry i.e., inlet of the nozzle we need to give one flow variable and another must be floated but at the exit the flow is supesonic, because of supersonic veloities of the flow the the characteristic is carried downstream. At the exit of the nozzle both of the characteristics are propagating out, so we must allow two flow field variables to float.

Conservative form:-

In this we will be doing the same simulation using conservative form instead of using non-conservative form of governing equations.

Above given equations are conservational form of Governing equations i.e., Continuity, Momentum and Energy.

For these equations we can\'t discretize them easily, we first convert these equations into other form by defining above terms as flux and source terms.

We can write the governing equations in quasi 1-D unsteady as below

Before directly discretizing the above equation we must express the primitive variables(rho,p,T) in terms of solution vectors(U1,U2 etc). Before we proceed we need stop and think about our flux vectors. If we don\'t express the flux vectors(F1,F2,F3) in terms of our solution vectors otherwise solution will blow up at subsonic section.

Main code:-

clear all
close all
clc
% inputs
n=481;
c=0.5;
x=linspace(0,3,n);
dx=x(2)-x(1);
gamma=1.4;
% area
a=1+2.2*(x-1.5).^2;
nt=8000;
th=(n-1)/2;

prompt=\'which form to be used? 1.Non-Conservative 2. Conservative\';
solver=input(prompt);
if solver==1;
  non_conservative_1d (n,c,x,dx,gamma,a,nt,th)
else if solver==2;
      conservative_1d (x,dx,n,gamma,a,nt,th,c)
 end
end

Non-Conservative code:-

function out=non_conservative_1d (n,c,x,dx,gamma,a,nt,th)

rho=1-0.3146*x;
t=1-0.2314*x;
v=(0.1+1.09*x).*t.^0.5;

for z=1:n
  dt(z)=(c*dx)/(t(z)^0.5+v(z));
end
dt=min(dt);
  
  
% outer time loop
for k=1:nt
  rho_old=rho;
  t_old=t;
  v_old=v;
  
 
  % predictor step
   
   
   for e=2:n-1
     dvdx=(v(e+1)-v(e))/dx;
     drhodx=(rho(e+1)-rho(e))/dx;
     dtdx=(t(e+1)-t(e))/dx;
     dadx=(log(a(e+1))-log(a(e)))/dx;
     
     % continuity equation
           
      drhodt_p(e)=-rho(e)*dvdx-rho(e)*v(e)*dadx-v(e)*drhodx;
      
      % momentum equation
      
      dvdt_p(e)=-v(e)*dvdx-(1/gamma)*(dtdx+((t(e)/rho(e))*drhodx));
      
      % Energy equation
      
      dtdt_p(e)=-v(e)*dtdx-(gamma-1)*t(e)*(dvdx+v(e)*dadx);
      
      % soution update
      rho(e)=rho(e)+drhodt_p(e)*dt;
      v(e)=v(e)+dvdt_p(e)*dt;
      t(e)=t(e)+dtdt_p(e)*dt;     
      
   
     end  
% corrector step

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;
     dadx=(log(a(j))-log(a(j-1)))/dx;
     
     % continuity equation
           
      drhodt_c(j)=-rho(j)*dvdx-rho(j)*v(j)*dadx-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)*dadx);
      
      
 end
% average time derivatives

      drhodt=0.5*(drhodt_p+drhodt_c);
      dvdt=0.5*(dvdt_p+dvdt_c);
      dtdt=0.5*(dtdt_p+dtdt_c);

      
      %final update
     for l=2:n-1
      rho(l)=rho_old(l)+drhodt(l)*dt;  
      v(l)=v_old(l)+dvdt(l)*dt;
      t(l)=t_old(l)+dtdt(l)*dt;
     end
     % Boundary conditions
 
   % inlet
  v(1)=2*v(2)-v(3);
  
  % outlet
  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);
  
  
  
    p=rho_old.*t_old;
    Mach_no= v_old./(t_old).^0.5;
    mass_flow_rate=rho_old.*a.*v_old;
  
  
 % calculating values at throat
     p_th(k)=p(th);
     Mach_no_th(k)=Mach_no(th);
     rho_th(k)=rho(th);
     mass_flow_rate_th(k)=mass_flow_rate(th);     
  
end
figure(3)
subplot(4,1,1)
plot(p_th,\'color\',\'r\')
xlabel(\'no. of time steps\')
ylabel(\'p/P_o\')

subplot(4,1,2)
plot(Mach_no_th,\'color\',\'b\')
xlabel(\'no. of time steps\')
ylabel(\'Mach Number\')

subplot(4,1,3)
plot(rho_th,\'color\',\'y\')
xlabel(\'no. of time steps\')
ylabel(\'Non dimensionalised density\')

subplot(4,1,4)
plot(mass_flow_rate_th,\'color\',\'c\')
xlabel(\'no. of time steps\')
ylabel(\'Non Dimensionalised Mass flow rate\');





end

Conservative form

function out=conservative_1d(x,dx,n,gamma,a,nt,th,c)
  
 % calculating initial profiles
for i=1:n
    % inputs
    if (x(i)>=0 && x(i)<=0.5)
        rho(i)=1;
        t(i)=1;
    end
    if (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);
    end
    
    if (x(i)>=1.5 && x(i)<=3.5)
        rho(i)=0.634-0.3879*(x(i)-1.5);
        t(i)=0.833-0.3507*(x(i)-1.5);
    end
end
 
v=0.59./(rho.*a);
p=rho.*t;
c=0.5;
for z=1:n
    dt(z)=(c*dx)/(t(z)^0.5+v(z));
end
dt=min(dt);

    U1= rho.*a;
    U2=rho.*a.*v;
    U3=rho.*((t./(gamma-1))+(0.5*gamma.*v.^2)).*a;
    
    
   
% outer time loop

for k=1:nt
   
    rho_old=rho;
    v_old=v;
    t_old=t;
    p_old=p;
    U1_old=U1;
    U2_old=U2;
    U3_old=U3;
   
    F1=U2;
    F2=((U2.^2)./U1)+((gamma-1)./gamma).*(U3-(0.5*gamma.*U2.^2.*(U1.^-1)));
    F3=(gamma.*U2.*U3.*(U1.^-1))-(0.5*gamma*(gamma-1).*U2.^3.*(U1.^-2));
    
    
    for i=2:n-1
        
        J2(i)=(1/gamma).*rho(i).*t(i).*((a(i+1)-a(i))/dx);
        
        %predictor step
        
        dU1_dt_p(i)=-(F1(i+1)-F1(i))/dx;
        
        dU2_dt_p(i)=((-(F2(i+1)-F2(i))/(dx))+J2(i));
        
        dU3_dt_p(i)=-((F3(i+1)-F3(i))/dx);
        
        % updating values
        
        U1(i)=U1(i)+(dU1_dt_p(i))*dt;
        
        U2(i)=U2(i)+(dU2_dt_p(i))*dt;
        
        U3(i)=U3(i)+(dU3_dt_p(i))*dt;    
    rho(i)=U1(i)./a(i);
    v(i)=U2(i)./U1(i);
    t(i)=(gamma-1)*((U3(i)./U1(i))-((gamma/2)*(v(i).^2)));
    
    
    
    F1(i)=U2(i);
    F2(i) = (U2(i).^2./U1(i)) + ((gamma-1)/gamma)*(U3(i) - (gamma*U2(i).^2)./(2*U1(i)));
    F3(i) = (gamma*U2(i).*U3(i)./U1(i)) - (gamma*(gamma-1)/2)*(U2(i).^3./U1(i).^2);
    
    
    end
  
    
    
 for i=2:n-1        
   
        J2(i)=(1/gamma).*rho(i).*t(i).*((a(i)-a(i-1))/dx);
        
        %corrector step
        
        dU1_dt_c(i)=-(F1(i)-F1(i-1))/dx;
        
        dU2_dt_c(i)=((-(F2(i)-F2(i-1))/(dx))+J2(i));
        
        dU3_dt_c(i)=-((F3(i)-F3(i-1))/dx);
        
        
        
 end
   
 
 % average values
 
        dU1_dt=0.5*(dU1_dt_p+dU1_dt_c);
        
        dU2_dt=0.5*(dU2_dt_p+dU2_dt_c);
        
        dU3_dt=0.5*(dU3_dt_p+dU3_dt_c);   
        
        
        for i=2:n-1
            U1(i)=dU1_dt(i).*dt+U1_old(i);
            U2(i)=dU2_dt(i).*dt+U2_old(i);
            U3(i)=dU3_dt(i).*dt+U3_old(i);
        end
   % applying boundary conditions
   
   U2(1)= 2*U2(2)-U2(3);
   
   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);
        
    % updation of values 
    rho=U1_old./a;
    v=U2_old./U1_old;
    t=(gamma-1)*((U3_old./U1_old)-(gamma/2)*v_old.^2);
    p=rho_old.*t_old;
    Mach_no= v_old./sqrt(t);
    mass_flow_rate=rho_old.*a.*v_old;
    
     % calculating values at throat
     p_th(k)=p(th);
     Mach_no_th(k)=Mach_no(th);
     rho_th(k)=rho(th);
     mass_flow_rate_th(k)=mass_flow_rate(th);  
    
  end
   
  
    figure(1)
    subplot(4,1,1)

    plot(x,rho,\'b\',\'linewidth\',2)
     grid on
 subplot(4,1,2)
    plot(x,v,\'r\',\'linewidth\',2)
      grid on
 subplot(4,1,3)
   plot(x,rho,\'g\',\'linewidth\',2)
     grid on
 subplot(4,1,4)
 plot(x,t,\'y\',\'linewidth\',2)
   grid on
   
  
  
   figure(2)
   grid on
     subplot(4,1,1)
     plot(Mach_no_th,\'color\',\'c\')
     xlabel(\'Mach Number\')
     ylabel(\'no. of time steps\')
      grid on
     subplot(4,1,2)
    plot(p_th,\'color\',\'r\')
    xlabel(\'non-dimensional pressure\')
    ylabel(\'no. of time steps\')
    
    
    grid on
    
    subplot(4,1,3)
    plot(rho_th,\'color\',\'b\')
    grid on
    
    subplot(4,1,4)
    plot(mass_flow_rate_th,\'color\',\'y\')
    grid on
    
    
    pause(0.0003)
  
         
end 

   
   

CONSERVATIVE-1 n=31

Non-Conservative-1 n=31

CONSERVATIVE-2 n=31

Non-Conservative-2 n=31

Conservative-3 n=31

Non-Conservative-3 n=31

Conservative-4 n=61

Non-Conservative-4 n=61

Conservative-5 n=121

Non-Conservative-5 n=121

• When we take no of grid points as 31 the simulation is taking nearly 600 time steps to achieve steady state and converge.
• In Non-Conservative code the simulation is taking approximately 100 timesteps less than the Conservative approach.
• If we take no. of grid points as 61, simulation is taking nearly 800 time steps in non-conservative form and 600 steps in conservative form.
• We achieved grid independence at 121 grid points. Grid independence is can be seen in above plots. It is the effect of unchanged results even after you increase the no of grid points.
• For n=31 with 1400 time steps, the Mach,Density and Pressure are not approximately equal to the analytical values.
• When n=61 with 1400 time steps, there\'s change in the results. This can be observed by looking the plot values closely.
• When n=121 with 2800 time steps, there\'s change in the results and are close to the analytical values such as Mach at throat is 1 but here in results we achieved 0.96. Which almost close.
• With these points taking into account we can conclude that grid independence is achieved at n= 121.
• Mass flow is coverging early in the non-conservative form than conservative this is because of the more complex code and more no. of variables to calculate in the conservative form.
• Overall non-conservative form is taking less time to converge than the conservative form of governing equations.