Week 7 - Simulation of a 1D Super-sonic nozzle flow simulation using Macormack Method
To simulate the isentropic flow through a quasi 1D subsonic-supersonic nozzle:
subsonic-supersonic nozzle to be solved using two form. They are conservation and non-conservation forms.
Conservation form equation:
∂(ρ’A’)/ ∂t’ + ∂(ρ’A’v’)/ ∂x’ = 0
∂(ρ’A’v’)/ ∂t’ + ∂(ρ’A’v2’)/ ∂x’ = -A’*∂ρ’/ ∂x’
∂(A’ρ’v2’/2)/ ∂t’ + ∂(A’ρ’v3’/2)/ ∂x’ = -A’v’*∂ρ’/ ∂x’
Non-conservation form equation:
continuety equation: ∂(ρ’A’)/ ∂t’ + v’*∂(ρ’A’)/ ∂x’ = -ρ’A’*∂v’/∂x’
momentum equation: ρ’∂v’/ ∂t’ + p’v’∂v’/ ∂x’ = -∂p’/ ∂x’
energy equation: ρ’∂e/ ∂t’ + ρ’v’∂e/ ∂x’ = p’∂v’/ ∂x’ – ρ’v’∂(lnA’)/ ∂x’ e = CvT
Governing equation:
continuety equation: ∂ρ’/dt’ = - ρ’*∂v’/∂x’ – ρ’*v’*∂(lnA’)/ ∂x’ – v’*∂ρ’/∂x
momentum equation: ∂v’/∂t’ = v’*∂v’/∂x’ - (1/γ)*( ∂T’/∂x’+(T’/ρ’)* ∂ρ’/∂x’)
energy equation: ∂T’/∂t’ = -v’*∂T’/∂x’ –(γ-1)*T*(∂v’/∂x’ + v’*∂ (lnA’)/ ∂x’)
The equations are changed to dimensionsless units;
x’ = x/L
v’ = v/a0
t’ = t/(L/a0)
A’ = A/A*
The boundary conditions and the input value like area, distance, etc which are in non-dimension units are converted using the simple equations mentioned in the john andersons book of CFD, in 7th chapter.
MacCormack's method:
the MacCormack's method is used to solve the governing equation. this method consist of two steps. first step is predictive step where the governing equation is solved using forward diferencing. THen the second step is corrected step where backward differencing is done. Then the soution is obtained from the avarage of both predictive and corrected step. this average is multipiled with the timestep so that this resulting value is added to the initial value to get results.
main Program:
clear all
close all
clc
n=31;
x=linspace(0,3,n);
dx=x(2)-x(1);
gamma=1.4;
C=0.5;
error1=1e9;
error2=1e9;
tol=1e-5;
a=1+2.2*(x-1.5).^2;
n_throat=(n-1)/2;
tic
[rho , v , t , p , Mach ,m ,rho_throat , t_throat,p_throat,Mach_throat,k]=NonConservative_Form_Fun(n,x,dx,gamma,a,n_throat,C,error1,tol);
time_nc = toc;
fprintf(' Minimum Iteration for Non-Conservative Form: %fn',k);
fprintf('Computation time for Non-Conservative Form: %f sn',time_nc);
Transient_Throat=figure(2);
subplot(4,1,1);
plot(rho_throat,'m');
title('Variations of the properties at the nozzle throat for Non-Conservative Form');
legend('Density ratio');
ylabel('rho/rho_{o}');
axis([0 k 0.5 1]);
grid minor;
subplot(4,1,2);
plot(t_throat,'r');
legend('Temperature');
ylabel('T/T_{o}');
axis([0 k 0.6 1]);
grid minor;
subplot(4,1,3);
plot(p_throat,'g');
legend('Pressure');
ylabel('p/p_{o}');
axis([0 k 0.4 1]);
grid minor;
subplot(4,1,4);
plot(Mach_throat,'b');
legend('Mach Number');
ylabel('M');
axis([0 k 0.8 1]);
grid minor;
xlabel('Number of iterations');
figure(3)
subplot(4,1,1);
plot(x,rho,'m');
title('Profiles of the properties at the Nozzle for Non-Conservative Form');
ylabel('rho/rho_{o}');
legend('Density ratio','location','southwest');
grid minor;
subplot(4,1,2);
plot(x,t,'r');
ylabel('T/T_{o}');
legend('Temperature ratio','location','southwest');
grid minor;
subplot(4,1,3);
plot(x,p,'g');
ylabel('p/p_{o}');
legend('Pressure ratio','location','southwest');
grid minor;
subplot(4,1,4);
plot(x,Mach,'b');
ylabel('M');
legend('Mach Number','location','northwest');
grid minor;
tic
[rho_c , v_c , t_c , p_c , Mach_c , m_c , rho_throat_c , t_throat_c , p_throat_c , Mach_throat_c ,k2]=Conservative_Form_Fun(n,x,dx,gamma,a,n_throat,C,error1,tol);
time_c = toc;
fprintf('Minimum Iterations for Convergence for Conservative Form: %fn',k2);
fprintf('Computation time for Conservative Form: %f sn',time_c);
Transient_Throat_con=figure(5);
subplot(4,1,1);
plot(rho_throat_c,'m');
title('Timewise variations of the properties at the nozzle throat for Conservative Form');
legend('Density ratio');
ylabel('rho/rho_{o}');
axis([0 k2 0.5 1]);
grid minor;
subplot(4,1,2);
plot(t_throat_c,'r');
legend('Temperature');
ylabel('T/T_{o}');
axis([0 k2 0.6 1]);
grid minor;
subplot(4,1,3);
plot(p_throat_c,'g');
legend('Pressure');
ylabel('p/p_{o}');
axis([0 k2 0.4 1]);
grid minor;
subplot(4,1,4);
plot(Mach_throat_c,'b');
legend('Mach Number');
ylabel('M');
axis([0 k2 0.8 1]);
grid minor;
xlabel('Number of iterations');
figure(6)
subplot(4,1,1);
plot(x,rho_c,'m');
title('Profiles of the properties at the Nozzle for Conservative Form');
ylabel('rho/rho_{o}');
legend('Density ratio','location','southwest');
grid minor;
subplot(4,1,2);
plot(x,t_c,'r');
ylabel('T/T_{o}');
legend('Temperature ratio','location','southwest');
grid minor;
subplot(4,1,3);
plot(x,p_c,'g');
ylabel('p/p_{o}');
legend('Pressure ratio','location','southwest');
grid minor;
subplot(4,1,4);
plot(x,Mach_c,'b');
ylabel('M');
legend('Mach Number','location','northwest');
grid minor;
figure(7)
plot(x,m)
hold on
plot(x,m_c)
title('Steady Mass-Flow Rate Profile')
legend('Non-Conservative','Conservative')
grid minor
function for conservative form:
function [rho_c , v_c , t_c , p_c , Mach_c , m_c , rho_throat_c , t_throat_c , p_throat_c , Mach_throat_c , k2]=Conservative_Form_Fun(n,x,dx,gamma,a,n_throat,C,error2,tol)
for i = 1:length(x)
if x(i)>=0 && x(i)<0.5
rho_c(i)=1;
t_c(i)=1;
else if x(i)>=0.5 && x(i)< 1.5
rho_c(i)=1.0-0.366*(x(i)-0.5);
t_c(i)=1.0-0.167*(x(i)-0.5);
else if x(i)>=1.5 && x(i)<=3.0
rho_c(i)=0.634-0.3879*(x(i)-1.5);
t_c(i)=0.833-0.3507*(x(i)-1.5);
end
end
end
end
v_c =0.59./(rho_c.*a);
p_c=rho_c.*t_c;
U1=rho_c.*a;
U2=rho_c.*v_c.*a;
U3=rho_c.*a.*((t_c./(gamma - 1)) + (gamma/2)*v_c.^2);
k2=0;
while error2>tol
k2=k2+1;
massflow_old_con=rho_c.*v_c.*a;
U1_old=U1;
U2_old=U2;
U3_old=U3;
F1= U2;
F2=(U2.^2./U1)+((gamma-1)/gamma)*(U3-(gamma*U2.^2)./(2*U1));
F3=(gamma*U2.*U3./U1)-(gamma*(gamma-1)/2)*(U2.^3./U1.^2);
a_con=sqrt(t_c);
TimeStep=C*dx./(a_con+v_c);
dt=min(TimeStep);
for j=2:n-1
J2(j)=(1/gamma).*rho_c(j).*t_c(j).*((a(j+1)-a(j))/dx);
dU1dt_p(j)=-(F1(j+1)-F1(j))/dx;
dU2dt_p(j)=-(F2(j+1)-F2(j))/dx+J2(j);
dU3dt_p(j)=-(F3(j+1)-F3(j))/dx;
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_c=U1./a;
v_c=U2./U1;
t_c=(gamma-1)*(U3./U1-(gamma/2)*v_c.^2);
p_c=rho_c.*t_c;
F1=U2;
F2=(U2.^2./U1)+((gamma-1)/gamma)*(U3-(gamma*U2.^2)./(2*U1));
F3=(gamma*U2.*U3./U1)-(gamma*(gamma-1)/2)*(U2.^3./U1.^2);
for j=2:n-1
J2(j)=(1/gamma)*rho_c(j).*t_c(j)*((a(j)-a(j-1))/dx);
dU1dt_c(j)=-(F1(j)-F1(j-1))/dx;
dU2dt_c(j)=-(F2(j)-F2(j-1))/dx+J2(j);
dU3dt_c(j)=-(F3(j)-F3(j-1))/dx;
end
dU1dt_avg=0.5*(dU1dt_p+dU1dt_c);
dU2dt_avg=0.5*(dU2dt_p+dU2dt_c);
dU3dt_avg=0.5*(dU3dt_p+dU3dt_c);
for i=2:n-1
U1(i)=U1_old(i)+dU1dt_avg(i)*dt;
U2(i)=U2_old(i)+dU2dt_avg(i)*dt;
U3(i)=U3_old(i)+dU3dt_avg(i)*dt;
end
rho_c(1)=1;
t_c(1)=1;
U1(1)=rho_c(1)*a(1);
U2(1)=2*U2(2)-U2(3);
U3(1)=U1(1)*(t_c(1)/(gamma-1)+(gamma/2)*v_c(1).^2);
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);
rho_c=U1./a;
v_c=U2./U1;
t_c=(gamma -1)*(U3./U1-(gamma/2)*v_c.^2);
p_c=rho_c.*t_c;
Mach_c=v_c./sqrt(t_c);
m_c=rho_c.*v_c.*a;
p_c=rho_c.*t_c;
rho_throat_c(k2)=rho_c(n_throat);
t_throat_c(k2)=t_c(n_throat);
p_throat_c(k2)=p_c(n_throat);
Mach_throat_c(k2)=Mach_c(n_throat);
error2=max(abs(m_c-massflow_old_con));
TransientMassFlow_Con=figure(4);
hold on;
if k2 == 100
plot(x,m_c,'g');
else if k2 == 200
plot(x,m_c,'r');
else if k2 == 700
plot(x,m_c,'-.b','linewidth',1.25);
end
end
end
end
title('Non-dimensional Mass Flow Rates of Conservative form');
xlabel('x/L');
ylabel('rhoVA/rho_{o}V_{o}A_{o}');
legend('100Deltat','200Deltat','700Deltat');
axis([0 3 0.5 0.8]);
grid minor;
end
Function for non-conservative form:
function [rho , v , t , p , Mach , m , rho_throat , t_throat , p_throat , Mach_throat ,k]=NonConservative_Form_Fun (n,x,dx,gamma,a,n_throat,C,error,tol)
rho=1-0.3146*x;
t=1-0.2314*x;
v=(0.1+1.09*x).*t.^0.5;
k=0;
while error>tol
k=k+1;
rho_old=rho;
v_old=v;
t_old=t;
a_old=sqrt(t_old);
m_old=rho_old.*v_old.*a;
tt=C*dx./(a_old+v_old);
dt=min(tt);
for j=2:n-1
dvdx=(v(j+1)-v(j))/dx;
dlogadx=(log(a(j+1))-log(a(j)))/dx;
drhodx=(rho(j+1)-rho(j))/dx;
dtdx=(t(j+1)-t(j))/dx;
drho_dt_p(j)= -rho(j)*dvdx -rho(j)*v(j)*dlogadx -v(j)*drhodx ;
dvdt_p(j)=-v(j)*dvdx -(1/gamma)*(dtdx + t(j)/rho(j)*drhodx);
dtdt_p(j)=-v(j)*dtdx -(gamma-1)*t(j)*(dvdx + v(j)*dlogadx);
rho(j)=rho(j)+drho_dt_p(j)*dt;
v(j)=v(j)+dvdt_p(j)*dt;
t(j)=t(j)+dtdt_p(j)*dt;
end
for j=2:n-1
dvdx=(v(j)-v(j-1))/dx;
dlogadx=(log(a(j))-log(a(j-1)))/dx;
drhodx=(rho(j)-rho(j-1))/dx;
dtdx=(t(j)-t(j-1))/dx;
drho_dt_c(j)= -rho(j)*dvdx -rho(j)*v(j)*dlogadx -v(j)*drhodx ;
dvdt_c(j)=-v(j)*dvdx -(1/gamma)*(dtdx + t(j)/rho(j)*drhodx);
dtdt_c(j)=-v(j)*dtdx -(gamma-1)*t(j)*(dvdx + v(j)*dlogadx);
end
drhodt_avg=0.5*(drho_dt_p+drho_dt_c);
dvdt_avg=0.5*(dvdt_p+dvdt_c);
dTdt_avg=0.5*(dtdt_p+dtdt_c);
for i=2:n-1
rho(i)=rho_old(i)+drhodt_avg(i)*dt;
v(i)=v_old(i)+dvdt_avg(i)*dt;
t(i)=t_old(i)+dTdt_avg(i)*dt;
end
rho(1)=1;
v(1)=2*v(2)-v(3);
t(1)=1;
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);
Mach=v./sqrt(t);
m=rho.*v.*a;
p=rho.*t;
rho_throat(k)=rho(n_throat);
t_throat(k)=t(n_throat);
p_throat(k)=p(n_throat);
Mach_throat(k)=Mach(n_throat);
error=max(abs(m-m_old));
TransientMassFlow=figure(1);
hold on;
if k== 1
plot(x,m,'--k');
else if k == 50
plot(x,m,'r');
else if k == 100
plot(x,m,'g');
else if k == 150
plot(x,m,'c');
else if k == 200
plot(x,m,'b');
else if k == 700
plot(x,m,'-.b','linewidth',1.25);
plot(x(end),m(end),'k*');
end
end
end
end
end
end
end
title('Non-dimensional Mass Flow Rates in Non-Conservative form');
xlabel('x/L');
ylabel('rhoVA/rho_{o}V_{o}A_{o}');
legend('0Deltat','50Deltat','100Deltat','150Deltat','200Deltat','700Deltat');
axis([0 3 0.3 2]);
grid minor;
end
result graphs:
conservative form:
The graph shows that when del(t)= 700, in the conservative form it reaches steady state with flow for the ratio of dimensionless dnesity, velocity and area with distance.
THis graph show the increase in mach number at exit of nozzle, where the expansion of gasses takes place decreasing the pressure, density and temperature.
In this graph, it shows the timewise variation in conservative form at the postion of throath area of the nozzle, where the mach number, density, pressure and temperature reaches a steady flow through a ceratain iteration
non-conservative form:
In this nonconservative form, the mass flow reaches a steady state del(t)=700 as timestep.
since this graph show the exit for of nozzle for non-consevative form, result are as aspected. the mach number eventually with the expansion of masses that flow throught nozzle , decreasing the pressure, temperature and density.
in this timewise graph, it states the number of iteration that took to reach steady state for temperature, pressure and density in its dimensionless state. the area is constant so that in this time changing iteration gives a time marching result at thhat area.
steady mass-flow graph:
On comparing, the conservative and non-conservative for with the mass flow graph and time taken to slove; the conservatives might fluctuate a lot to reach the answer but both form produces the necessary result as seen in the graph. But the time takeen and iteration shows the effeicient one to solve.
conservative form:
number of iteration: 951 time taken: 3.36s
non-conservative form:
number of iteration: 767 time taken: 3.03s
Thus the Non-conservative form takes less time and less numbe of iteration the derive the solution from the governing equation while comparing it with conservative form. this show that non-conservative form is very efficient than conservative form.