Week 7 - Simulation of a 1D Supersonic nozzle flow simulation using Macormack Method
Objective: To simulate the isentropic flow through a quasi 1D subsonic-supersonic nozzle using both the conservation and non-conservation forms of the governing equations and sovle them using the MacCormack's technique and comparing their solutions.
Description: We consider the study, isentropic flow through a convergent divergent nozzle. The flow at the inlet to the nozzle comes from reservoir. The cross section area of the reservoir is very large and hence velocity is very small. The flow is locally subsonic the convergent section of the nozzle, Sonic at the throat and Supersonic at the divergent section.
Let us consider a small control volume of the nozzle as shown below.
We assume that at a given section, where the cross section area is A, the flow properties are uniform across the section. Hence, although the area of the nozzle changes as a function of distance along the nozzle, X, and therefore in reality the flow field is two dimensional(the flow varies in the two dimensional x-y space), We make the assumption that the flow properties vary only with X; this is tantamount to asuming uniform flow properties across any given cross-section. Such flow is defined as Quasi one dimensional flow.
The governing continuty, momentum and energy equations in terms of non dimensional variables for this quassi one dimensional staedy and isentropc flow can be expressed respectively as:
MATLAB CODE :
Matlab Code with explaination:
close all
clear all
clc
L=3; %length of nozzle
n=[31 62 93]; %number of grid points
nt=1:1400; %number of time steps
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),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,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(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_throat_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('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')
endNon-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));
end
dt=min(dt);
for i = 1:length(nt)
%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;
drho_dt_p(j) = -rho(j)*dv_dx - rho(j)*v(j)*da_dx - v(j)*drho_dx; %continuty
dv_dt_p(j) = -v(j)*dv_dx - (1/g)*(dT_dx + (T(j)/rho(j))*drho_dx); %momentum
dT_dt_p(j) = -v(j)*dT_dx - (g-1)*T(j)*(dv_dx + v(j)*da_dx); %energy
%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; %continuity
dv_dt_c(j) = -v(j)*dv_dx - (1/g)*(dT_dx + (T(j)/rho(j))*drho_dx); %momentum
dT_dt_c(j) = -v(j)*dT_dx - (g-1)*T(j)*(dv_dx + v(j)*da_dx); %energy
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;
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==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(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 the primitive variables at throat point i=16')
end
mfr_nc=mfr;
mfr_throat_nc=mfr_throat;
endConservative 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, 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
a= 1+2.2*(x-1.5).^2; %area
g = 1.4; %gamma value
% solution vectors u1, u2, u3
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); %mach number
mfr = rho.*a.*v; %mass flow rate
p=rho.*T; %pressure
for z=1:31 %to find dt value
u=(T(z))^(1/2);
dt(z) = cfl*dx/(u + v(z));
end
dt=min(dt);
for i = 1:length(nt)
%storing 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 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
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-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
du2_dt_c(j) = -(f2(j)-f2(j-1))/dx + j2; %momentum
du3_dt_c(j) = -(f3(j)-f3(j-1))/dx; %energy
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 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(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(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 the primitive variables at throat point i=16')
end
mfr_c=mfr;
mfr_throat_c=mfr_throat;
endOutput:
Non-conservative form:
fig(1)
fig(2)
fig(3)
conservative form:
fig(4)
fig(5)
fig(6)
Non-conservative form vs conservative form:
fig(7)
fig(7) shows comparision of mass flow rate using both conservative and non-conservative forms.
from subplot 1 : Both the forms reached the stable and exact solution by the end of time steps. But the conservative form obatined the stability at lesser time steps.
from subplot 2 : The difference between both the forms is very less in the nozzle. But the mass flow rate in the conservative form is constant where as in non-conservative form it is varying with 0.582±0.003.
From fig(7) we can conclude that the final solution is not much changing when using conservative and non-conservative forms. Although there is a sligth change at the throat area alone in the conservative form, which is not considerably significant. Though the stability is obtained at lesser time steps the computational time is more in case of conservative form than the other.
Grid independency test for both Non-conservative form & Conservative form:
fig(8)
fig(8) and fig(9) shows grid independence test in both the forms.
fig(8) - mass flow rate varying with time step at throat point - subplot(1) - conservative form & subplot(2) - non-conservative form.
At 3 different grid sizes, though the values are different at the starting of time marching process the solution obatined is stable and almost same by the end of time marching.
fig(9) - mass flow rate in the nozzle - subplot(1) - conservative form & subplot(2) - non-conservative form.
The solution obtained using 3 different grid sizes resulting 3 different solutions within the range of 0.589±0.006. i.e. the difference between them is very small which can be neglected/ignored.
From fig(8) and fig(9) we can conclude that the plot is not much changing when the number of grid points are doubled up and nearly tripled up. Although there is a sligth change at the throat area alone in the conservative form, which is not considerably significant.