Simulation of a 1D Super-sonic nozzle flow simulation using Macormack Method
OBJECTIVE
To simulate the isentropic flow through a quasi 1D subsonic-supersonic nozzle for both conservative and non-conservative forms of governing equations using the MacCormack's technique.
THEORY
Consider a steady state isentropic flow through a convergent-divergent nozzle as shown in figure 1.
Figure 1: Schematic representation for subsonic-supersonic nozzle flow
The inlet is connected to a large reservoir, with initial parameters as pressure (Po) and temperature (To). Since the area of the reservoir is assumed to be ∞, initial velocity Vo is very small, almost 0. Also, it is a subsonic region, hence, mach number M is less than 1.The flow undergoes isentropc expansion, and reaches supersonic condition (M>1) at the exit with Pe,Te,Ve and Ae as the parameters.At the throat section, sonic flow is assumed (M=1), and hence, the flow velocity is equal to the speed of sound at this region, denoted by a∗, area A∗, pressure P∗ and temperature T∗.
There are mainly two assumptions considered in this case:
- At a given cross-sectional area A, the flow properties are uniform.
- Even though the area changes with distance x and the flow properties are 2 dimensional, we assume 1-directional flow properites, thus called as quasi-one-dimensional flow.
Figure 2: Control volume for deriving partial differential equations for unsteady quasi-1D-flow
Assume a small control volume with parameters as shown in figure 2.
The governing equations that define the flow are continuity, momentum and energy equations.
The NON-CONSERVATIVE FORM of governing equations, which are non-dimensional are
- continuity equation
- momentum equation
- energy equation
where γ is the ratio of specefic heats Cp/Cv, and is 1.4 for standard air conditions.
Initial conditions:
The CONSERVATIVE FORM of governing equations, which are non-dimensional are
- continuity equation
- momentum equation
- energy equation
These equations can be written as
where U1,U2 and U3 are called as the solution vectors; F1,F2 and F3 are the flux vectors; and J is the source term.
Inorder to avoid difficulties during computation and calculations, the flux vectors are represented in terms of solution vectors, i.e
F1=F1(U1,U2,U3)
F2=F2(U1,U2,U3)
F3=F3(U1,U2,U3)
J2=J2(U1,U2,U3)
where,
Initial conditions:
- Calculation of Time Step-Adaptive Time Step Method
These equations are solved by the MacCormack's technique.
Assume f(x) as the function to be solved. Solving this using MacCormack's technique involves
- Initial Profile : Here, certain initial assumptions are considered for the flow field variables.
- Predicted step: Forward differencing scheme is applied on the spatial terms, and obtaining the predicted derivatives.
With these predicted time derivatives, we update the initial conditions.
- Corrected step :Here, backward differencing is employed for the spatial terms.
It is of immense importance to use the updated values of initial conditions from the predicted step.
- Average step : Here, we compute the average values of time derivatives from predicted and corrected time steps.
- Final Solution Update:
CODE
MAIN PROGRAM
clear all
close all
clc
%creating x mesh
n=31;
x=linspace(0,3,n);
dx=x(2)-x(1);
CFL=0.5;
throat=(n+1)/2;
%area
a=1+2.2*(x-1.5).^2;
%initial profile
[m_dot_throat_nc,M_throat_nc,T_throat_nc,P_throat_nc,rho_throat_nc,t,M_nc,P_nc,rho,v,m_dot_nc]=non_conservative(n,x,dx,throat,a,CFL);
[m_dot_throat,M_throat,T_throat,P_throat,rho_throat,t,M,P,rho,v,m_dot]=conservative(n,x,dx,throat,a,CFL)
%COMPARISON GRAPHS
figure(1)
hold on
plot(x,m_dot,'-*')
plot(x,m_dot_nc,'-*')
ylim([0.575 0.62])
set(gca,'YTick',(0.575:0.005:0.62))
xlabel('x')
ylabel('mass flow rate')
legend('conservative','non-conservative')
title('mass flow rate- conservative vs non-conservative')
grid on
figure(2)
hold on
plot(M_throat,'linewidth',0.95)
plot(M_throat_nc,'linewidth',0.95)
xlabel('timesteps')
ylabel('mach number')
legend('conservative','non-conservative')
title('mach number variation at throat for conservative and non-conservative')
figure(3)
hold on
plot(T_throat,'linewidth',0.95)
plot(T_throat_nc,'linewidth',0.95)
xlabel('timesteps')
ylabel('Temperature')
legend('conservative','non-conservative')
title('Temperature variation at throat for conservative and non-conservative')
figure(4)
hold on
plot(P_throat,'linewidth',0.95)
plot(P_throat_nc,'linewidth',0.95)
xlabel('timesteps')
ylabel('Pressure')
legend('conservative','non-conservative')
title('Pressure variation at throat for conservative and non-conservative')
figure(5)
hold on
plot(rho_throat,'linewidth',0.95)
plot(rho_throat_nc,'linewidth',0.95)
xlabel('timesteps')
ylabel('Density')
legend('conservative','non-conservative')
title('Density variation at throat for conservative and non-conservative')
%Non conservative form
figure(6)
subplot(4,1,1)
plot(x,t,'-*')
xlabel('x-distance')
ylabel('temperature ratio')
title('Flow parameters -non-conservative form')
subplot(4,1,2)
plot(x,M_nc,'-*')
xlabel('x-distance')
ylabel('mach number')
subplot(4,1,3)
plot(x,P_nc,'-*')
xlabel('x-distance')
ylabel('pressure ratio')
subplot(4,1,4)
plot(x,rho,'-*')
xlabel('x-distance')
ylabel('density ratio')
figure(7)
hold on
plot(x,t,'linewidth',0.95)
plot(x,rho,'linewidth',0.95,'color','r')
plot(x,v,'linewidth',0.95,'color','g')
plot(x,m_dot_nc,'linewidth',0.95,'color','b')
plot(x,M_nc,'linewidth',0.95,'color','c')
plot(x,P_nc,'linewidth',0.95,'color','y')
xlabel('x-distance')
ylabel('flow parameters')
legend('Temperature','density','velocity','mass flow rate','mach number','pressure')
title('Flow parameters -non-conservative form')
ylim([0 2])
set(gca,'YTick',(0:0.2:2))
figure(8)
plot(T_throat_nc,'linewidth',0.95)
xlabel('timesteps')
ylabel('temperature ratio at throat')
title('Temperature vs time-step -non-conservative form')
figure(9)
plot(P_throat_nc,'linewidth',0.95)
xlabel('timesteps')
ylabel('pressure ratio at throat')
title('pressure vs time-step -non-conservative form')
figure(10)
plot(rho_throat_nc,'linewidth',0.95)
xlabel('timesteps')
ylabel('density ratio at throat')
title('density vs time-step -non-conservative form')
figure(11)
plot(M_throat_nc,'linewidth',0.95)
xlabel('timesteps')
ylabel('mach number at throat')
title('mach number vs time-step -non-conservative form')
figure(12)
plot(m_dot_throat_nc,'linewidth',0.95)
xlabel('timesteps')
ylabel('mass flow rate at throat')
title('mass flow rate vs time-step -non-conservative form')
%concervative graphs
figure(13)
hold on
plot(x,t,'linewidth',0.95)
plot(x,rho,'linewidth',0.95,'color','r')
plot(x,v,'linewidth',0.95,'color','g')
plot(x,m_dot,'linewidth',0.95,'color','b')
plot(x,M,'linewidth',0.95,'color','c')
plot(x,P,'linewidth',0.95,'color','y')
xlabel('x-distance')
ylabel('flow parameters')
legend('Temperature','density','velocity','mass flow rate','mach number','pressure')
title('Flow parameters -conservative form')
ylim([0 2])
set(gca,'YTick',(0:0.2:2))
figure(14)
plot(T_throat,'linewidth',0.95)
xlabel('timesteps')
ylabel('temperature ratio at throat')
title('Temperature vs time-step -conservative form')
figure(15)
plot(P_throat,'linewidth',0.95)
xlabel('timesteps')
ylabel('pressure ratio at throat')
title('Pressure vs time-step -conservative form')
figure(16)
plot(rho_throat,'linewidth',0.95)
xlabel('timesteps')
ylabel('density ratio at throat')
title('Density vs time-step -conservative form')
figure(17)
plot(M_throat,'linewidth',0.95)
xlabel('timesteps')
ylabel('mach number at throat')
title('Mach number vs time-step -conservative form')
figure(18)
plot(m_dot_throat,'linewidth',0.95)
xlabel('timesteps')
ylabel('mass flow rate at throat')
title('mass flow rate vs time-step -conservative form')
FUNCTION
function [m_dot_throat_nc,M_throat_nc,T_throat_nc,P_throat_nc,rho_throat_nc,t,M_nc,P_nc,rho,v,m_dot_nc]=non_conservative(n,x,dx,throat,a,CFL);
clear all
close all
clc
%creating x mesh
n=31;
x=linspace(0,3,n);
dx=x(2)-x(1);
throat=(n+1)/2;
a=1+2.2*(x-1.5).^2; %area
CFL=0.5;
rho=1-0.3146*x; %pressure
t=1-0.2314*x; %temperature
v=(0.1+1.09*x).*t.^0.5; %velocity
gamma=1.4;
%calculating the time step based on cfl number
for d=1:n
A=sqrt(t(d)); %non-dimensional eqn for speed of sound
dt_i(d)=CFL*(dx/(v(d)+A));
end
dt=min(dt_i);
nt=1400; %number of iterations
%NON-CONSERVATIVE FORM
for k=1:nt %outer time loop
%creating a copy of the initial conditions
rho_old=rho;
v_old=v;
T_old=t;
%step 1: predicted step-forward differencing for spatial terms
for i=2:n-1
%defining spatial terms
dv_dx=(v(i+1)-v(i))/dx;
drho_dx=(rho(i+1)-rho(i))/dx;
dlog_A_dx=(log(a(i+1))-log(a(i)))/dx;
dt_dx=(t(i+1)-t(i))/dx;
%continuity equation
drho_dt_p(i)=-rho(i)*dv_dx-rho(i)*v(i)*dlog_A_dx-v(i)*drho_dx;
%momentum equation
dv_dt_p(i)=-v(i)*dv_dx-(1/gamma)*(dt_dx+(t(i)/rho(i))*drho_dx);
%energy equation
dT_dt_p(i)=-v(i)*dt_dx-(gamma-1)*t(i)*(dv_dx+v(i)*dlog_A_dx);
%first solution update-updating the initial conditions
rho(i)=rho_old(i)+drho_dt_p(i)*dt;
v(i)=v_old(i)+dv_dt_p(i)*dt;
t(i)=T_old(i)+dT_dt_p(i)*dt;
end %predicted step end
%step 2: corrected step-backward differencing for spatial terms
for i=2:n-1
%defining spatial terms
dv_dx=(v(i)-v(i-1))/dx;
drho_dx=(rho(i)-rho(i-1))/dx;
dlog_A_dx=(log(a(i))-log(a(i-1)))/dx;
dt_dx=(t(i)-t(i-1))/dx;
%continuity equation
drho_dt_c(i)=-rho(i)*dv_dx-rho(i)*v(i)*dlog_A_dx-v(i)*drho_dx;
%momentum equation
dv_dt_c(i)=-v(i)*dv_dx-(1/gamma)*(dt_dx+(t(i)/rho(i))*drho_dx);
%energy equation
dT_dt_c(i)=-v(i)*dt_dx-(gamma-1)*t(i)*(dv_dx+v(i)*dlog_A_dx);
end %corrected step 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 i=2:n-1
rho(i)=rho_old(i)+drho_dt_avg(i)*dt;
v(i)=v_old(i)+dv_dt_avg(i)*dt;
t(i)=T_old(i)+dT_dt_avg(i)*dt;
end
%applying boundary conditions- should be done after final solution
%update
%inlet
v(1)=2*v(2)-v(3);
rho(1)=1;
t(1)=1;
%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);
%finding primitive variables
m_dot_nc= rho.*a.*v; %mass flow rate
M_nc=v./sqrt(t); %mach number
P_nc=rho.*t;
%taking values at throat
v_throat_nc(k)=v(throat);
rho_throat_nc(k)=rho(throat);
T_throat_nc(k)=t(throat);
P_throat_nc(k)=P_nc(throat);
M_throat_nc(k)=M_nc(throat);
m_dot_throat_nc(k)=m_dot_nc(throat);
if k==50
figure(19)
plot(x,m_dot_nc,'linewidth',0.95)
hold on
elseif k==100
figure(19)
plot(x,m_dot_nc,'linewidth',0.95)
hold on
elseif k==200
figure(19)
plot(x,m_dot_nc,'linewidth',0.95)
hold on
elseif k==400
figure(19)
plot(x,m_dot_nc,'linewidth',0.95)
hold on
elseif k==800
figure(19)
plot(x,m_dot_nc,'linewidth',0.95)
hold on
elseif k==1400
figure(19)
plot(x,m_dot_nc,'linewidth',0.95)
xlabel('non-dimensional x')
ylabel('mass flow rate')
legend('t=50','t=100','t=200','t=400','t=800','t=1400')
title('mass flow rate at different time steps-NON-CONSERVATIVE')
hold on
end
end
figure(6)
subplot(4,1,1)
plot(x,t,'-*')
xlabel('x-distance')
ylabel('temperature ratio')
title('Flow parameters -non-conservative form')
subplot(4,1,2)
plot(x,M_nc,'-*')
xlabel('x-distance')
ylabel('mach number')
subplot(4,1,3)
plot(x,P_nc,'-*')
xlabel('x-distance')
ylabel('pressure ratio')
subplot(4,1,4)
plot(x,rho,'-*')
xlabel('x-distance')
ylabel('density ratio')
figure(7)
hold on
plot(x,t,'linewidth',0.95)
plot(x,rho,'linewidth',0.95,'color','r')
plot(x,v,'linewidth',0.95,'color','g')
plot(x,m_dot_nc,'linewidth',0.95,'color','b')
plot(x,M_nc,'linewidth',0.95,'color','c')
plot(x,P_nc,'linewidth',0.95,'color','y')
xlabel('x-distance')
ylabel('flow parameters')
legend('Temperature','density','velocity','mass flow rate','mach number','pressure')
title('Flow parameters -non-conservative form')
ylim([0 2])
set(gca,'YTick',(0:0.2:2))
figure(8)
plot(T_throat_nc,'linewidth',0.95)
xlabel('timesteps')
ylabel('temperature ratio at throat')
title('Temperature vs time-step -non-conservative form')
figure(9)
plot(P_throat_nc,'linewidth',0.95)
xlabel('timesteps')
ylabel('pressure ratio at throat')
title('pressure vs time-step -non-conservative form')
figure(10)
plot(rho_throat_nc,'linewidth',0.95)
xlabel('timesteps')
ylabel('density ratio at throat')
title('density vs time-step -non-conservative form')
figure(11)
plot(M_throat_nc,'linewidth',0.95)
xlabel('timesteps')
ylabel('mach number at throat')
title('mach number vs time-step -non-conservative form')
figure(12)
plot(m_dot_throat_nc,'linewidth',0.95)
xlabel('timesteps')
ylabel('mass flow rate at throat')
title('mass flow rate vs time-step -non-conservative form')
end
function [m_dot_throat,M_throat,T_throat,P_throat,rho_throat,t,M,P,rho,v,m_dot]=conservative(n,x,dx,throat,a,CFL)
n=31;
x=linspace(0,3,n);
dx=x(2)-x(1);
throat=(n+1)/2;
CFL=0.5;
%area
a=1+2.2*(x-1.5).^2;
for i=1:n
if x(i)>=0 && x(i)<=0.5
rho(i)=1;
t(i)=1;
elseif 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);
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.59./(rho.*a);
gamma=1.4;
%adaptive time step method
for d=1:n
A=sqrt(t(d));
dt_i(d)=CFL*(dx/(v(d)+A));
end
dt=min(dt_i);
nt=1400;
%defining the solution vectors
u_1=rho.*a;
u_2=u_1.*v;
u_3=u_1.*((t./(gamma-1))+((gamma/2)*v.^2));
%MACORMACK'S METHOD
for k=1:nt
%predictor step
%initialising the flux terms
term_1=(u_2.^2./u_1);
term_2=((gamma-1)/gamma)*(u_3-(gamma/2)*term_1);
term_3=(gamma*u_2.*u_3)./u_1;
term_4=0.5*gamma*(gamma-1)*u_2.^3./u_1.^2;
F_1=u_2;
F_2=term_1+term_2;
F_3=term_3-term_4;
%creating a copy of the solution vectors
u_1_old=u_1;
u_2_old=u_2;
u_3_old=u_3;
for j=2:n-1
%defining spatial terms-forward differencing in spatial terms
dF_1_dx(j)= (F_1(j+1)-F_1(j))/dx;
dF_2_dx(j)= (F_2(j+1)-F_2(j))/dx;
dF_3_dx(j)= (F_3(j+1)-F_3(j))/dx;
%defining the source term-forward differencing
dA_dx=((a(j+1))-(a(j)))/dx;
J_2=(1/gamma).*rho(j).*t(j).*dA_dx;
%continuity equation
du_1_dt_p(j)=-dF_1_dx(j);
%momentum equation
du_2_dt_p(j)=-dF_2_dx(j)+J_2;
%energy equation
du_3_dt_p(j)=-dF_3_dx(j);
%first solution update
u_1(j)=u_1_old(j)+du_1_dt_p(j)*dt;
u_2(j)=u_2_old(j)+du_2_dt_p(j)*dt;
u_3(j)=u_3_old(j)+du_3_dt_p(j)*dt;
end %end of predicted step
%corrected step
term_1=(u_2.^2./u_1);
term_2=((gamma-1)/gamma)*(u_3-(gamma/2)*term_1);
term_3=(gamma*u_2.*u_3)./u_1;
term_4=0.5*gamma*(gamma-1)*u_2.^3./u_1.^2;
F_1=u_2;
F_2=term_1+term_2;
F_3=term_3-term_4;
rho=u_1./a;
v=u_2./u_1;
t=(gamma-1)*((u_3./u_1)-(0.5*gamma*v.^2));
for j=2:n-1
%defining spatial terms-backward differencing in spatial terms
dF_1_dx(j)= (F_1(j)-F_1(j-1))/dx;
dF_2_dx(j)= (F_2(j)-F_2(j-1))/dx;
dF_3_dx(j)= (F_3(j)-F_3(j-1))/dx;
%defining the source term-backward differencing
dA_dx=((a(j))-(a(j-1)))/dx;
J_2=(1/gamma).*rho(j).*t(j).*dA_dx;
%continuity equation
du_1_dt_c(j)=-dF_1_dx(j);
%momentum equation
du_2_dt_c(j)=-dF_2_dx(j)+J_2;
%energy equation
du_3_dt_c(j)=-dF_3_dx(j);
end %end of corrected step
%computing average
du_1_dt_avg=0.5*(du_1_dt_p+du_1_dt_c);
du_2_dt_avg=0.5*(du_2_dt_p+du_2_dt_c);
du_3_dt_avg=0.5*(du_3_dt_p+du_3_dt_c);
%final solution update
for j=2:n-1
u_1(j)=u_1_old(j)+du_1_dt_avg(j)*dt;
u_2(j)=u_2_old(j)+du_2_dt_avg(j)*dt;
u_3(j)=u_3_old(j)+du_3_dt_avg(j)*dt;
end
%initial and boundary condition
%subsonic inflow-point 1
u_1(1)=rho(1).*a(1);
u_2(1)=2*u_2(2)-u_2(3);
u_3(1)=u_1(1)*((t(1)/(gamma-1)+0.5*gamma*v(1).^2));
%supersonic outflow-point n
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);
%calculating the primitive variables from solution vectors
rho=u_1./a;
v=u_2./u_1;
t=(gamma-1)*((u_3./u_1)-(0.5*gamma*v.^2));
P=rho.*t;
M=v./sqrt(t);%mach number
m_dot= rho.*a.*v; %mass flow rate
%taking values at throat
v_throat(k)=v(throat);
rho_throat(k)=rho(throat);
T_throat(k)=t(throat);
P_throat(k)=P(throat);
M_throat(k)=M(throat);
m_dot_throat(k)=m_dot(throat);
if k==50
figure(20)
plot(x,m_dot,'linewidth',0.95)
hold on
elseif k==100
figure(20)
plot(x,m_dot,'linewidth',0.95)
hold on
elseif k==200
figure(20)
plot(x,m_dot,'linewidth',0.95)
hold on
elseif k==400
figure(20)
plot(x,m_dot,'linewidth',0.95)
hold on
elseif k==800
figure(20)
plot(x,m_dot,'linewidth',0.95)
hold on
elseif k==1400
figure(20)
plot(x,m_dot,'linewidth',0.95)
xlabel('non-dimensional x')
ylabel('mass flow rate')
legend('t=50','t=100','t=200','t=400','t=800','t=1400')
title('mass flow rate at different time steps-CONSERVATIVE')
hold on
end
end
PLOTS
The flow parameters (non-dimensional) at the throat are compare
CFL Number Dependency Test
CONCLUSION
- The flow parameters in non-conservative form become constant after the 600th time step, whereas the same occur in conservative form after the 800th time step
- It can be seen from the CFL number test that,CFL number doesn't have a significant impact on conservative form, whereas, in non-conservative form, the parameter values increases with CFL number, howerver, the solution blows-off when it hits 1.2