Week 7 - Simulation of a 1D Super-sonic nozzle flow simulation using Macormack Method
Solution of a 1D Super sonic nozzle flow simulation using MacCormack Method
MacCormack's method is a variant of the Lax -Wendroff approach but is very simpler in the application. It is an explicit finite-difference technique that is second-order accurate in both space and time.It is very much useful in solving hyperbolic equations.
Physical Problem :-
We consider a steady,isentropic flow through a convergent -divergent nozzle .The flow at the inlet to the nozzle comes from a reservoir where the pressure and temperature are denoted by p0 and T0 ,respectively .The cross -sectional area of the reservoir is very large so velcoity can be considered .The flow expands isentropically to supersonic speeds where (Mach number >1) at the nozzle exit .The flow is subsonic at the convergent section of the nozzle ,sonic at the throat (minimum area ) and supersonic at the divergent section .
Schematic Diagram :-
Assumptions :-
We assume that at a given section ,the flow properties are uniform across that section .So ,although the area of the nozzle changes as a function of distance along the nozzle in x direction ,therefore in reality the flow field is two -dimensional (flow properties vary along x-y direction ) but we make the assumption that the flow properties vary only in x direction .Such flow is called as quasi one dimensional flow .
Governing equations :-
The governing equations for 1D steady ,isentropic flow can be expressed as :-
since there is a twisted physics so we apply the integral forms of the conservation equations
and
Control Volume :-
The governing equations for the unsteady state are :-
To reduce the number of variables we non -dimensionalize the variables so that they can be solved analytically ( For equations to be solved ,number of unkowns = number of equations)
The governing equations can be derived in two forms
1) Conservative form
2) Non-conservative form
Matlab Code with Explaination
Main Program:-
close all
clear all
clc
% Solving Quasi 1D Supersonic Nozzle Flow Equations
L=3;%Length of nozzle
n=[31,62,93];% Variation in number of nodes
n_t=1:1400;% Time step
cfl = 0.5;% CFL number
for i = 1:length(n)
x=linspace(0,L,n(i));
[mfr_nc,mfr_throat_nc]=non_conservative(L,n(i),n_t,cfl);
[mfr_c,mfr_throat_c]=conservative(L,n(i),n_t,cfl);
figure(7)
subplot(2,1,1)
plot(n_t,mfr_throat_c,'b',n_t,mfr_throat_nc,'r','linewidth',1)
xlabel('Time steps (n_t)')
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,2)
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(n_t,mfr_throat_c,'linewidth',1)
hold on
xlabel('Time steps (n_t)')
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(n_t,mfr_throat_nc,'linewidth',1)
hold on
xlabel('Time steps (n_t)')
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('Length 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)
plot(x,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')
endCode for Non-Conservative Case :-
function [mfr_nc,mfr_throat_nc]=non_conservative(L,n,n_t,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));
end
dt=min(dt);
for i = 1:length(n_t)
%storing 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;
% Continuty 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);
% Updating primitives
v(j) = v_old(j) + dv_dt_p(j)*dt;
rho(j) = rho_old(j) + drho_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;
drho_dt_c(j) = -rho(j)*dv_dx - rho(j)*v(j)*da_dx - v(j)*drho_dx;
dv_dt_c(j) = -v(j)*dv_dx - (1/g)*(dT_dx + (T(j)/rho(j))*drho_dx);
dT_dt_c(j) = -v(j)*dT_dx - (g-1)*T(j)*(dv_dx + v(j)*da_dx);
end
% average time derivatives
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 corected values
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;
end
% Entry boundary conditions
v(1) = 2*v(2)-v(3);
% Exit boundary conditions
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 update of primitives
p = rho.*T;%Pressure
M= v./T.^0.5;%Mach Number
error(k) = abs(v(k)-v(k-1));
mfr= rho.*a.*v;% Mass Flow Rate
% 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==i
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 the 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(n_t,rho_throat)
hold on
plot(n_t,v_throat,'linewidth',1)
plot(n_t,p_throat,'linewidth',1)
plot(n_t,T_throat,'linewidth',1)
plot(n_t,M_throat,'linewidth',1)
plot(n_t,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 the primitive variables at throat point i=16')
end
mfr_nc=mfr;
mfr_throat_nc=mfr_throat;
endCode for Conservative Case :-
function [mfr_c,mfr_throat_c]=conservative(L,n,n_t,cfl)
x=linspace(0,L,n);
dx=x(2)-x(1);
% Initial conditions rho, a, v, T, M, mfr, 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.167*(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
% area
a= 1+2.2*(x-1.5).^2;
% gamma value
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);
% Mach number
M= v./(T.^0.5);
% Mass flow rate
mfr = rho.*a.*v;
%pressure
p=rho.*T;
% To find the dt value
for z=1:31
u=(T(z))^(1/2);
dt(z) = cfl*dx/(u + v(z));
end
dt=min(dt);
for i = 1:length(n_t)
%storing initial values
u1_old = u1;
u2_old = u2;
u3_old = u3;
for j=2:n-1
% for 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 term
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 Equation
du2_dt_p(j)= -(f2(j+1)-f2(j))/dx + j1;% Momentum Equation
du3_dt_p(j)= -(f3(j+1)-f3(j))/dx;% Energy Equation
% 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-0.5*g*(g-1).*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 Equation
du2_dt_c(j) = -(f2(j)-f2(j-1))/dx + j2;% Momentum Equation
du3_dt_c(j) = -(f3(j)-f3(j-1))/dx;% Energy Equation
end
% average time derivatives
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 corected values
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;
end
% Entry boundary conditions
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);
% Exit boundary conditions
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 update 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 the throat point of the nzzle
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==i
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 the 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(n_t,rho_throat)
hold on
plot(n_t,v_throat,'linewidth',1)
plot(n_t,p_throat,'linewidth',1)
plot(n_t,T_throat,'linewidth',1)
plot(n_t,M_throat,'linewidth',1)
plot(n_t,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 the primitive variables at throat point i=16')
end
mfr_c=mfr;
mfr_throat_c=mfr_throat;
endOutput :-
Non Conservative Form :-
Figure 1
Figure 2
Figure 3
Figure 1 shows the variation of mass flow rate with respect to the length of the nozzle at various time steps. It is found that
mass flow rate at the steady-state is found to be 0.61.
Figure 2 shows the variation of non-dimensional parameters with respect to the length of the nozzle. It is in accordance with the theory.
Figure 3 shows the variation of primitive variables at the throat. Steady State is reached in between time steps of 400 to 600.
Conservative form:-
Figure 4
Figure 5
Figure 6
Figure 4 shows the variation of mass flow rate with respect to the length of the nozzle at different time steps during the time marching process. It is found that the steady state mass flow rate is around 0.58.
Figure 5 shows the steady state distribution of primitive variables inside the nozzle in conservative form and it is found that mass flow rate at the steady state is around 0.52.
Figure 6 shows the variation of non -dimensional parameters with respect to the length of the nozzle at the throat point. Steady state values are reached at the time step of 600.
Non -Conservative v/s Conservative form :-
Figure 7 shows the variation of mass flow rate in both conservative form and non -conservative form at the throat point.
From plot 1, both conservative form and non conservative form reached steady state at the end of the time step. But conservative form reached the steady state early i.e at fewer time steps.
From plot 2, it can be seen that difference between mass flow rate in conservative form and in non conservative form is not much. But the mass flow rate in conservative form has achieved a steady state value whereas in the non conservative form it is slightly varying i.e (0.582 +_ small error ).
Thus it can be concluded that both in non conservative and in conservative form values are not changing much and the stability is achieved more quickly in the case of conservative form.
Grid Independency test for both Conservative and Non -Conservative Form :-
Figure 8
Figure 9
Figure 8 and Figure 9 shows grid independency test in both conservative form and in non -conservative form.
From figure 8 , plot 1 shows the variation of mass flow rate at the different number of nodes for grid independency test at the throat in conservative form and plot 2 shows the same in non conservative form.
It can be concluded that the mass flow rate at the end of the time step is same and has achieved stability at the throat point.
Figure 9 shows the variation of mass flow rate along the length of the nozzle at different number of nodes for grid independency test . Plot 1 shows the variation in conservative form and Plot 2 shows the variation in non-conservative form.
It can be concluded that 3 different solutions are obtained for mass flow rate with variation in the number of nodes but the variation in values is not much even though the number of nodes is tripled. So the variation can be considered insignificant.