Week 7 - Simulation of a 1D Super-sonic nozzle flow simulation using Macormack Method
AIM- To simulate 1D-super sonic flow using macormack method for both conservative and non-conservative form of governing equations.
EXPLANATION-
we are considering a nozzle which is firstly converging then diverging. The inlet of nozzle is attached to reservoir and due to which area at inlet will tends to infinity and hence velocity will be small at inlet.The flow will be subsonic at convergent or reservoir section,then sonic at throat section and supersonic at outlet or divergent section as shown in figure 1.
Fig 1.
We assme that at a given section,where cross-sectional area is A ,the flow properties are uniform across that section. Hence,although the area of nozzle changes as a function of distance along the nozzle,x, and therefore,in reality the flow field is 2-D,we make the assumption that flow properties vary only with x:this is tantamount to assuming uniform flow properties across any given cross-section. Such flow is defined as quasi-one dimensional flow.
The governing continuity,momentum and energy equations for this quasi-one dimensional,steady,isentropic flow can be expressed ,respectively ,as
Non-conservative continuity equation:-
`(dp')/(dt')=-p'*(dv')/(dx')-p'*v'*(d(ln A'))/dx'-v'*(dp')/(dx')`
Non-conservative momentum equation:-
`(dv')/(dt')=-v'*(dv')/(dx')-(1/gamma)*((dt')/(dx')+((T')/(p'))*(dp')/(dx'))`
Non-conservative energy equation:-
`(dT')/(dt')=-v'*(dT')/(dx')-(gamma-1)*T'*((dv')/(dx')+v'*((d(ln A'))/(dx'))`
Conservative continuity equation:-
`(d(p'A'))/dt+(d(p'A'v'))/(dx')=0`
conservative momentum equation:-
`(d(p'A'v'))/(dt')+d((p'A'(v')^2+(1/gamma)*p'A')/(dx')=(1/gamma)*p'(dA')/(dx')`
Conservative energy equation:-
`d((p'((e')/(gamma-1)+(gamma/2)*(V')^2)A'))/(dt')+ d((p'((e')/(gamma-1)+(gamma/2)*(V')^2)*V'A'+p'A'V'))/(dx')=0`
where `gamma=`ratio of specific heats=Cp/Cv. For air at standard conditions,`gamma=1.4`
These equations are solved by macormack method using time marching process.
Macormack method:-
Let f(x) is the equation to be solved. These are the steps to be followed:-
1. Initial profile:-
These are the guess values for all the flow field variables which are velocity,density and temperature.
2. Predictor step:-
`((df)/(dt))_p=(f(i+1)-f(i))/(Delta(x))`
3. Solution update:-
`(f)_p=f(Initial)+((df)/(dt))_p*(Delta t)`
4. Corrector step:-
`((df_p)/(dt))_c=(f(i)-f(i-1))/(Delta x)`
`((df)/(dt))_avg=0.5(((df)/(dt))_p+((df_p)/(dt))_c)`
5. Final solution update:-
`f_(sol)=f_(Initial)+((df)/(dt))_avg*(Delta t)`
Code:
clear all
close all
clc
L=3,
n=[31 62 93]
nt=1:1400;
cfl=0.5;
for i=1:length(n)
x=linspace(0,L,n(i));
[mfr_nc,mfr_throat_nc]=non_conservative(L,n(i),nt,cfl);
[mfr_c,mfr_throat_c]=conservative(L,n(i),nt,cfl);
figure(7)
subplot(2,1,1)
plot(nt,mfr_throat_c,'b',nt,mfr_throat_nc,'r','linewidth',1)
xlabel('Time steps(nt)')
ylabel('Mass flow rate')
legend('conservative form','Non conservative form')
title('Mass flow rate using "Non-conservative form vs conservative form" at throat point i=16')
subplot(2,1,1)
plot(x,mfr_c,'b',x,mfr_nc,'r','linewidth',1)
xlabel('length of nozzle')
ylabel('Mass flow rate')
legend('conservative form','Non conservative form')
title('Mass flow rate using "Non-conservative form vs conservative form" in the nozzle')
figure(8)
subplot(2,1,1)
plot(nt,mfr_throat_c,'linewidth',1)
hold on
xlabel('Time steps(nt)')
ylabel('Mass flow rate')
legend('@31 nodes','@62 nodes','@93 nodes')
title('conservative form : Grid independency test at throat point i=16')
subplot(2,1,2)
plot(nt,mfr_nc,'linewidth',1)
hold on
xlabel('Time steps(nt)')
ylabel('Mass flow rate')
legend('@31 nodes','@62 nodes','@93 nodes')
title('Non-conservative form : Grid independency test at throat point i=16')
figure(9)
subplot(2,1,1)
plot(x,mfr_c,'linewidth',1)
hold on
xlabel('lenth of nozzle')
ylabel('Mass flow rate')
legend('@31 nodes','@62 nodes','@93 nodes')
title('conservative form : Grid independency test along the nozzle')
subplot(2,1,2)
gnuplot(nt,mfr_nc,'linewidth',1)
hold on
xlabel('length of nozzle')
ylabel('Mass flow rate')
legend('@31 nodes','@62 nodes','@93 nodes')
title('Non-conservative form : Grid independency test along the nozzle')
end
Non-conservative function:-
function [mfr_nc,mfr_throat_nc]=non_conservative(L,n,nt,cfl)
x=linspace(0,L,n);
dx=x(2)-x(1);
%initial conditions rho,a,v,T,mfr,M,p
rho=1-0.3146*x;
T=1-0.2314*x;
v=(0.1+1.09*x).*T.^0.5;
a=1+2.2*(x-1.5).^2;
M=v./(T.^0.5);
mfr=rho.*a.*v;
p=rho.*T;
g=1.4;
for z=1:31
u=(T(z))^(1/2);
dt(z)=cfl*dx/(u+v(z));
endfor
dt=min(dt);
for i=1:length(nt)
%initial values
rho_old=rho;
T_old=T;
v_old=v;
for j=2:n-1
%predictor step
drho_dx=(rho(j+1)-rho(j))/dx;
dv_dx=(v(j+1)-v(j))/dx;
dT_dx=(T(j+1)-T(j))/dx;
da_dx=(log(a(j+1))-log(a(j)))/dx;
%continuity equation
drho_dt_p(j)=-rho(j)*dv_dx-rho(j)*v(j)*da_dx-v(j)*drho_dx;
%Momentum equation
dv_dt_p(j)=-v(j)*dv_dx-(1/g)*(dT_dx+(T(j)/rho(j))*drho_dx);
%energy equation
dT_dt_p(j)=-v(j)*dT_dx-(g-1)*T(j)*(dv_dx+v(j)*da_dx);
%solution update
rho(j)=rho_old(j)+(drho_dt_p(j)*dt);
v(j)=v_old(j)+(dv_dt_p(j)*dt);
T(j)=T_old(j)+(dT_dt_p(j)*dt);
%corrector step
drho_dx=(rho(j)-rho(j-1))/dx;
dv_dx=(v(j)-v(j-1))/dx;
dT_dx=(T(j)-T(j-1))/dx;
da_dx=(log(a(j))-log(a(j-1)))/dx;
%continuity equation
drho_dt_c(j)=-rho(j)*dv_dx-rho(j)*v(j)*da_dx-v(j)*drho_dx;
%Momentum equation
dv_dt_c(j)=-v(j)*dv_dx-(1/g)*(dT_dx+(T(j)/rho(j))*drho_dx)
%energy equation
dT_dt_c(j)=-v(j)*dT_dx-(g-1)*T(j)*(dv_dx+v(j)*da_dx);
end
%compute the average time derivative
drho_dt_avg=0.5*(drho_dt_p+drho_dt_c);
dv_dt_avg=0.5*(dv_dt_p+dv_dt_c);
dT_dt_avg=0.5*(dT_dt_p+dT_dt_c);
%final solution update
for k=2:n-1
rho(k)=rho_old(k)+drho_dt_avg(k)*dt;
v(k)=v_old(k)+dv_dt_avg(k)*dt;
T(k)=T_old(k)+dT_dt_avg(k)*dt;
endfor
%apply the 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);
% Final updates of primitives
p=rho.*T;
M=v./T.^0.5;
error(k)=abs(v(k)-v(k-1));
mfr=rho.*a.*v;
%to find values at throat point
T_throat(i)=T(16);
rho_throat(i)=rho(16);
v_throat(i)=v(16);
M_throat(i)=M(16);
p_throat(i)=p(16);
error_throat(i)=error(16);
mfr_throat(i)=mfr(16);
if n==31
for ii=[1:50:201 700]
if ii==1
figure(1)
plot(x,mfr,'linewidth',1)
hold on
title('Non-conservative form-Variation of Mass flow rate distribution inside the nozzle at different time steps during time marching process')
xlabel('length of nozzle')
ylabel('Mass flow rate')
elseif i==1400
plot(x,mfr,'*r')
legend('@1','@51','@101','@151','@201','@700','@1400')
end
end
end
end
if n==31
figure(2)
plot(x,v,'b',x,p,'r',x,T,'k',x,rho,'m',x,M,'g',x,mfr,'c')
ylim([-0.1 3.7])
title(sprintf('Time step=%d',i))
xlabel('length of nozzle')
ylabel('Non dimensional parameters')
legend('Velocity','pressure','temperature','density','mach number','mass flow rate','location','northwest')
title('Non-conservative form-steady-state distribution of primitive variables inside the nozzle')
figure(3)
plot(nt,rho_throat)
hold on
plot(nt,v_throat,'linewidth',1)
plot(nt,p_throat,'linewidth',1)
plot(nt,T_throat,'linewidth',1)
plot(nt,M_throat,'linewidth',1)
plot(nt,mfr_throat,'linewidth',1)
xlabel('Time steps(nt)')
ylabel('Non dimensional parameters')
legend('Density','Velocity','pressure','temperature','mach number','mass flow rate')
title('Non conservative form-Time wise variation of primitive variables at throat point i=16')
end
mfr_nc=mfr;
mfr_throat_nc=mfr_throat;
end
Conservative function:-
function [mfr_c,mfr_throat_c]=conservative(L,n,nt,cfl)
x=linspace(0,L,n);
dx=x(2)-x(1);
%initial conditions rho,a,v,T,mfr,M,p
for ii=1:length(x)
if (x(ii)<=0.5)
rho(ii)=1;
T(ii)=1;
elseif (x(ii)>0.5)&(x(ii)<=1.5)
rho(ii)=1-0.366*(x(ii)-0.5);
T(ii)=1-0.617*(x(ii)-0.5);
elseif(x(ii)>1.5)&(x(ii)<=3.5)
rho(ii)=0.634-0.3879*(x(ii)-1.5);
T(ii)=0.833-0.3507*(x(ii)-1.5);
end
end
a=1+2.2*(x-1.5).^2;
g=1.4;
%solution vectors
u1=rho.*a;
u2=0.59*ones(1,n);
v=u2./u1;
u3=u1.*((T./(g-1))+(g/2).*v.^2);
M=v./(T.^0.5);
mfr=rho.*a.*v;
p=rho.*T;
for z=1:31
u=(T(z))^(1/2);
dt(z)=cfl*dx/(u+v(z));
endfor
dt=min(dt);
for i=1:length(nt)
%initial values
u1_old=u1;
u2_old=u2;
u3_old=u3;
for j=2:n-1
%flux vectors
f1=u2;
f2=((u2.^2)./u1)+((g-1)/g)*(u3-((g-2)*(u2.^2)./u1));
f3=((g.*u2.*u3)./u1)-((g*(g-1)/2)*(u2.^3)./(u1.^2));
%source terminal_size
j1=(1/g)*rho(j)*T(j)*((a(j+1)-a(j))/dx);
%predictor step
du1_dt_p(j)=-(f1(j+1)-f1(j))/dx; %continuity
du2_dt_p(j)=-(f2(j+1)-f2(j))/dx+j1; %momentum
du3_dt_p(j)=-(f3(j+1)-f3(j))/dx; %energy
%updating solution vectors
u1(j)=u1_old(j)+du1_dt_p(j)*dt;
u2(j)=u2_old(j)+du2_dt_p(j)*dt;
u3(j)=u3_old(j)+du3_dt_p(j)*dt;
%updating primitives
rho=u1./a;
T=(g-1).*(u3./u1-(g-2).*(u2./u1).^2);
%updating flux terms
f1=u2;
f2=((u2.^2)./u1)+((g-1)/g)*(u3-((g-2)*(u2.^2)./u1));
f3=((g.*u2.*u3)./u1)-((g*(g-1)/2)*(u2.^3)./(u1.^2));
%source term for corrector step
j2=(1/g)*rho(j)*T(j)*((a(j)-a(j-1))/dx);
%corrector step
du1_dt_c(j)=-(f1(j)-f1(j-1))/dx; %continuity
du2_dt_c(j)=-(f2(j)-f2(j-1))/dx+j2; %momentum
du3_dt_c(j)=-(f3(j)-f3(j-1))/dx; %energy
end
%compute the average time derivative
du1_dt_avg=0.5*(du1_dt_p+du1_dt_c);
du2_dt_avg=0.5*(du2_dt_p+du2_dt_c);
du3_dt_avg=0.5*(du3_dt_p+du3_dt_c);
%final solution update
for k=2:n-1
u1(k)=u1_old(k)+du1_dt_avg(k)*dt;
u2(k)=u2_old(k)+du2_dt_avg(k)*dt;
u3(k)=u3_old(k)+du3_dt_avg(k)*dt;
endfor
%apply the boundary conditions
%inlet
u1(1)=a(1);
u2(1)=2*u2(2)-u2(3);
u3(1)=u1(1)*(T(1)/(g-1)+(g*(u2(1)/u1(1))^2)/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);
% Final updates of primitives
rho=u1./a;
v=u2./u1;
T=(g-1).*(u3./u1-(g-2).*(u2./u1).^2);
p=rho.*T;
M=v./T.^0.5;
error(k)=abs(v(k)-v(k-1));
mfr=u2;
%to find values at throat point
T_throat(i)=T(16);
rho_throat(i)=rho(16);
v_throat(i)=v(16);
M_throat(i)=M(16);
p_throat(i)=p(16);
error_throat(i)=error(16);
mfr_throat(i)=mfr(16);
if n==31
for ii=[1:50:201 700]
if ii==1
figure(4)
plot(x,mfr,'linewidth',1)
hold on
title('conservative form-Variation of Mass flow rate distribution inside the nozzle at different time steps during time marching process')
xlabel('length of nozzle')
ylabel('Mass flow rate')
elseif i==1400
plot(x,mfr,'*r')
legend('@1','@51','@101','@151','@201','@700','@1400')
end
end
end
end
if n==31
figure(5)
plot(x,v,'b',x,p,'r',x,T,'k',x,rho,'m',x,M,'g',x,mfr,'c')
ylim([-0.1 3.7])
title(sprintf('Time step=%d',i))
xlabel('length of nozzle')
ylabel('Non dimensional parameters')
legend('Velocity','pressure','temperature','density','mach number','mass flow rate','location','northwest')
title('conservative form-steady-state distribution of primitive variables inside the nozzle')
figure(6)
plot(nt,rho_throat)
hold on
plot(nt,v_throat,'linewidth',1)
plot(nt,p_throat,'linewidth',1)
plot(nt,T_throat,'linewidth',1)
plot(nt,M_throat,'linewidth',1)
plot(nt,mfr_throat,'linewidth',1)
xlabel('Time steps(nt)')
ylabel('Non dimensional parameters')
legend('Density','Velocity','pressure','temperature','mach number','mass flow rate')
title('conservative form-Time wise variation of primitive variables at throat point i=16')
end
mfr_c=mfr;
mfr_throat_c=mfr_throat;
end