Week 7 - Simulation of a 1D Super-sonic nozzle flow simulation using Macormack Method
Objective:
To simulate the isentropic flow through a 1D subsonic - supersonic nozzle using both conservative and non conservative forms of the governing equations and solving them using the Maccormack technique, and analyse the results.
Solution
Theory
We conside the isentropic flow through a convergent - divergent nozzle. The flow at the inlet comes from a reservoir. The cross sectional area of the reservoir is large and the velocity is correspondingly small.
The flow is subsonic (Mach number <1) at the convergent section,sonic at the throat (Mach numer = 1) and supersonic at the divergent section (M>1)
Considering a small control volume as shown below
The cross sectional area is A, and it is assumed that the flow properties are uniform across the section.
Although the area of the nozzle changes as a function of distance along X-direction, the flow field in reality is two dimensional. The flow varies along X and Y directions. We assume that the flow varies only along X - direction.
Non conservative form
The governing continuity, momentum and energy equations can be expressed in non-conservative form as
Continuity equation
`(delrho')/(delt') = -rho'((delV')/(delx')) - rho'V'((del(ln A'))/(delx')) - V'((delrho')/(delx'))`
Momentum equation
`(delV')/(delt') = -V'((delV')/(delx')) -(1/gamma)((delT')/(delx') + (T')/(rho')((delrho')/(delx')))`
Energy equation
`(delT')/(delt') = -V'((delT')/(delx')) - (gamma -1)T'[(delV')/(delx') + V'(del(lnA'))/(delx')]`
Conservative form
Continuity equation
`(del(rho'A'))/(delt') + (del(rho'A'V'))/(delx') = 0`
Momentum equation
`(del(rho'A'V'))/(delt') + (del[rho'A'(V')^2 +(1/gamma)rho'A'])/(delx') = (1/gamma)rho'((delA')/(delx'))`
Energy equation
`(del[rho'((e')/(gamma-1) + (gamma/2)(V')^2)A'])/(delt') + (del[rho'((e')/(gamma-1) + (gamma/2)(V')^2)V'A' + rho'A'V'])/(delx') = 0`
For air at standard conditions, `gamma = 1.4`
Maccormack method
The steps involved are
1. Initial profile
The initial profile for all the flow field variables - velocity, density and temperature has to be determined.
2. Predictor step
`((delf)/(dt))_p = (f(i+1) - f(i))/(Deltax`
3. Solution update
`f_p = f_(Initial) + ((delf)/(delt))_p * Deltat`
`f_p =` Flow field variables after predictor step
`f_(Initial) =` flow field variable in the initial step
`Deltat =` time step size
4. Corrector step
`((delf_p)/(dt))_c = (f(i) - f(i-1))/(Deltax`
`((delf)/(dt))_(avg) = 0.5(((delf)/(dt))_p +((delf_p)/(dt))_c )`
5. Final solution update
`f_(sol) = f_(Initial) + ((delf)/(delt))_(avg) * Deltat`
`f_(sol) =` Final solution of the flow field variabes for first time step
Matlab code
Main program
clear all
close all
clc
L = 3;
n = [31 62 93];
nt = 1:1400;
cfl = 0.5;
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 - Conservative vs non conservative form at throat point');
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 - Conservative vs non 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');
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');
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)
% inputs
x = linspace(0,3,n);
dx = x(2) - x(1);
% initial profiles
rho = 1-0.3146*x;
t = 1-0.2314*x;
v = (0.1+0.9*x).*t.^0.5;
% area
a = 1+2.2*(x-1.5).^2;
M = v./(t.^0.5);
mfr = rho.*a.*v;
p = rho.*t;
% Time steps
for z = 1:31
u = (t(z))^(1/2);
dt(z) = cfl*dx/(u+v(z));
end
dt = min(dt);
gamma = 1.4;
% Outer loop - time
for k = 1:length(nt)
rho_old = rho;
v_old = v;
t_old = t;
% Predictor
for j = 2:n-1
dv_dx = (v(j+1) - v(j))/dx;
drho_dx = (rho(j+1)-rho(j))/dx;
dt_dx = (t(j+1) - t(j))/dx;
dloga_dx = (log(a(j+1)) - log(a(j)))/dx;
% Continuity eqn
drho_dt_p(j) = -rho(j)*dv_dx - rho(j)*v(j)*dloga_dx - v(j)*drho_dx;
% Momentum eqn
dv_dt_p(j) = -v(j)*dv_dx - (1/gamma)*(dt_dx + (t(j)/rho(j))*drho_dx);
% Energy equation
dt_dt_p(j) = -v(j)*dt_dx - (gamma-1)*t(j)*(dv_dx + v(j)*dloga_dx);
% Solution update
rho(j) = rho_old(j) + drho_dt_p(j)*dt;
v(j) = v_old(j) + dv_dt_p(j)*dt;
t(j) = t_old(j) + dt_dt_p(j)*dt;
% Corrector method
dv_dx = (v(j) - v(j-1))/dx;
drho_dx = (rho(j)-rho(j-1))/dx;
dt_dx = (t(j) - t(j-1))/dx;
dloga_dx = (log(a(j)) - log(a(j-1)))/dx;
% Continuity eqn
drho_dt_c(j) = -rho(j)*dv_dx - rho(j)*v(j)*dloga_dx - v(j)*drho_dx;
% Momentum eqn
dv_dt_c(j) = -v(j)*dv_dx - (1/gamma)*(dt_dx + (t(j)/rho(j))*drho_dx);
% Energy equation
dt_dt_c(j) = -v(j)*dt_dx - (gamma-1)*t(j)*(dv_dx + v(j)*dloga_dx);
end
% Average time derivative
drho_dt = 0.5*(drho_dt_p + drho_dt_c);
dv_dt = 0.5*(dv_dt_p + dv_dt_c);
dt_dt = 0.5*(dt_dt_p + dt_dt_c);
% Final soltion update
for i = 2:n-1
rho(i) = rho_old(i) + drho_dt(i)*dt;
v(i) = v_old(i) + dv_dt(i)*dt;
t(i) = t_old(i) + dt_dt(i)*dt;
end
% Boundary condition
% Inlet
v(1) = 2*v(2) - v(3);
% 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);
%Primitives update
p = rho.*t;
M = v./t.^0.5;
error(i) = abs(v(i) - v(i-1));
mfr = rho.*a.*v;
%values at throat point
t_throat(k) = t(16);
rho_throat(k) = rho(16);
v_throat(k) = v(16);
p_throat(k) = p(16);
M_throat(k) = M(16);
error_throat(k) = error(16);
mfr_throat(k) = mfr(16);
if n==31
for p = [1:50:201 700]
if p == k
figure(1)
plot(x,mfr,'linewidth',2);
hold on
title('Non conservative form - Varitation of mass flow rate at different time steps');
xlabel('Length of nozzle');
ylabel('Mass flow rate');
elseif k == 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');
title('Non conservative form - Steady state distriubution 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 primitive variables at throat 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);
%Inputs
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;
gamma = 1.4;
% Solution vectors u1, u2 and u3
u1 = rho.*a;
u2 = 0.59*ones(1,n);
v = u2./u1;
u3 = u1.*((T./(gamma-1))+(gamma/2).*v.^2);
M = v./(T.^0.5);
mfr = rho.*a.*v;
p = rho.*T;
for z = 1:31
u = (T(z))^(1/2);
dt(z) = cfl*dx/(u+v(z));
end
dt = min(dt);
for k = 1:length(nt)
%Initial values
u1_old = u1;
u2_old = u2;
u3_old = u3;
for j = 2:n-1
% Flux vectors
f1 = u2;
f2 = ((u2.^2)./u1) + ((gamma-1)/gamma)*(u3 - ((gamma/2)*(u2.^2)./u1));
f3 = ((gamma.*u2.*u3)./u1) - ((gamma*(gamma-1)/2)*(u2.^3)./(u1.^2));
% Source term
j1 = (1/gamma)*rho(j)*T(j)*((a(j+1) - a(j))/dx);
% Predictor step
du1_dt_p(j) = -(f1(j+1) - f1(j))/dx;
du2_dt_p(j) = -(f2(j+1) - f2(j))/dx + j1;
du3_dt_p(j) = -(f3(j+1) - f3(j))/dx;
% 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 = (gamma - 1).*(u3./u1-(gamma/2).*(u2./u1).^2);
% Updating flux terms
f1 = u2;
f2 = u2.^2./u1 + ((gamma-1)/gamma).*(u3 - (gamma/2)*u2.^2./u1);
f3 = gamma.*u2.*u3./u1 - 0.5*gamma*(gamma-1).*u2.^3./u1.^2;
% Source term for corrected step
j2 = (1/gamma)*(rho(j)*T(j)*(a(j) - a(j-1))/dx);
% Corrector step
du1_dt_c(j) = -(f1(j) - f1(j-1))/dx;
du2_dt_c(j) = -(f2(j) - f2(j-1))/dx + j2;
du3_dt_c(j) = -(f3(j) - f3(j-1))/dx;
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 corrected values
for i = 2:n-1
u1(i) = u1_old(i) + du1_dt_avg(i)*dt;
u2(i) = u2_old(i) + du2_dt_avg(i)*dt;
u3(i) = u3_old(i) + du3_dt_avg(i)*dt;
end
% Boundary conditions
% Inlet
u1(1) = a(1);
u2(1) = 2*u2(2) - u2(3);
u3(1) = u1(1)*(T(1)/(gamma-1) + (gamma*(u2(1)/u1(1))^2)/2);
% Outlet
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 = (gamma-1).*(u3./u1 - (gamma/2).*(u2./u1).^2);
p = rho.*T;
M = v./T.^0.5;
error(i) = abs(v(i) - v(i-1));
mfr = u2;
%values at throat point
t_throat(k) = T(16);
rho_throat(k) = rho(16);
v_throat(k) = v(16);
M_throat(k) = M(16);
p_throat(k) = p(16);
error_throat(k) = error(16);
mfr_throat(k) = mfr(16);
if n==31
for p = [1:50:201 700]
if p == k
figure(4)
plot(x,mfr,'linewidth',2);
hold on
title('conservative form - Varitation of mass flow rate at different time steps');
xlabel('Length of nozzle');
ylabel('Mass flow rate');
elseif k == 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');
title('Conservative form - Steady state distriubution 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 primitive variables at throat i = 16');
end
mfr_c = mfr;
mfr_throat_c = mfr_throat;
endResults
Non Conservative form
1. Mass flow rate
It can be observed that the mass flow rate stabilizes as the number of time steps are increased
2. Non dimensional parameters
The figure above shows the variation of all the calculated flow field variables over the length of the nozzle.
3. Primitive variables - Time wise variation
The primitive variables calculated in the flow field stabilize and achieve convergence at around 500 steps. It can be seen, that the values remain constant thereafter.
Conservative form
1. Mass flow rate
It can be observed that the mass flow rate stabilizes as the number of time steps are increased
2. Non dimensional parameters
The figure above shows the variation of all the calculated flow field variables over the length of the nozzle.
3. Primitive variables - Time wise variation
The primitive variables calculated in the flow field stabilize and achieve convergence at around 500 steps. It can be seen, that the values remain constant thereafter.
Comparison of mass flow rate
It can be seen, that the non conservative form takes a longer time to reach stable convergence, compared to conservative form at the throat point.
Over the length of the nozzle, the mass flow rate is constant for conservative form. There is some variation in the mass flow rate of the non conservative form, but it is very minor (order of 10^-3).
Grid independence test
The simulations are carried out with three different grid sizes - 31,62 and 93 points in order to illustrate the independence of the solution on the grid size.
The results can be seen in the images below.
Nozzle
It can be seen, that a stable result is obtained for all three grid sizes. There is some variation in the actual mass flow rate. As the size of the grid increases, the mass flow rate reduces slightly for both conservative and non conservative flows
Throat point
The mass flow rates for all three grid sizes reach the same constant value, therefore ensuring stability. The simulation takes more time steps to reach this convergence, as the grid size increases. Overall, a stable solution is reached for all three grid sizes.
Conclusion
Thus, the 1D supersonic flow over a nozzle has been simulated for both conservative and non conservative forms, and the results have been analyzed and compared.
The matlab m files have been attached below