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.
Answer:
Let us consider a small control volume of the nozzle ,
We assume that the flow parameters are consistent over a particular segment with a cross section area of A. As a result, even though the area of the nozzle changes as a function of distance along the nozzle, X, and thus the flow field is two dimensional (the flow varies in the two dimensional x-y space), we assume that the flow properties vary only with X; this is equivalent to assuming uniform flow properties across any given cross-section. This type of flow is known as quasi one-dimensional flow.
In terms of non-dimensional variables, the governing continuity, momentum, and energy equations for this quassi one-dimensional staedy and isentropc flow are as follows:
These equations are solved by Macormack method using time marching process.
Macormack Method:
Let f(x) is the equation to be solved.Then the following steps are followed,
i. inital profile
Before the calculation, the initial profile for all the flow field variables which are velocity, density and temperature has to be determined. These are guess values but should be specified with more awareness on the problem.
ii. predector step
iii. solution update
iv. Corrector step
v. Final solution update
Function code:
Non-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;
end
Conservative 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;
end
Matlab Code
OUTPUT:
Non conservative form:
1.
2.
3.
Figure (1) depicts the steady state distribution of nondimensional primitive variables at the conclusion of a time step. These flow field variable settings compel the precise values shown in John D. Anderson's book.
Figure (2) depicts the time-dependent fluctuation of basic variables at throat point i =16. The flow field variables reached a stable solution at time step 500.
Figure (3) depicts the change in mass flow rate distribution inside the nozzle at various time steps throughout the time marching procedure. The mass flow rate eventually stabilised at 0.58.
Conservative form:
4.
5.
6.
Figure (4) depicts the steady state distribution of nondimensional primitive variables at the conclusion of a time step. These flow field variable settings compel the precise values shown in John D. Anderson's book.
Figure (5) depicts the time-dependent fluctuation of basic variables at throat point i =16. The flow field variables reached a stable solution at time step 600.
The fluctuation of mass flow rate distribution inside the nozzle at different time steps throughout the time marching operation is seen in fig(6). The mass flow rate eventually stabilised at 0.56.
Non Conservative form vs Conservative form:
7.
Figure (7) compares mass flow rate using conservative and non-conservative formulations.
Subplot 1: By the conclusion of the time steps, both forms had reached a stable and perfect solution. However, the conservative variant preserved stability at smaller time steps.
Subplot 2 shows that the difference between the two shapes is quite small in the nozzle. However, in the conservative form, the mass flow rate is constant, but in the non-conservative version, it varies by 0.5820.003.
We may conclude from fig(7) that the final solution does not differ much when conservative and non-conservative forms are used. Although there is a slight shift in the conservative form in the neck area alone, it is not statistically significant.
Grid independency test for both Non-conservative form & Conservative form:
8.
9.
Figures (8) and (9) depict the grid independence test in both formats.
fig(8) - mass flow rate varies with time step at neck point - conservative form (subplot(1)) & non-conservative form (subplot(2)).
Though the values are varied at the start of the time marching procedure, the solution obtained is stable and almost identical by the completion of the time marching.
fig(9) - mass flow rate in the nozzle - conservative form (subplot(1)) & non-conservative form (subplot(2)).
The solution found using three different grid widths yielded three distinct solutions within the range of 0.5890.006. That is, the difference between them is so minor that it may be overlooked.
We may deduce from figs. 8 and 9 that the plot does not change significantly when the number of grid points is doubled and almost tripled. Although there is a slight shift in the conservative form in the neck area alone, it is not statistically significant.