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

We consider the steady isoentropic flow through the convergent -divergent nozzle as given in diagram below.

   

The flow at the inlet to the nozzle comes from the reservoir is large (therotically $A\to \infty$0 and hence the velocity is very small $(V\to 0)$ .Thus ,$p_o$ and $T_o$ are the stagnantation values or total pressure and total temperature,respectively.The flow expands isoentropically to supersonic speeds at the nozzle exit,where the exit pressure,temperature where exit parameters are reperesented with subscript e.

The flow is locally subsonic in the convergent section of the nozzle,sonic at the throat (minimum area) ,and supersonic at the divergent section . The sonic flow (M=1)at  the throat means the local velocity at this section is equal to the local speed of the sound.We can repersent the throat values with our own representation either with * or any other.

 

Some of the useful results obtained.

 

${\left(\frac{A}{A^{\ast }}\right)}^2=\frac{1}{M^2}\cdot {(\frac{2}{\gamma +1}\cdot (1+\frac{\gamma -1}{2}\cdot M^2))}^{\frac{\gamma +1}{\gamma -1}}$

 

where $\gamma$ is the ratio of specific heat.For air at standard conditons it is equal to 1.4.

Some other equations include 

$\left(\frac{p}{p_o}\right)={(1+\frac{\gamma -1}{2}\cdot M^2)}^{-\frac{\gamma }{\gamma -1}}$

$\left(\frac{\rho }{{\rho }_o}\right)={(1+\frac{\gamma -1}{2}\cdot M^2)}^{-\frac{1}{\gamma -1}}$

$\left(\frac{T}{T_o}\right)={(1+\frac{\gamma -1}{2}\cdot M^2)}^{-1}$

 

Mach number:

Mach number  is a dimensionless quantity in fluid dynamics representing the ratio of flow velocity past a boundry to the local speed of sound

       $M=\frac{u}{c},$

where:

M  is the local Mach number,

u  is the local flow velocity with respect to the boundaries (either internal, such as an object immersed in the flow, or external, like a channel), and

c  is the speed of sound in the medium, which in air varies with the square root of the thermodynamics temeprature.

 

Assuming air to be an ideal gas, the formula to compute Mach number in a subsonic compressible flow is derived from Bernoullis equation for M < 1>

$M=\sqrt{\frac{2}{\gamma +1}\cdot [{(\frac{q_c}{p}+1)}^{\frac{\gamma -1}{\gamma }}-1]}$

 

and the  speed of sound varies with teperature as follows:

$c=\sqrt{\gamma \cdot R_{\ast }T}$

 

where$q_c$ is imact pressure

$p$ is the static pressure.

$\gamma$ is the ratio of specific heat of a gas at a constant pressure to heat at a constant volume (1.4 of air).

$R_{\ast }$ is the specific gas constant.

Set Up

MaCcormack mathod:

In computational fluid dynamics, the MacCormack method is a widely used discretization scheme for the numerical solution of hyperbolic partial differential equations.

 

Maccormacks method provides second order accuracy without the invlovement of second order derivative thus reducing the complexity  of the calculation.It is a explicit method.

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

$F^{actual}=F^{Initial}+\frac{\partial \left(f\right)}{{(\partial t)}_{avg}}\cdot \partial t$

 

The average time derivative can be found as follows ,which is the  half of sum of predicotor value and correcotr value .

${\left(\frac{\partial v}{\partial t}\right)}_{avg}={\left(\frac{\partial v}{\partial t}\right)}_p+{\left(\frac{\partial v}{\partial t}\right)}_c$

 

So basically before that we need to set up equations for solving the predictor and corrector method.

Governing equations:

The super sonic flow through nozzle is governed by Navier stokes equations,which could be derived for both the cases of conservative and Non conservative methods.

The continuity equation can be given as follows:

$\frac{\partial {\rho }^{\prime }}{\partial t^{\prime }}=-\left({\rho }^{\prime }\right)\cdot \frac{\partial V^{\prime }}{\partial x^{\prime }}-{\rho }^{\prime }\cdot V^{\prime }\cdot \frac{\partial \left({\ln A}^{\prime }\right)}{\partial x\text{'}}-V^{\prime }\cdot \left(\frac{\partial {\rho }^{\prime }}{\partial x^{\prime }}\right)$

 

The momentum equation:

$\frac{\partial V^{\prime }}{{(\partial t)}^{\prime }}=-V^{\prime }\cdot \frac{\partial V^{\prime }}{\partial x^{\prime }}-\frac{1}{\gamma }\cdot (\frac{\partial T^{\prime }}{\partial x^{\prime }}+\frac{T^{\prime }}{{\rho }^{\prime }}\cdot \frac{\partial {\rho }^{\prime }}{\partial x^{\prime }})$

 

The energy equation is as follows:

$\frac{\partial T^{\prime }}{\partial t^{\prime }}=-V^{\prime }\cdot \frac{\partial T^{\prime }}{\partial x^{\prime }}-(\gamma -1)\cdot T^{\prime }\cdot [\frac{\partial V^{\prime }}{\partial x^{\prime }}+V^{\prime }\cdot \frac{\partial \left({\ln A}^{\prime }\right)}{\partial x^{\prime }}]$

 

These equation can be derived from the 

 

 

The $T_o$ and ${\rho }_o$ arethe reservoir temperature and density.We define the non dimensional temperature and density respectively as follows:

 

 

 

Next we define the initial condition :

• The shape of nozzle ,A=A( x) is specified and held fixed ,independent of time,For this case we take the parabolic area distribution 
• $A=1+2.2\cdot {(x-1.5)}^2$     $0\le x\le 3$
• For density ,Temperature and velocity is given by
• $\rho =1-0.3146\cdot x$
• $T=1-0.2314x$
• $V=(0.1+1.09x)\cdot \sqrt{T}$
• Pressure  $P=\rho \cdot T$
• Mach number= $\frac{u}{\sqrt{T}}$
• The above mentioned are the initial conditions at t=0

Next we come to predictor step:

This is a kind of descritisation step.Where we descretize with forward difference equation.

We first set up the finite difference expressions as follows:

We set up the spatial derivatives as forward differences .Also to reduce the complexity of the notation we will drop the prime to denote the dimensionless variable.

 

 

And some of the other equations including for other variables like momentum and energy as well  .

For grid generation we  divide the length of the nozzle in number of points along the axis .The characteristics of the flow changes from subsonic to super sonic from right to left .The distance between the two grid points is dx.Smaller the mesh size the error also decreases.

In discretisation part we apply thr MaCcormach technique to solve the the conservative and non conservative forms of governing equations.

 

 

From these we get the predicted equation as follows: These are denoted by barred quantities.

 

 

In these equations ${\rho }_i^t,V_i^t$ and $T_i^t$ are the known values at time t.

 

Next step is  the corrcetor method:

Here we use the backward differencing scheme and rest remains same as that of the predictor method.

 

Continuity equation

Momentum equation

 

 

Energy equation

 

 

Calculation of Average:

Its basically the half of the sum of corrector and predictor value.

1)

2)

3)

 

Then finally the corrected values at time t+ $\Delta t$ are:

 

 

  

 

Calculation of time step:

$dt=C\cdot \frac{dx}{u+a}$

where  C is the courant number itts value is 0.5

u is the local velocity 

a is the local speed of sound

Boundry condtions at subsonic section

 

 

 

 

By using linear extrapolation 

$V_1=2\cdot V_2-V_3$

Density and Temperature at point  1 is held fixed.

 

Supersonic outflow boundry conditions:

$V_N=2\cdot V_{N-1}-V_{N-2}$

${\rho }_N=2\cdot {\rho }_{N-1}-{\rho }_{N-2}$

$T_N=2\cdot T_{N-1}-T_{N-2}$

 

 

Conservative form

The equations  in  Conservative form in dimension less form is given as follows:

The continuity Equation:

• $\frac{\partial \left({\rho }^{\prime }{A}^{\prime }\right)}{\partial t^{\prime }}+\frac{\partial \left({\rho }^{\prime }{A}^{\prime }{V}^{\prime }\right)}{\partial x^{\prime }}=0$

The momentum equation

 

The energy equation :

 

 

 

The above are the non dimensional conservation form of the coninuity ,momentum and energy equations for quasi-one -dimensional flow ,respectively.The equations for quasi one dimensional flow can be expressed in generic form.Let us define the elements of the solution vector U ,flux vector F and the source term J as follows 

 

From the above equations we can write that 

These are the equations we wish to  solve numerically by using the MacCormacks technique.

To obtain the primitve variables we must decode the elements U1,U2,U3 as follows.From the defination of U1,U2 and U3 given above we have 

 

 

 

 

Expressing the flux and source term in form of solution  vectors

 

 

 

Applying the boundry conditions:

Inlet boundry conditions:

Outlet boundry condition:The flow properties at the downstream ,supersonic outflow boundry are obtained by linear extrapolation from the two adjacent internal points.

If N denotes the grid point at the outflow boundry ,then:

 

 

The values of  F1,F2 and F3 at grid point i=N  are obtained from the values of U1,U2 and U3 at point i=N using the equations.

 

Initial conditions:

 

Matlab Code :

close all
clear all
clc

%input parameters 
n=31;
x=linspace(0,3,n);
dx=x(2)-x(1);
gamma=1.4;
CFL=0.5;
Area=1+2.2*(x-1.5).^2;
throat=((n+1)/2);
nt=1400;

%main body
tic;
[V,rho,T,m,M,P,V_t,rho_t,T_t,m_t,M_t,P_t]=non_conservative(throat,x,gamma,dx,Area,n,nt,CFL);
non_conservative=toc;

tic
[Vc,rhoc,Tc,cons_mass_flow_throat,Mc,pc,Vc_throat,rho_throat,T_throat,cons_mass_flow,Mc_throat,pc_throat]=conservative(throat,x,gamma,dx,Area,n,nt,CFL);
conservative=toc;

%plotting section 
figure(3)
plot(x,m,'linewidth',0.95);
hold on;
plot(x,cons_mass_flow,'*');
%exact analytical solution rho*A*v=0.579
plot(x,linspace(0.579,0.579,n),'--','linewidth',1.5);
xlabel('Non Dimensional distance of nozzle (x)');
ylabel('mass flow rate');
legend('non conservative form','conservative form','exact analytical solution');
title('Comparision of mass flow rate(conventional form vs non conventional');
grid on ;

figure(4)
plot(x,M,'linewidth',1.5);
hold on
plot(x,Mc,'*','linewidth',0.65,'color','g');
grid on 
xlabel('distance of nozzle(x)');
ylabel('Mach number m');
legend('Non conservative form','conservative form');
title('Mach number profile');
grid on 
figure(5)

plot(x,P,'linewidth',1.25);
hold on 
plot(x,pc,'g','linewidth',1.25);
grid on 
xlabel('distance of nozzle(x)');
ylabel('pressure');
legend('non conservative form','conservative form');
title('Pressure Profile');
grid on 

figure(6)
plot(x,rho,'linewidth',1.5)
hold on 

plot(x,rhoc,'*','linewidth',1.25);
grid on 
xlabel('distance of nozzle(x)');
ylabel('density');
legend('non conservative form','conservative form');
title('density profile');
grid on 
figure(7)

plot(x,T,'linewidth',1.5);
hold on 
plot(x,Tc,'*','linewidth',0.65,'color','g');
grid on
xlabel('distance of nozzle(x) ');
ylabel('Temperature');
legend('non conservative form','conservative form');
title('Temperature change');

figure(8)
plot(x,V,'linewidth',1.0);
hold on ;
plot(x,Vc,'*','linewidth',0.65,'color','g');
grid on 
xlabel('distance of nozzle(x)');
ylabel('velocity');
legend('non conservative form','conservative form');
title('Velocity profile');

grid on

figure(9)

plot(1:nt,rho_throat,'linewidth',1.5);
hold on 
plot(1:nt,rho_t,'c','linewidth',1.65);
grid on 
ylabel('density');
legend('Conservative form','non conservative form');
title('density variation at various time step');
grid on 

figure(10)
plot(1:nt,T_throat,'linewidth',0.75);
hold on 
plot(1:nt,T_t,'linewidth',0.95)
grid on 
ylabel('temeperature');
legend('Conservative form','non conservative form');
title('Temperature variation at various time step at throat');
grid on 


figure(11)

plot(1:nt,pc_throat,'linewidth',0.85)
hold on 
plot(1:nt,P_t,'linewidth',0.95)
grid on 
ylabel('Presssure');
legend('conservative form','non conservative form');
title('pressure variatiuon at various time step at throat');
grid on

figure(12)
plot(1:nt,Mc_throat,'linewidth',0.85);
hold on 
plot(1:nt,M_t,'linewidth',0.90)
grid on 
xlabel('Time step');
ylabel('Mach number');
legend('Conservative form','non conservative form');
title('Mach number variation at various time step at throat');

The code for Conservative function:

%function for non-conservative form 

function [V,rho,T,m,M,P,V_t,rho_t,T_t,m_t,M_t,P_t]=non_conservative(throat,x,gamma,dx,A,n,nt,C)
  
�lculations

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

dt=min(C.*(dx./(sqrt(T)+V)));

for it=1:nt
    rho_old=rho;
    V_old=V;
    T_old=T;
    
    for j=2:n-1
        
        drhodx_p=(rho(j+1)-rho(j))/dx;
        dvdx_p=(V(j+1)-V(j))/dx;
        dtdx_p=(T(j+1)-T(j))/dx;
        dlogA_dx_p=(log(A(j+1))-log(A(j)))/dx;
        
        %continuity equation
        drhodt_p(j)=-rho(j)*dvdx_p-rho(j)*V(j)*dlogA_dx_p-V(j)*drhodx_p;
        %momentum equation
        dvdt_p(j)=-V(j)*dvdx_p-(1/gamma)*(dtdx_p+(T(j)/rho(j))*drhodx_p);
        %energy equation
        dtdt_p(j)=-V(j)*dtdx_p-(gamma-1)*T(j)*(dvdx_p+V(j)*dlogA_dx_p);
        %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 step.
    for j=2:n-1
        drhodx_c=(rho(j)-rho(j-1))/dx;
        dvdx_c=(V(j)-V(j-1))/dx;
        dtdx_c=(T(j)-T(j-1))/dx;
        dlogA_dx_c=(log(A(j))-log(A(j-1)))/dx;
        
        %continuity equation
        drhodt_c(j)=-(rho(j)*dvdx_c)-rho(j)*V(j)*dlogA_dx_c-V(j)*drhodx_c;
        %momentum equation       
        dvdt_c(j)=-(V(j)*dvdx_c)-(1/gamma)*(dtdx_c+(T(j)/rho(j))*drhodx_c);
        %energy equation
        dtdt_c(j)=-(V(j)*dtdx_c)-(gamma-1)*T(j)*(dvdx_c+V(j)*dlogA_dx_c);
    end
        
    %average values
       dvdt=0.5*(dvdt_c+dvdt_p);
       drhodt=0.5*(drhodt_p+drhodt_c);
       dtdt=0.5*(dtdt_c + dtdt_p);
       
       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
       V(1)=2*V(2)-V(3);
       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);
       
       m=rho.*A.*V;
       P=rho.*T;
       M=V./sqrt(T);
       
       
       rho_t(it)=rho(throat);
       V_t(it)=V(throat);
       T_t(it)=T(throat);
       P_t(it)=P(throat);
       M_t(it)=M(throat);
       m_t(it)=m(throat);
       
       rho_exact=0.634*ones(1,1400);
       t_exact=0.833*ones(1,1400);
       p_exact = 0.528*ones(1,1400);
       mach_exact=1*ones(1,1400);
      
       figure(1) 
       if it==50
           plot(x,m,'b','linewidth',1.5);
           hold on
        elseif it==100
              plot(x,m,'r','linewidth',1.5);
              hold on 
        elseif it==150
               plot(x,m,'c','linewidth',1.5);
               hold on 
        elseif  it==300
                plot(x,m,'g','linewidth',1.5);
                hold on 
        elseif  it==600
                plot(x,m,'y','linewidth',1.5);
                hold on 
        elseif it==900
                plot(x,m,'k','linewidth',1.5);
                hold on 
        elseif it==1400
                plot(x,m,'m','linewidth',1.5);
                grid on
       end
end

  xlabel('Non dimensional distance through nozzle (x)','fontsize',10);
  ylabel('Non dimensional mass flow rate through nozzle x','fontsize',10);
  title('Variation of mass flow rate (non conseravtive form) along x at different steps');
  legend({'50 time steps','100 time steps','150 time steps','300 time steps','600 time steps','900 time steps','1400 time steps'});
     
       
end

%conservative form
function [Vc,rhoc,Tc,cons_mass_flow_throat,Mc,pc,Vc_throat,rho_throat,T_throat,cons_mass_flow,Mc_throat,pc_throat]=conservative(throat,x,gamma,dx,A,n,nt,C)

R=8.314;
for i=1:n
    if (x(i)<=0.5)
        rho_in (i)=1;
        T_in(i)=1;
    elseif (x(i)>=0.5) && (x(i) <= 1.5)
            rho_in(i)=1-(0.366*(x(i)-0.5));
            T_in(i)=1-0.167*(x(i)-0.5);
        else 
                rho_in(i)=0.634-0.3879*(x(i)-1.5);
                T_in(i)=0.833-0.3507*(x(i)-1.5);
    end
end

 u_in=0.59./(rho_in.*A);
 rhoc = rho_in;
 Tc=T_in;
 Vc=u_in;
        
    



dt=0.0206;
t_sim=(nt-1)*dt;
t=0:dt:t_sim;

U1=rhoc.*A;
U2=rhoc.*A.*Vc;
U3=rhoc.*A.*((Tc/(gamma-1))+(0.5*gamma.*Vc.*Vc));

%Macormack method
for j=1:length(t)
    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 method
    for k=2:n-1
        dA_dx(j)=(A(k+1)-A(k))/dx;
        dlogA_dx = (log(A(k+1))-log(A(k)))/dx;
        J2(k)=dlogA_dx*((gamma-1)/gamma)*(U3(k)-(0.5*gamma*U2(k)^2)/(U1(k)));
        
        
        %continuity step
        
        dF1_dx_p=(F1(k+1)-F1(k))/dx;
        dU1_dt_p(k)=-dF1_dx_p;
        %momentum equation
        dF2_dx_p=(F2(k+1)-F2(k))/dx;
        dU2_dt_p(k)=-dF2_dx_p+J2(k);
        
        dF3_dx_p=(F3(k+1)-F3(k))/dx;
        dU3_dt_p(k)=-dF3_dx_p;
        
        U1(k)=(dU1_dt_p(k).*dt)+U1_old(k);
        U2(k)=(dU2_dt_p(k).*dt)+U2_old(k);
        U3(k) = (dU3_dt_p(k).*dt)+U3_old(k);
    end
    %applying boundry condition
    %input
    U1(1)=A(1);
    U2(1)=2*U2(2)-U2(3);
    U3(1)=U1(1)*((Tc(1)/(gamma-1)+((gamma/2)*Vc(1)^2)));
    %output
    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);    
    
    rhoc=U1./A;
    Tc=(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);
    
    %corrector method
    for k=2:n-1
        dA_dx(j)=(A(k)-A(k-1))/dx;
        dlogA_dx=(log(A(k))-log(A(k-1)))/dx;
        J2(k)=dlogA_dx*((gamma-1)/gamma)*(U3(k)-(0.5*gamma*U2(k)^2)/(U1(k)));
        
        %continuity step
        dF1_dx_c=(F1(k)-F1(k-1))/dx;
        dU1_dt_c(k)=-dF1_dx_c;
        
        %momentum equation
        dF2_dx_c=(F2(k)-F2(k-1))/dx;
        dU2_dt_c(k)=-dF2_dx_c +J2(k);
        
        %energy equation
        dF3_dx_c=(F3(k)-F3(k-1))/dx;
        dU3_dt_c(k)=-dF3_dx_c;
    end
    for pc=2:n-1
        dU1_dt_avg(pc)=0.5*(dU1_dt_c(pc)+dU1_dt_p(pc));
        dU2_dt_avg(pc)=0.5*(dU2_dt_c(pc)+dU2_dt_p(pc));
        dU3_dt_avg(pc)=0.5*(dU3_dt_c(pc)+dU3_dt_p(pc));
        
        U1(pc)=(dU1_dt_avg(pc)*dt)+U1_old(pc);
         U2(pc)=(dU2_dt_avg(pc)*dt)+U2_old(pc);
          U3(pc)=(dU3_dt_avg(pc)*dt)+U3_old(pc);
    end
    
    %boundry 
    %input
    U1(1)=A(1);
    U2(1)=2*U2(2)-U2(3);
    U3(1)=U1(1)*((Tc(1)/(gamma-1)+((gamma/2)*Vc(1)^2)));
    %input
    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);
    
    
    rhoc=U1./A;
    Vc=U2./U1;
    Tc=(gamma-1)*((U3./U1)-((gamma/2).*(Vc.^2)));
    
    pc=rhoc.*Tc;
    Mc=Vc./(Tc.^0.5);
    cons_mass_flow=rhoc.*Vc.*A;
    %variables at throat
    rho_throat(j)=rhoc(throat);
    T_throat(j)=Tc(throat);
    pc_throat(j)=pc(throat);
    Mc_throat(j)=Mc(throat);
    cons_mass_flow_throat(j)=cons_mass_flow(throat);
    Vc_throat(j)=Vc(throat);
    
    figure(2)
    if j==50
        plot(x,cons_mass_flow,'b','linewidth',1.5);
        hold on 
    elseif j==100
        plot(x,cons_mass_flow,'r','linewidth',1.5)
        hold on 
    elseif j==150
        plot(x,cons_mass_flow,'c','linewidth',1.5)
        hold on 
    elseif j==300
        plot(x,cons_mass_flow,'g','linewidth',1.5)
        hold on 
    elseif j==600
        plot(x,cons_mass_flow,'y','linewidth',1.5)
        hold on 
    elseif j==900
        plot(x,cons_mass_flow,'*','linewidth',1.5)
        hold on 
    elseif j==1400
        plot(x,cons_mass_flow,'m','linewidth',1.5)
        hold on 
    
    
    end
end
xlabel('Non dimensioanl distance through nozzle(x)','fontsize',10)
ylabel('Non dimensional mass flow rate through nozzle x','fontsize',10);
title('variation of mass flow rate (conservation form) along x at different time steps');
  legend({'50 time steps','100 time steps','150 time steps','300 time steps','600 time steps','900 time steps','1400 time steps'});
end

e plot for conservative:

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

Now suppose if the no of grid points is found to be 61 then,

At n=61