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

1D Super sonic nozzle flow:

assumptions:

since the we are taking only one dimension we assume that the flow properties vary only in one direction that is xdirection.

assuming uniform flow properties across any given cross section .such flow is defined as quasi one dimensional flow

steady state:

the governing equations are

continuity:

ρ1⋅v1⋅a1=ρ2⋅v2⋅a2">ρ1⋅v1⋅a1=ρ2⋅v2⋅a2ρ1⋅v1⋅a1=ρ2⋅v2⋅a2

Momentum :

p1⋅a1+ρ1⋅v12⋅a1⋅+∫p⋅da=p2⋅a2+ρ2⋅a2⋅v22">p1⋅a1+ρ1⋅v12⋅a1⋅+∫p⋅da=p2⋅a2+ρ2⋅a2⋅v22p1⋅a1+ρ1⋅v12⋅a1⋅+∫p⋅da=p2⋅a2+ρ2⋅a2⋅v22

energy:

h1+v122=h2+v222">h1+v122=h2+v222h1+v122=h2+v222

up on solving the equations we get

area

(AA⋅)2=1M2[2γ+1(1+γ-12⋅M2)]γ+1γ-1">(AA⋅)2=1M2[2γ+1(1+γ−12⋅M2)]γ+1γ−1(AA⋅)2=1M2[2γ+1(1+γ-12⋅M2)]γ+1γ-1

pressure:

ppo=(1+γ-12⋅M2)-γγ-1">ppo=(1+γ−12⋅M2)−γγ−1ppo=(1+γ-12⋅M2)-γγ-1

density:

ρρ0=(1+γ-12⋅M2)-1γ-1">ρρ0=(1+γ−12⋅M2)−1γ−1ρρ0=(1+γ-12⋅M2)-1γ-1

temperature:

TTo=(1+γ-12⋅M2)-1">TTo=(1+γ−12⋅M2)−1TTo=(1+γ-12⋅M2)-1

solving by using Maccormack's technique it is a predictor-corrector method

it consists of the following steps

1.solving the governing equations

2.applying the predictor method by forward differencing

3.applying the corrector method by backword differencing

4.calculating the average time values

5.updating the solution and applying the boundary conditions

solving the 1D super sonic flow by using the non conservative form of governing equations:

we get the

continuity equation:

dρdt=-ρ⋅dvdx-ρ⋅v⋅d(Ina)dx-v⋅(dρdx)">dρdt=−ρ⋅dvdx−ρ⋅v⋅d(Ina)dx−v⋅(dρdx)dρdt=-ρ⋅dvdx-ρ⋅v⋅d(Ina)dx-v⋅(dρdx)

momentum equation:

dvdt=-v⋅dvdx-1γ⋅(dTdx+Tρ⋅dρdx)">dvdt=−v⋅dvdx−1γ⋅(dTdx+Tρ⋅dρdx)dvdt=-v⋅dvdx-1γ⋅(dTdx+Tρ⋅dρdx)

energy equation

dTdt=-v⋅dTdx-(γ-1)⋅T⋅[dvdx+v⋅d(Ina)dx]">dTdt=−v⋅dTdx−(γ−1)⋅T⋅[dvdx+v⋅d(Ina)dx]dTdt=-v⋅dTdx-(γ-1)⋅T⋅[dvdx+v⋅d(Ina)dx]

conservative form of governing equations

continuity :

du1dt=-dF1dx">du1dt=−dF1dxdu1dt=-dF1dx

momentum :

du2dt=-dF2dx+J2">du2dt=−dF2dx+J2du2dt=-dF2dx+J2

energy :

du3=-dF3dx">du3=−dF3dxdu3=-dF3dx

where

u1=ρ⋅a">u1=ρ⋅au1=ρ⋅a

u2=ρ⋅a⋅v">u2=ρ⋅a⋅vu2=ρ⋅a⋅v

u3=ρ⋅(Tγ-1+γ2⋅v2)⋅a">u3=ρ⋅(Tγ−1+γ2⋅v2)⋅au3=ρ⋅(Tγ-1+γ2⋅v2)⋅a

F1=ρ⋅a⋅v">F1=ρ⋅a⋅vF1=ρ⋅a⋅v

F2=ρ⋅a⋅v+(1γ)⋅p⋅a">F2=ρ⋅a⋅v+(1γ)⋅p⋅aF2=ρ⋅a⋅v+(1γ)⋅p⋅a

F3=u3⋅v+p⋅a⋅v">F3=u3⋅v+p⋅a⋅vF3=u3⋅v+p⋅a⋅v

J2=1γ⋅p⋅dadx">J2=1γ⋅p⋅dadxJ2=1γ⋅p⋅dadx

a=area

T=temperature

v=velocity

rho=density

p=pressure

code for conservative and non conservative forms with maccormack's method:

code for non consrvative form

function[t_n_throat,rho_throat,p_n_throat,M_n_throat,rho,t,v,m_n,M_n,p_n]=nonconservativenozzle(x,n,nt,dx)
%inputs
throat=((n-1)/2)+1;
m_a=0.578651;%mass flow rate
%calculations of initial profiles
rho=1-0.3146*x;
t=1-0.2314*x;
v=(0.1+1.09*x).*t.^0.5;
gamma=1.4;
%area
a=1+2.2*(x-1.5).^2;
%time steps
c=0.5;
dt=min(c*(dx./(sqrt(t)+v)))
%outer time loop
for k=1:nt
  rho_old=rho;
  v_old=v;
  t_old=t;
  %predictor method forward diferencing
   
  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;
    %continuty equation
    drhodt_p(j)=-rho(j)*dvdx-rho(j)*v(j)*dlogadx-v(j)*drhodx;
    %momentum
    dvdt_p(j)=-v(j)*dvdx-(1/gamma)*(dtdx+(t(j)/rho(j))*drhodx);
    %energy 
    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
 %correction method backword method
  for l=2:n-1
    dvdx=(v(l)-v(l-1))/dx;
    drhodx=(rho(l)-rho(l-1))/dx;
    dtdx=(t(l)-t(l-1))/dx;
    dlogadx=(log(a(l))-log(a(l-1)))/dx;
    %continuty equation 
    drhodt_c(l)=-rho(l)*dvdx-rho(l)*v(l)*dlogadx-v(l)*drhodx;
    %momentum
    dvdt_c(l)=-v(l)*dvdx-(1/gamma)*(dtdx+t(l)/rho(l)*drhodx);
    %energy 
    dtdt_c(l)=-v(l)*dtdx-(gamma-1)*t(l)*(dvdx+v(l)*dlogadx);
  end
  %compute the average time drivatives
  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
   %apply boundary conditions
   v(1)=2*v(2)-v(3);
  %outlet
  rho(n)=2*rho(n-1)-rho(n-2);
  v(n)=2*v(n-1)-v(n-2);
  t(n)=2*t(n-1)-t(n-2);
  m_n=rho.*v.*a; %mass flow rate
  p_n=rho.*t; % pressure
  M_n=v./sqrt(t); %macron number.
  rho_throat(k)=rho(throat);
  p_n_throat(k)=p_n(throat);
  M_n_throat(k)=M_n(throat);
  t_n_throat(k)=t(throat);
  figure(2)
  if k==50
    plot(x,m_n,'b','linewidth',1.5)
    hold on
    elseif k==100
       plot(x,m_n,'g','linewidth',1.5)
       hold on
     else if k==150
       plot(x,m_n,'r','linewidth',1.5)
       hold on
        elseif k==200
       plot(x,m_n,'y','linewidth',1.5)
       hold on
        elseif k==500
       plot(x,m_n,'k','linewidth',1.5)
       hold on
        elseif k==750
       plot(x,m_n,'m','linewidth',1.5)
       hold on
        elseif k==1400
       plot(x,m_n,'c','linewidth',1.5)
       plot(x,m_a,'0','linewidth',1.5)
       xlabel('nozile length')
       ylabel('mass flowrate')
       legend('time50','time100','time150','time200','time500','time750','time1400','exact')
       title('Non conservative mass flowrate along length of nozle at different times& exact value')
   end
 end
end
end

code for conservative form

function[T_throat,p_throat,M_throat,rho_throat,rho,T,v,m,M,p]=conservativenozzle(x,n,nt,dx)
%inputs
gamma=1.4;
m_a=0.578651;%mass flow rate
throat=((n-1)/2)+1;
%initial clculations
a=1+2.2*(x-1.5).^2;
for i=1:n
  if(x(i)>=0 && x(i)<=0.5)
     rho(i)=1.0;
     T(i)=1.0;
   elseif(x(i)>=0.5 && x(i)<=1.5)
   rho(i)=1.0-0.366*(x(i)-0.5);
   T(i)=1.0-0.167*(x(i)-0.5);
   elseif(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.590)./(rho.*a);
mexact=0.579;
u1=rho.*a;
u2=rho.*a.*v;
u3=rho.*((T./(gamma-1))+(gamma.*v.^2)/2).*a;
%time steps
c=0.5;
dt=min(c*(dx./(sqrt(T)+v)))
%outer time loop
for k=1:nt
   u1_old=u1;
   u2_old=u2;
   u3_old=u3;
   %predictor method forward differencing
   for i=1:n-1
  J1(i)=(1/gamma).*rho(i).*T(i).*(a(i+1)-a(i))/dx;
end
   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));
    for j=2:n-1
   du1dt_p(j)=-(f1(j+1)-f1(j))/(dx);
   du2dt_p(j)=-(f2(j+1)-f2(j))/(dx)+J1(j);
   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
 %decoding the primitive values
rho=u1./a;
v=u2./u1;
T=(gamma-1)*((u3./u1)-(v.^2*gamma)/2);
u1=rho.*a;
u2=rho.*a.*v;
u3=u1.*((T./(gamma-1))+(gamma.*v.^2)/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));
   for i=2:n-1
  J2(i)=(1/gamma).*rho(i).*T(i).*(a(i)-a(i-1))/dx;
  end
   % corrector method backward differencing
 for i=2:n-1 
   du1dt_c(i)=-(f1(i)-f1(i-1))/(dx);
   du2dt_c(i)=-(f2(i)-f2(i-1))/(dx)+J2(i);
   du3dt_c(i)=-(f3(i)-f3(i-1))/(dx);
end
%compute average time derivative
du1dt_av=0.5*(du1dt_p+du1dt_c);
du2dt_av=0.5*(du2dt_p+du2dt_c);
du3dt_av=0.5*(du3dt_p+du3dt_c);
%final update of solution
for j=2:n-1
  u1(j)=u1_old(j)+du1dt_av(j)*dt;
  u2(j)=u2_old(j)+du2dt_av(j)*dt;
  u3(j)=u3_old(j)+du3dt_av(j)*dt;
end
%computing the values;
rho=u1./a;
v=u2./u1;
T=(gamma-1)*((u3./u1)-(v.^2*gamma)/2);
%initial conditions
u1(1)=rho(1).*a(1);
u2(1)=2*u2(2)-u2(3);
v(1)=u2(1)./u1(1);
u3(1)=u1(1).*(((T(1))./(gamma-1))+(gamma/2).*v(1).^2);
%outlet boundary condiions
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);

%Decoding the values;
rho=u1./a;
v=u2./u1;
T=(gamma-1)*((u3./u1)-(v.^2*gamma)/2);
%calculating the values
 m=rho.*v.*a;%mass flow rate
 p=rho.*T; % pressure
 M=v./sqrt(T); %macron number 
 rho_throat(k)=rho(throat);
 p_throat(k)=p(throat);
 M_throat(k)=M(throat);
 T_throat(k)=T(throat);
 figure(1)
  if k==50
    plot(x,m,'b','linewidth',1.5)
    hold on
    elseif k==100
       plot(x,m,'g','linewidth',1.5)
       hold on
     else if k==150
       plot(x,m,'r','linewidth',1.5)
       hold on
        elseif k==200
       plot(x,m,'y','linewidth',1.5)
       hold on
        elseif k==500
       plot(x,m,'k','linewidth',1.5)
       hold on
        elseif k==750
       plot(x,m,'m','linewidth',1.5)
       hold on
        elseif k==1400
       plot(x,m,'c','linewidth',1.5)
        plot(x,m_a,'0','linewidth',1.5)
       xlabel('nozile length')
       ylabel('mass flowrate')
       legend('time50','time100','time150','time200','time500','time750','time1400','exact')
       title('conservative mass flowrate along length of nozle at different times& exact value')
   end
 end
end
end

main program

clear all
close all
clc
n=61;
x=linspace(0,3,n);
dx=x(2)-x(1)
nt=1400;
%analytical results
A=1;%area
p_a=0.528;%pressure
rho_a=0.63393;%density
t_a=0.8333;%temperature
v_a=0.9128;% velocity 
m_a=0.578651;%mass flow rate
M_a=1;%macnumber
tic;
[T_throat,p_throat,M_throat,rho_throat,rho_c,T_c,v_c,m_c,M_c,p_c]=conservativenozzle(x,n,nt,dx);
toc;
tic;
[t_n_throat,rho_n_throat,p_n_throat,M_n_throat,rho_n,t_n,v_n,m_n,M_n,p_n]=nonconservativenozzle(x,n,nt,dx);
toc;
figure(3)
subplot(4,1,1)
plot(x,rho_c)
ylabel('density')
title('conservativeform,density,Temperature,velocity,Macnumber,variation along length of nozele')
subplot(4,1,2)
plot(x,T_c)
ylabel('Temperature')
subplot(4,1,3)
plot(x,v_c)
ylabel('velocity')
subplot(4,1,4)
plot(x,M_c)
ylabel('Macnumber')
xlabel('length of nozele')
figure(4)
subplot(4,1,1)
plot(x,rho_n)
ylabel('density')
title('nonconservativeform,density,Temperature,velocity,Macnumber,variation along length of nozele')
subplot(4,1,2)
plot(x,t_n)
ylabel('Temperature')
subplot(4,1,3)
plot(x,v_n)
ylabel('velocity')
subplot(4,1,4)
plot(x,M_n)
ylabel('Macnumber')
xlabel('length of nozele')
figure(5)
plot(x,m_c,'r','linewidth',1.5)
hold on
plot(x,m_n,'g','linewidth',1.5)
hold on
plot(x,m_a,'b','linewidth',1.5)
xlabel('length of nozele')
ylabel('mass floe rate')
title('mass flow rate in conservative & nonconservativeform & analytic')
legend('conservativeform','nonconservativeform','analytic')
figure(6)
subplot(4,1,1)
plot(1:nt,rho_throat,'-')
hold on
plot(1:nt,rho_n_throat)
plot(1:nt,rho_a)
xlabel('time step')
ylabel('density')
legend('conservativeform','nonconservativeform','analytic')
subplot(4,1,2)
plot(1:nt,p_throat,'-')
hold on
plot(1:nt,p_n_throat)
plot(1:nt,p_a)
xlabel('time step')
ylabel('pressure')
legend('conservativeform','nonconservativeform','analytic')
subplot(4,1,3)
plot(1:nt,T_throat,'-')
hold on
plot(1:nt,t_n_throat)
plot(1:nt,t_a)
xlabel('time step')
ylabel('Temperature')
legend('conservativeform','nonconservativeform','analytic')
subplot(4,1,4)
plot(1:nt,M_throat,'-')
hold on
plot(1:nt,M_n_throat)
plot(1:nt,M_a)
xlabel('time step')
ylabel('Macnumber')
legend('conservativeform','nonconservativeform','analytic')

plots for different time steps:

n=31

n=61

Time step variations:

n=31

n=61

convergence:

At time step 500 the values of the mass flow rate are not changing more and up on increasing the time steps at 750 time steps the values of mass flow rate are constant and almost same as the exact value and even after increasing the time steps there is no change in the values as we have seen for the 1400 time step.from the plots we observe that the conservative form at 200 time steps are more nearer to the exact solution when compared with the non conservative part.hence the conservative form convergeces more faster than the non conservative form and the minimum number of time steps required are 200 time steps for conservative.and for non conservative

mass flowrate

n=31

n=61

mass flow rate:

From the graph we observe that the conservative form is more stable than the non conservative form.we observe that the mass flow rate in non conservative form is unstady at the starting it is increasing and then it decreasing at the throat and than it is increasing at the end of the nozzle.

in conservative flow the mass flow remains almost constant in the nozzle and this gives use the better results

length vs properties:

n=31

n=61

grid dependence:

The plots that are drawn for both n=31&n=61 are same and there is no much change in the values hence we can say that the 31 grid points can be the dependence point .since after the 31 grid points even if we increase the grid number there is no change in the values.

faster convergence:

conservative form elapsed time is 7.2seconds

nonconservative form elapsed time is 10.5seconds

by this we can say that conservative form convergeces faster and at less time steps than the non conservative form.