Week 7 - Simulation of a 1D Super-sonic nozzle flow simulation using Macormack Method
AIM:
1. Steady-state distribution of primitive variables inside the nozzle
2. Time-wise variation of the primitive variables
3. Variation of Mass flow rate distribution inside the nozzle at different time steps during the time-marching process
4. Comparison of Normalized mass flow rate distributions of both forms
To perform a grid independence test on the solution
MAIN CODE
clear all
close all
clc
%inputs
n=31;
N=linspace(1,1400,1400);
x=linspace(0,3,n);
dx=x(2)-x(1);
a=1+2.2*(x-1.5).^2;
gamma=1.4;
nt=1400;
[x1,a1,v1,rho1,t1,mass1,Mach1]=nonconservativequasi(n,x,dx,a,gamma,N,nt)
[x2,a2,v2,rho2,t2,mass2,Mach2]=conservative(n,x,dx,a,gamma,N,nt)
figure(9);
plot(x1,mass1);
hold on
plot(x2,mass2);
hold on
title(sprintf('Comparison of mass variation of x-Nonconservative and conservative'))
xlabel('time step')
ylabel('Mass flow rate')
legend('mass1','mass2')
axis([0 3 0.55 0.65])
grid on
pause(0.03)
x1=x1';
rho1=rho1';
rho2=rho2';
v1=v1';
v2=v2';
t1=t1';
t2=t2';
Mach1=Mach1';
Mach2=Mach2';
tis=table(x1,Mach1,Mach2,rho1,rho2,t1,t2,v1,v2)NON CONSERVATIVE FORM
function [x,a,v,rho,t,mass,mach]=nonconservativequasi(n,x,dx,a,gamma,N,nt)
%inputs
x=linspace(0,3,n);
dx=x(2)-x(1);
dx=round(dx,3);
rho=1-(0.3146*x);
rho=round(rho,3);
t=1-0.2134*x;
t=round(t,3);
v=(0.1+1.09*x).*t.^0.5;
v=round(v,3);
a=1+2.2*(x-1.5).^2;
a=round(a,3);
gamma=1.4;
c=0.5;
%time step
for i=1:n
dt(i)=c*dx/(t(i)^0.5+v(i));
end
dt=min(dt);
dt=round(dt,3);
%outer loop
for k=1:nt
%predictor method
vold=v;
rhoold=rho;
told=t;
vold=round(vold,3);
rhoold=round(rhoold,3);
told=round(told,3);
%continnuity eq
for j=2:n-1
%drho_dt_p(j)=-rho(j)*(v(j+1)-v(j))/dx-rho(j)*v(j)*(log(a(j+1))-log(a(j)))/dx-v(j)*(rho(j+1)-rho(j))/dx;
%forward defferencing method
dvdx=(v(j+1)-v(j))/dx;
dvdx=round(dvdx,3);
drhodx=(rho(j+1)-rho(j))/dx;
drhodx=round(drhodx,3);
dtdx=(t(j+1)-t(j))/dx;
dtdx=round(dtdx,3);
dlogadx=(log(a(j+1))-log(a(j)))/dx;
dlogadx=round(dlogadx,3);
%continuity eq
drhodt_p(j)=-rho(j)*dvdx-rho(j)*v(j)*dlogadx-v(j)*drhodx;
drhodt_p(j)=round(drhodt_p(j),3);
%momentum eq
dvdt_p(j)=-v(j)*dvdx-(1/gamma)*(dtdx+(t(j)/rho(j)*drhodx));
dvdt_p(j)=round(dvdt_p(j),3);
%energy equation
dtdt_p(j)=-v(j)*dtdx-(gamma-1)*t(j)*(dvdx+v(j)*dlogadx)';
dtdt_p(j)=round(dtdt_p(j),3);
%solution update
v(j)=v(j)+dvdt_p(j)*dt;
rho(j)=rho(j)+drhodt_p(j)*dt;
t(j)=t(j)+dtdt_p(j)*dt;
v=round(v,3);
rho=round(rho,3);
t=round(t,3);
end
%corrector method
for j=2:n-1
%backward differencing method
dvdx=(v(j)-v(j-1))/dx;
dvdx=round(dvdx,3);
drhodx=(rho(j)-rho(j-1))/dx;
drhodx=round(drhodx,3);
dtdx=(t(j)-t(j-1))/dx;
dtdx=round(dtdx,3);
dlogadx=(log(a(j))-log(a(j-1)))/dx;
dlogadx=round(dlogadx,3);
%continuity eq
drhodt_c(j)=-rho(j)*dvdx-rho(j)*v(j)*dlogadx-v(j)*drhodx;
drhodt_c(j)=round(drhodt_c(j),3);
%momentum eq
dvdt_c(j)=-v(j)*dvdx-(1/gamma)*(dtdx+(t(j)/rho(j)*drhodx));
dvdt_c(j)=round(dvdt_c(j),3);
%energy equation
dtdt_c(j)=-v(j)*dtdx-(gamma-1)*t(j)*(dvdx+v(j)*dlogadx)';
dtdt_c(j)=round(dtdt_c(j),3);
end
%compute average time derivative
drhodt=0.5*(drhodt_p+drhodt_c);
drhodt=round(drhodt,3);
dvdt=0.5*(dvdt_p+dvdt_c);
dvdt=round(dvdt,3);
dtdt=0.5*(dtdt_p+dtdt_c);
dtdt=round(dtdt,3);
%final sol
for i=2:n-1;
rho(i)=rhoold(i)+drhodt(i)*dt;
v(i)=vold(i)+dvdt(i)*dt;
t(i)=told(i)+dtdt(i)*dt;
rho=round(rho,3);
v=round(v,3);
t=round(t,3);
end
%apply bound cond
%inlet cond
v(1)=2*v(2)-v(3);
v=round(v,3);
%outlet cond
v(n)=2*v(n-1)-v(n-2);
v=round(v,3);
rho(n)=2*rho(n-1)-rho(n-2);
rho=round(rho,3);
t(n)=2*t(n-1)-t(n-2);
t=round(t,3);
figure(1);
subplot(3,1,1)
plot(x,v);
title(sprintf('Steady-state distribution of primitive variables of NonConservative %d',k))
xlabel('x')
ylabel('v')
axis([0 3 0 2])
grid on
subplot(3,1,2)
plot(x,t);
xlabel('x')
ylabel('t')
axis([0 3 0 2])
grid on
subplot(3,1,3)
plot(x,rho);
xlabel('x')
ylabel('rho')
axis([0 3 0 2])
grid on
%pause(0.03)
m=v./(t.^0.5);
figure(2)
plot(x,m);
xlabel('x')
ylabel('mach')
axis([0 3 0 5])
grid on
title(sprintf('Steady-state distribution of primitive variables of NonConservative %d',k))
pause(0.03)
ax=a(16);
p=rho.*t;
rhoo=rho(1);
ao=(t(1))^0.5;
to=t(1);
po=p(1);
x=x';
a=a';
v=v';
p=p';
rho=rho';
t=t';
m=m';
arearatio=a/ax;
densityratio=rho/rhoo;
velocitybyao=v/ao;
tempratio=t/to;
pressratio=p/po;
mach=m;
mass=rho.*a.*v;
figure(3)
if(k==50)
plot(x,mass,'o');
xlabel('x')
ylabel('massflow rate')
grid on
hold on
elseif(k==100)
plot(x,mass,'r');
elseif(k==200)
plot(x,mass,'g');
elseif(k==300)
plot(x,mass,'b');
elseif(k==400)
plot(x,mass,'c');
elseif(k==500)
plot(x,mass,'m');
elseif(k==600)
plot(x,mass,'y');
elseif(k==700)
plot(x,mass,'k');
elseif(k==1400)
plot(x,mass,'*');
title(sprintf('Variation of mass flow rate with x-NonConservative'))
legend('50','100','200','300','400','500','600','700','1400')
end
T1= table(x,a,v,rho,t);
densityratio16(k)=densityratio(16);
tempratio16(k)=tempratio(16);
pressratio16(k)=pressratio(16);
end
densityratio16=densityratio16';
tempratio16=tempratio16';
pressratio16=pressratio16';
figure(4);
subplot(3,1,1)
plot(N,densityratio16);
title(sprintf('Density ratio variation with x-Nonconservative'))
xlabel('time step')
ylabel('densityratio')
axis([0 1400 0.5 1])
grid on
subplot(3,1,2)
plot(N,tempratio16);
title(sprintf('Temperature ratio variation with x-Nonconservative'))
xlabel('time step')
ylabel('tempratio')
axis([0 1400 0.5 1])
grid on
subplot(3,1,3)
plot(N,pressratio16);
title(sprintf('Pressure ratio variation with x-Nonconservative'))
xlabel('time step')
ylabel('pressureratio')
axis([0 1400 0 1])
grid on
endCONSERVATIVE FORM
function [x,a,v,rho,t,mass,mach]=conservative(n,x,dx,a,gamma,N,nt)
%inputs
x=linspace(0,3,n);
dx=x(2)-x(1);
a=1+2.2*(x-1.5).^2;
gamma=1.4;
for i=1:length(x)
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);
pr=rho.*t;
u1=rho.*a;
u2=rho.*a.*v;
u3=u1.*((t./(gamma-1))+((gamma/2).*(v.^2)));
c=0.5;
dt=c.*(dx./((t.^0.5)+v));
dt=min(dt);
dt=round(dt,3);
for k=1:nt
%timestep
u1old=u1;
u2old=u2;
u3old=u3;
f1=u2;
term=(((gamma-1)/gamma)*(u3-((gamma/2)*((u2.^2)./u1))));
f2=((u2.^2)./u1)+term;
f3=((gamma.*u2.*u3)./u1)-(gamma*(gamma-1)*(u2.^3))./(2*u1.^2);
for i=2:n-1
%predicted value
df1dx(i)=(f1(i+1)-f1(i))/dx;
df2dx(i)=(f2(i+1)-f2(i))/dx;
df3dx(i)=(f3(i+1)-f3(i))/dx;
dlogadx(i)=(log(a(i+1))-log(a(i)))/dx;
dadx1(i)=(a(i+1)-a(i))/dx;
j2(i)=(1/gamma)*(rho(i)*t(i))*dadx1(i);
du1dt_p(i)=-df1dx(i);
du2dt_p(i)=-df2dx(i)+j2(i);
du3dt_p(i)=-df3dx(i);
%updated value
u1(i)=u1(i)+du1dt_p(i)*dt;
u2(i)=u2(i)+du2dt_p(i)*dt;
u3(i)=u3(i)+du3dt_p(i)*dt;
end
rho=(u1./a);
v=(u2./u1);
t=(gamma-1)*((u3./u1)-((gamma*v.^2)/2));
f1=u2;
term1=(((gamma-1)/gamma)*(u3-(gamma*u2.^2)./(2*u1)));
f2=((u2.^2)./u1)+term1;
f3=((gamma.*u2.*u3)./u1)-(gamma*(gamma-1)*(u2.^3))./(2*u1.^2);
for i=2:n-1
%corrected value
df1dxc(i)=(f1(i)-f1(i-1))/dx;
df2dxc(i)=(f2(i)-f2(i-1))/dx;
df3dxc(i)=(f3(i)-f3(i-1))/dx;
dlogadx(i)=(log(a(i))-log(a(i-1)))/dx;
dadx1(i)=(a(i)-a(i-1))/dx;
j3(i)=(1/gamma)*(rho(i)*t(i))*dadx1(i);
du1dt_c(i)=-df1dxc(i);
du2dt_c(i)=-df2dxc(i)+j3(i);
du3dt_c(i)=-df3dxc(i);
end
%average value
du1dt=0.5*(du1dt_p+du1dt_c);
du2dt=0.5*(du2dt_p+du2dt_c);
du3dt=0.5*(du3dt_p+du3dt_c);
%solution value
for i=2:n-1
u1(i)=u1old(i)+du1dt(i)*dt;
u2(i)=u2old(i)+du2dt(i)*dt;
u3(i)=u3old(i)+du3dt(i)*dt;
end
%boundary condition
%inlet condition
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 condition
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));
%deducing the values
rho=(u1./a);
v=(u2./u1);
t=(gamma-1).*((u3./u1)-((gamma*v.^2)/2));
p=rho.*t;
mass=rho.*v.*a;
mach=v./((t).^0.5);
figure(5);
subplot(3,1,1)
plot(x,v);
axis([0 3 0 2])
title(sprintf('Steady-state distribution of primitive variables of Conservative %d',k))
xlabel('x')
ylabel('v')
grid on
subplot(3,1,2)
plot(x,t);
axis([0 3 0 2])
xlabel('x')
ylabel('t')
grid on
subplot(3,1,3)
plot(x,rho);
axis([0 3 0 2])
xlabel('x')
ylabel('rho')
grid on
pause(0.03)
figure(6)
plot(x,mach);
xlabel('x')
ylabel('mach')
axis([0 3 0 5])
grid on
title(sprintf('Steady-state distribution of primitive variables of Conservative %d',k))
%pause(0.03)
figure(7)
if(k==50)
plot(x,mass,'o');
xlabel('x')
ylabel('massflow rate')
grid on
hold on
elseif(k==100)
plot(x,mass,'r');
elseif(k==200)
plot(x,mass,'g');
elseif(k==300)
plot(x,mass,'b');
elseif(k==400)
plot(x,mass,'c');
elseif(k==500)
plot(x,mass,'m');
elseif(k==600)
plot(x,mass,'y');
elseif(k==700)
plot(x,mass,'k');
elseif(k==1400)
plot(x,mass,'*');
title(sprintf('Variation of mass flow rate with Conservative'))
legend('50','100','200','300','400','500','600','700','1400')
end
T1= table(x,a,v,rho,t);
ax=a(16);
rhoo=rho(1);
ao=(t(1))^0.5;
to=t(1);
po=p(1);
arearatio=a/ax;
densityratio=rho/rhoo;
velocitybyao=v/ao;
tempratio=t/to;
pressratio=p/po;
densityratio16(k)=densityratio(16);
tempratio16(k)=tempratio(16);
pressratio16(k)=pressratio(16);
end
densityratio16=densityratio16';
tempratio16=tempratio16';
pressratio16=pressratio16';
figure(8);
subplot(3,1,1)
plot(N,densityratio16);
title(sprintf('Density ratio variation with conservative'))
xlabel('time step')
ylabel('densityratio')
axis([0 1400 0.5 1])
grid on
subplot(3,1,2)
plot(N,tempratio16);
title(sprintf('Temperature ratio variation with conservative'))
xlabel('time step')
ylabel('tempratio')
axis([0 1400 0.5 1])
grid on
subplot(3,1,3)
plot(N,pressratio16);
title(sprintf('Pressure ratio variation with conservative'))
xlabel('time step')
ylabel('pressureratio')
axis([0 1400 0 1])
grid on
end
RESULTS
N=31
The above plots provide the variation of velocity ,Temperature and Density along the length of the domain of the nozzle for both non conservative method and conservative method.Both The methods provide nearly exact results.
The velocity increases as the length of the nozzle increases.
Temperature and Density decreases along the length of the nozzle.
The Mach number increases along the length of the nozzle.
The mass flow rate at differt time steps are plotted for both Non conservative and Conservative method.The massflow rate is not linear untill 700 iterations after that it stabilizes and provides linear results.
The above plots depict the variation of pressure ratio,Temperature ratio and density ratio with respect to increase in time step.It is found that in Non conservative method the values of time derivative are very large upto 200 time step and then it becomes smaller at about 400 iteration and then stabilizes.
In conservative the values of time derivatives are small initially itself and also becomes stable at 400 time step
The conservative form is slightly stabler and more accurate than non conservative form while conserving mass flow of the solution as the difference between the results are not that big.
N=61 Grid dependence Test
The timestep is calculated again as dx changes it is approximately equal to twice the previous value.
So While comparing a value of primitive variable for two grid size.Ex- At a particular position in the nozzle N=31 at 500 time step value should be same as value at N=61 ata1000 time step.
The plots obtained at N=31 are nearly similar to the plots obtained at N=61.There
Inference :
1.Complexity in solving conservative form is higher compared to Non conservative form.The primitive variables need to be found from flux variables in conservative form.
2.The conservative form does a better job in conserving mass than non conservative form
3.The Nonconservative form does a better job in finding the values of primitive variables as it has higher quality index and smaller residue.
4.As grid size increases time step increases and time required for solution also increases.