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

Non Conservation form Governing Equation for 1D Supersonic Nozzle

1) Continuity Equation

`del/(delt)(rhoA) +del/(delx)(rhoAV)=0`

2) Momentum Equation

`rho(delV)/(delt) +rhoV(delV)/(delx)=-(delp)/(delx)`

3) Energy Equation

`rhoc_v(delT)/(delt) +rhoVc_v(delT)/(delx)=-p(delV)/(delx)-pV(del(lnA))/(delx)`

We define non dimensional variables as follows

1) `T'=T/T_o`

2) `rho'=rho/rho_o`

3) `x'=x/L`

4) `a_o=sqrt(gammaRT_o)`

5) `V'=V/a_o`

where T_o, rho_o,a_o are the temperature, density and speed of sound at the reservoir

Writing the governing equation in terms of non dimensional temperature,density and velocity.

1) Continuity Equation

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

2) Momentum Equation

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

3) Energy Equation

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

Macormack Method

1) Predictor Step

In this step we solve the spatial derivatives in the governing equations as forward differences

a) Continuity Equation

`((delrho)/(delt))^t=-(rho_i^t)((V_(i+1)^t-V_i^t)/(dx))-rho_i^tV_i^t((lnA_(i+1)-lnA_i)/dx)-V_i^t(rho_(i+1)^t-rho_(i)^t)/dx`

b) Momentum Equation

`((delV)/(delt))^t=-V_i^t(V_(i+1)^t-V_i^t)/(Deltax) -1/gamma((T_(i+1)^t-T_i^t)/(Deltax) +(T_i^t/rho_i^t)(rho_(i+1)^t-rho_i^t)/(Deltax))`

c) Energy Equation

`((delT)/(delt))_i^t=-V_i^t((T_(i+1)^t-T_i^t)/(Deltax))-(gamma-1)T_i^t((V_(i+1)^t-V_i^t)/(Deltax) + V_i^t(lnA_(i+1)-lnA_(i))/(Deltax))`

Predicted Values can be obtained as follows

`barrho_i^(t+Deltat)=rho_i^(t) +((delrho)/(delt))_i^t (Deltat)`

`barV_i^(t+Deltat)=V_i^t +((delV)/(delt))_i^t (Deltat)`

`barT_i^(t+Deltat)=T_i^t +((delT)/(delt))_i^t (Deltat)`

Corrector Step

In this step we write the spatial derivatives in the governing equations with rearward differences using the predicted (barred quantities)

a) Continuity Equation

`((bar(delrho)/(delt)))_i^t=-(barrho_i^(t+Deltat))((barV_(i)^(t+Deltat)-barV_(i-1)^(t+Deltat))/(dx))-barrho_i^(t+Deltat)barV_i^(t+Deltat)((lnbarA_(i)-lnA_(i-1))/dx)-V_i^(t+Deltat)(rho_(i)^(t+Deltat)-rho_(i-1)^(t+Deltat))/dx`

b) Momentum Equation

`((bar(delV)/(delt))_i^t=-barV_i^(t+Deltat)(barV_(i)^(t+Deltat)-barV_(i-1)^(t+Deltat))/(Deltax) -1/gamma((barT_(i)^(t+DeltaT)-barT_(i-1)^(t+Deltat))/(Deltax) +(barT_i^(t+Deltat)/barrho_i^(t+Deltat))(barrho_(i)^(t+Deltat)-barrho_(i-1)^(t+Deltat))/(Deltax))`

c) Energy Equation

`((bar(delT))/(delt))_i^(t+Deltat)=-barV_i^(t+Deltat)((barT_(i)^(t+DeltaT)-barT_(i-1)^(t+Deltat))/(Deltax))-(gamma-1)barT_i^(t+Deltat)((barV_(i)^(t+Deltat)-barV_(i-1)^(t+Deltat))/(Deltax) + barV_i^(t+Deltat)(lnbarA_(i)-lnbarA_(i-1))/(Deltax))`

Computation of Average Time Derivatives

`((delrho)/(delt))_(av)=0.5*[((delrho)/(delt))_(predicter) +(bar(delrho)/(delt))_(crrectr)]`

`((delV)/(delt))_(av)=0.5*[((delV)/(delt))_(predicter) +(bar(delV)/(delt))_(crrectr)]`

`((delT)/(delt))_(av)=0.5*[((delT)/(delt))_(predicter) +(bar(delT)/(delt))_(crrectr)]`

Corrected Values of Flow field Variables

`rho_i^(t+Deltat)=rho_i^(t) +((delrho)/(delt))_(av) (Deltat)`

`V_i^(t+Deltat)=V_i^t +((delV)/(delt))_ (av)(Deltat)`

`T_i^(t+Deltat)=T_i^t +((delT)/(delt))_(av)(Deltat)`

Calculation of Time Step

We calculate 'Deltat _i^t'at all grid points i=1 to i=N and choose the minimum value for use

Deltat=minimum('Deltat_1^t','Deltat_2^t',..........'Deltat_N^t')

Boundary Conditions

1 Subsonic Inflow Boundary Conditions (Point 1)

V1 is the inflow condition, ie point 1 on the grid. V2 and V3 are the next interior nodes. V1 is calculated as follows

`V_1= V_2-(V_3-V_2)/(Deltax) Deltax`

'V_1=2V_2-V_3

'rho_1=1' (Fixed independent of time)

'T_1=1'(Fixed independent of time)

Supersonic outflow boundary (Point N)

`V_N=2V_(N-1)-V_(N-2)`

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

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

Nozzle Shape

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

`rho=1-0.3146x`

`T=1-0.2314x`

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

Results

1) Flow Variables for Non Conservative Equation

From the above plot we note that pressure,temperature and density decrease along +ve X direction as the flow transitions from sonic to supersonic across the throat. The transition from sonic to supersonic is indicated by positive Mach number.

1) Timewise Variation of Density, Pressure and Temperature and Mach Number at nozzle throat.

From the above plot we note that the largest change in flow property takes place in the initial stage of the solution after which the steady state value is approached. At early times time derivatives are large and oscillate in value. At later stage these derivatives grow small rapidly.

Variation of Mass flow rate with time

The time wise Variation can also be seen from the above plot which shows the non dimensionalized mass flow rates at various time steps. We note a drastic change in profiles in the initial stages of the solution i.e 0 to 500 timesteps. However after 500 time steps the profile begins to settle down eventually coming to to straight line across the plot. This corresponds to the steady state solution where in mass flow rate is constant.

Demonstration of Grid Independence:- We tabulate the value of primitive variables at the throat for 31 and 61 grid points

No of Grid Point	Rho	T	p	M
31	0.6387	0.8365	0.5342	0.9994
61	0.6373	0.8353	0.5323	0.9997

We note from the above table that the grid points do not change appreciably when the number of grid points change from 31 to 61. Therefore it can be said that the solution is grid independent at 31 points.

Conservative form of Governing Equations

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'])/(delt') +(del[rho'((e')/(gamma-1) +(gamma')/2V'^2)A'V' +p'A'V'])/(delx')=0`

The above equations can be written in generic form in terms of 3 Vectors namely

1. Solution Vector U

2. Flux Vector F

3. Source Vector J

We define the above Variables as below

Pure Form of Flux terms

The modified Governing Equation is as follows

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

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

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

From the definition of U1,U2 and U3, we have

`rho'=U_1/(A')`

`V'=U_2/U_1`

`T'=e'=(gamma-1)(U_3/U_1 -gamma/2V'^2)`

`p'=rho'T'`

Boundary Conditions

Inlet Conditions (i=1)

`U_(1(i=1)) =(rho'A')_(i=1) =A'_(i=1)=fixed value`

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

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

Outflow Boundary Conditions

`(U_1)_N=2(U_1)_(N-1) -(U_1)_(N-2)`

`(U_2)_N=2(U_2)_(N-1) -(U_2)_(N-2)`

`(U_3)_N=2(U_3)_(N-1) -(U_3)_(N-2)`

Initial Condtions

Variations for rho' and T'

For initial conditions only we define the Velocity as

Initial Conditions for U1,U2 and U3

Macormack Method is applied to the below equations in a similar manner as it was done for non conservation equations

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

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

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

Results

1) Variation of Flow variables for Conservative Equation

2) Timewise variation of density, pressure , Temperature and Mach Number at nozzle throat

Demonstration of Grid Independence;- We tabulate the value of primitive variables at the throat for 31 and 61 grid points

No of Grid Points	rho	T	p	M
31	0.6383	0.8351	0.533	0.9992
61	0.6513	0.8404	0.5474	0.98

There is a slight variation between the two results when compared to the grid independence test performed for non conservative form.

Comparison of Conservation and Non Conservation Form

1. Conservation form yields better mass flow distribution when compared to non conservation form. It can be seen that the system converges to a steady state much earlier when compared to non conservation form

2. Conservation Form requires marginally more computational effort. We need to convert primitive variables into Flux vector, Solution Vector and source vector.

3. Non conservative form yields better results in the grid independence test.

%Main Program

clc
clear all
close all

n=31;
x=linspace(0,3,n);
dx=x(2)-x(1);
gamma=1.4;

nt=1500
for i=1:nt
NT(i)=i;
end
courant=0.5;
anc=1+(2.2*(x-1.5).^2)
[rhonc,Tnc,Vnc,pnc,Mnc,Densnc,Tempnc,Pressnc, Machnc,mflnc]=non_conservative(nt,courant,gamma,dx,x,n);
[rhoc,Tc,vc,pc,Mc,Densc,Tempc,Pressc,Machc,mflc]=conservative(nt,courant,gamma,dx,x,n);

a=rhonc(15)
b=Tnc(15)
c=pnc(15)
d=Mnc(15)
e=rhoc(16)
f=Tc(16)
g=pc(16)
h=Mc(16)

   figure(1)
subplot(4,1,1)
plot(x,Mnc)
 title('Flow Variables for Non Conservative Equation')
xlabel('Grid points')
ylabel('Mach Number')
grid on
subplot(4,1,2)
plot (x,rhonc)
xlabel('Grid points')
ylabel('Density')
grid on
subplot(4,1,3)
plot (x,Tnc)
xlabel('Grid points')
ylabel('Temperature')
grid on
subplot(4,1,4)
plot (x,pnc)
xlabel('Grid points')
ylabel('Pressure')
grid on


 figure(2)

subplot(4,1,1)
plot(x,Mc)
xlabel('Grid points')
ylabel('Mach Number')
title('Flow Variables for Conservative Equation')
grid on
subplot(4,1,2)
plot (x,rhoc)
xlabel('Grid points')
ylabel('Density')
grid on
subplot(4,1,3)
plot (x,Tc)
xlabel('Grid points')
ylabel('Temperature')
grid on
subplot(4,1,4)
plot (x,pc)
xlabel('Grid points')
ylabel('Pressure')
grid on



figure(3)
subplot(4,1,1)
plot(NT,Densnc)
xlabel('No of Time Steps')
ylabel('density')
title('Timewise variation of density, pressure and temperature and Mach number at the nozzle throat-Non Conservative')
grid on
subplot(4,1,2)
plot (NT,Tempnc)
xlabel('No of Time Steps')
ylabel('Temperature')
grid on
subplot(4,1,3)
plot (NT, Pressnc)
xlabel('No of Time Steps')
ylabel('Pressure')
grid on
subplot(4,1,4)
plot (NT,Machnc)
xlabel('No of Time Steps')
ylabel('Mach No')
grid on

figure(4)
subplot(4,1,1)
plot(NT,Densc)
xlabel('No of Time Steps')
ylabel('density')
title('Timewise variation of density, pressure and temperature and Mach number at the nozzle throat-Conservative')
grid on
subplot(4,1,2)
plot (NT,Tempc)
xlabel('No of Time Steps')
ylabel('Temperature')
grid on
subplot(4,1,3)
plot (NT, Pressc)
xlabel('No of Time Steps')
ylabel('Pressure')
grid on
subplot(4,1,4)
plot (NT,Machc)
xlabel('No of Time Steps')
ylabel('Mach No')
grid on


figure(5)
plot(x,mflc(1,:),'k')
title('Variation of mass flow rate with time- Conservative')
hold on
plot(x,mflc(50,:),'k','LineStyle','--')
hold on
plot(x,mflc(100,:),'k','LineStyle','-.')
hold on
plot(x,mflc(150,:),'k','Marker','*')
hold on
plot(x,mflc(200,:),'k','Marker','o')
hold on
plot(x,mflc(500,:),'k','LineStyle','--','Marker','o')
hold on
plot(x,mflc(1000,:),'k','LineStyle','-','Marker','d')
legend('Timestep 0','Timestep 50','Timestep 100','Timestep 150','Timestep 200','Timestep 500','Timestep 1000')
xlabel ('Non Dimensionalized distance in x direction')
ylabel('Mass Flow rate')
figure(6)
plot(x,mflnc(1,:),'k')
title('Variation of mass flow rate with time- Non Conservative')

hold on
plot(x,mflnc(50,:),'k','LineStyle','--')

hold on
plot(x,mflnc(100,:),'k','LineStyle','-.')

hold on
plot(x,mflnc(150,:),'k','Marker','*')

hold on
plot(x,mflnc(200,:),'k','Marker','o')

hold on
plot(x,mflnc(500,:),'k','LineStyle','--','Marker','o')

hold on
plot(x,mflnc(1000,:),'k','LineStyle','-','Marker','d')
legend('Timestep 0','Timestep 50','Timestep 100','Timestep 150','Timestep 200','Timestep 500','Timestep 1000')
xlabel ('Non Dimensionalized distance in x direction')
ylabel('Mass Flow rate')

figure (7)
plot(x,rhonc,'Linestyle','-','Marker','o')
hold on
plot(x,rhoc,'Linestyle','-','Marker','d')
legend('Non Conservative','Conservative')

%Non Conservative Part Code
function [rhonc,Tnc,Vnc,pnc,Mnc,Densnc,Tempnc,Pressnc, Machnc,mflnc]=non_conservative(nt,courant,gamma,dx,x,n)



%initial profiles
rhonc=1-0.3146*x;
Tnc=1-0.2134*x;
Vnc=(0.1+1.09*x).*Tnc.^0.5;
a=1+2.2*(x-1.5).^2;




for k=1:nt
    dt=min(courant.*(dx./(Tnc.^0.5+Vnc)));
    rho_old=rhonc;
    T_old=Tnc;
    V_old=Vnc;
  for j=2:n-1
      %Predictor Step
      %Continuity Equation
      dVdx=(( Vnc(j+1)-Vnc(j))/dx);
      dlogadx=(log(a(j+1))-log(a(j)))/dx;
      drhodx=(rhonc(j+1)-rhonc(j))/dx;
      dTdx=(Tnc(j+1)-Tnc(j))/dx;
      
     Term1 =-rhonc(j)*(dVdx)  ;
     Term2=-rhonc(j)*Vnc(j)*(dlogadx);
     Term3=-Vnc(j)*(drhodx) ;
     drhodt_p(j)= Term1+Term2+Term3;
     
      %Momentum Equation
      Term4=-Vnc(j)*(dVdx)  ;
      Term5=-(1/gamma)*(dTdx);
      Term6=-(1/gamma)*(Tnc(j)/rhonc(j))*(drhodx);
      dVdt_p(j)=Term4+Term5+Term6;
     
      %Energy Equation
      Term7=-Vnc(j)*(dTdx)  ;
      Term8=-(gamma-1)*Tnc(j)*(dVdx);
      Term9=-(gamma-1)*Tnc(j)*Vnc(j)*(dlogadx);
      dTdt_p(j)= Term7+Term8+Term9;
      
      Vnc(j)=Vnc(j)+dVdt_p(j)*dt;
      rhonc(j)=rhonc(j)+ drhodt_p(j)*dt;
      Tnc(j)= Tnc(j)+ dTdt_p(j)*dt;
      
  end
  for j=2:n-1
      %Corrector Step
      %Continuity Equation
      dVdx=(( Vnc(j)-Vnc(j-1))/dx);
      dlogadx=(log(a(j))-log(a(j-1)))/dx;
      drhodx=(rhonc(j)-rhonc(j-1))/dx;
      dTdx=(Tnc(j)-Tnc(j-1))/dx;
      
     Term1 =-rhonc(j)*(dVdx)  ;
     Term2=-rhonc(j)*Vnc(j)*(dlogadx);
     Term3=-Vnc(j)*(drhodx) ;
     drhodt_c(j)= Term1+Term2+Term3;
     
      %Momentum Equation
      Term4=-Vnc(j)*(dVdx)  ;
      Term5=-(1/gamma)*(dTdx);
      Term6=-(1/gamma)*(Tnc(j)/rhonc(j))*(drhodx);
      dVdt_c(j)=Term4+Term5+Term6;
     
      %Energy Equation
      Term7=-Vnc(j)*(dTdx)  ;
      Term8=-(gamma-1)*Tnc(j)*(dVdx);
      Term9=-(gamma-1)*Tnc(j)*Vnc(j)*(dlogadx);
      dTdt_c(j)= Term7+Term8+Term9;
      
  end
  % Computing Average Time derivative
  drhodt=0.5*(drhodt_c+drhodt_p);
  dVdt= 0.5*(  dVdt_c+ dVdt_p);
   dTdt=0.5*(dTdt_c+ dTdt_p);
   
   for i=2:n-1
    % Final Solution Update
   Vnc(i)=V_old(i)+dVdt(i)*dt;
   rhonc(i)=rho_old(i)+drhodt(i)*dt;
   Tnc(i)=T_old(i)+dTdt(i)*dt;
   pnc=rhonc.*Tnc;
 Mnc=Vnc./(Tnc.^0.5);
   if i==16
        Densnc(k)=rhonc(i);
 Tempnc(k)=Tnc(i);
 Pressnc(k)=pnc(i);
 Machnc(k)=Mnc(i);
   end
   
   
   
   
   end
  %Apply Boundary Conditions
  Vnc(1)=2*Vnc(2)-Vnc(3);%Inlet
  
  %Outlet
  
 Vnc(n)=2*Vnc(n-1)-Vnc(n-2);
 rhonc(n)=2*rhonc(n-1)-rhonc(n-2);
 Tnc(n)=2*Tnc(n-1)-Tnc(n-2);
 
 pnc=rhonc.*Tnc;
 Mnc=Vnc./(Tnc.^0.5);
 mflnc(k,:)=rhonc.*Vnc.*a;

end
end

%Conservative

function [rhoc,Tc,vc,pc,Mc,Densc,Tempc,Pressc,Machc,mflc]=conservative(nt,courant,gamma,dx,x,n)

%initial profiles
a=1+(2.2*(x-1.5).^2);%area

for i=1:n
    if(x(i)>=0&&x(i)<=0.5)
        rhoc(i)=1;
        Tc(i)=1;
    elseif(x(i)>=0.5&&x(i)<=1.5)
        rhoc(i)=1-0.366*(x(i)-0.5);
        Tc(i)=1-0.167*(x(i)-0.5);
    elseif(x(i)>=1.5 && x(i)<=3.0)
        rhoc(i)=0.634-0.3879*(x(i)-1.5);
        Tc(i)=0.833-0.3507*(x(i)-1.5);
    end
end
vc=0.59./(rhoc.*a);
    
  


 u1=rhoc.*a;
   u2=rhoc.*a.*vc;
   u3=u1.*((Tc./(gamma-1))+((gamma/2).*(vc.^2)));

for k=1:nt
  u1_old=u1;
  u2_old=u2;
  u3_old=u3;
    dt=min(courant.*(dx./(Tc.^0.5+vc)));
  %flux terms 
  f1 = u2;
    f2 = ((u2.^2)./u1) + (((gamma-1)/gamma)*(u3 - ((gamma*0.5*(u2.^2))./u1)));
    f3 = ((gamma*u2.*u3)./u1) - ((gamma*0.5*(gamma-1)*(u2.^3))./(u1.^2));
  

  for j=2:n-1
      dadx=(a(j+1)-a(j))/dx ;
      j2=(rhoc(j)*Tc(j)*dadx)/gamma;
      
      %Predictor Step
      
      du1dt_p(j)=-(f1(j+1)-f1(j))/dx ;%Continuity Equation
      du2dt_p(j)=-(f2(j+1)-f2(j))/dx + j2;%Momentum Equation
      du3dt_p(j)=-(f3(j+1)-f3(j))/dx; %Energy Equation
      
      %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_p=u1./a;
  v_p=u2./u1;
  T_p=(gamma-1)*((u3./u1)-(gamma/2)*(v_p.^2));
 f1 = u2;
    f2 = ((u2.^2)./u1) + (((gamma-1)/gamma)*(u3 - ((gamma*0.5*(u2.^2))./u1)));
    f3 = ((gamma*u2.*u3)./u1) - ((gamma*0.5*(gamma-1)*(u2.^3))./(u1.^2));
  for j=2:n-1
       dadx=(a(j)-a(j-1))/dx ;
      j2=(rhoc(j)*Tc(j)*dadx)/gamma;
      
      %Corrector Step
            du1dt_c(j)=-(f1(j)-f1(j-1))/dx ;%Continuity Equation
      du2dt_c(j)= -(f2(j)-f2(j-1))/dx + j2;%Momentum Equation
      du3dt_c(j)=-(f3(j)-f3(j-1))/dx; %Energy Equation
     
      
  end
  % Computing Average Time derivative
  du1dt=0.5*(du1dt_p+du1dt_c);
  du2dt=0.5*(du2dt_p+du2dt_c);
  du3dt=0.5*(du3dt_p+du3dt_c);
   for i=2:n-1
    % Final Solution Update
    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
 u1(1)=rho_p(1)*a(1) ;
 u2(1)=2*u2(2)-u2(3);
 v_p(1)=u2(1)/u1(1);
 u3(1)=u1(1)*((T_p(1)/(gamma-1))+((gamma/2)*(v_p(1)^2)));
  
  %Outlet
  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)*(v_p.^2));
 mass(k,:)=rhoc.*vc.*a;
 pc=rhoc.*Tc;
 Mc=vc./(Tc.^0.5);
 mflc(k,:)=rhoc.*vc.*a;
   
        Densc(k)=rhoc(16);
 Tempc(k)=Tc(16);
 Pressc(k)=pc(16);
 Machc(k)=Mc(16);
   
 
end
end