Simulation of Quasi 1-D Supersonic flow through a nozzle
OBJECTIVE:
(1) Solve the Conservative and Non-conservative form of Governing equations of Quasi 1-D Supersonic Nozzle flow using the Maccormack method.
(2) Compare the results of the Conservative and Non-Conservative form.
(3) Perform grid dependency test.
THEORY:
- The quasi 1-D Flow through a nozzle has been modelled using MacCormack\'s method which is an explicit downstream marching method where in the in the case of this experiment we march along the length of the nozzle from the inlet to the outlet.
- The values have been normalized and are equated and plotted in terms of non-dimensional forms for simplicity of the solution. The Field variables are normalized using their reservoir values, the velocity using the speed of sound at the reservoir and Area with the throat cross sectional area.
- The following 2 case setups studied, using 2 methods of Navier-Stokes equation gives us the corresponding comparison of its effectiveness.
- The MacCormack method uses a Predictor-Corrector approach to solve the unsteady equations in an explicit manner.
- The Converging-Diverging (CD) nozzle is defined with non-dimensional area of `A = 1+2.2(x-1.5)^2`, is Axi-symmetric in nature
- The results obtained are concurrent with data provided in the \'Computational Fluid Dynamics - Basics with Application by John D Anderson Jr. This allows us to effectively understand the quantifiable differences between the Conservation and Non-Conservation forms.
MATLAB CODE:
- MAIN FUNCTION
close all
clear all
clc
n = 100;
gamma =1.4;
x = linspace(0,3,n);
dx = x(2) - x(1);
nt = 1400;
CFL = 0.5;
A = 1+2.2*(x-1.5).^2;
%Calculating throat location
th = (n+1)/2;
tic
[v_c,rho_c,T_c,m_c,M_c,p_c] = conservativeform(A,gamma,n,dx,nt,x,CFL);
t1 = toc;
tic
[v_nc,rho_nc,T_nc,m_nc,M_nc,p_nc] = Nonconservativeform(A,gamma,n,dx,nt,x,CFL);
t2 = toc;
figure(3)
plot(x,m_nc)
hold on
plot(x,m_c)
title(sprintf(\'Non-dimensional mass flow rates comparison\'));
legend(\'Non-conservative form\',\'Conservative form\');
figure(4)
subplot(4,1,1)
hold on
plot(x,M_nc)
xlabel(\'Length of nozzle\')
ylabel(\'Mach Number\')
grid minor
subplot(4,1,2)
plot(x,rho_nc)
xlabel(\'Length of nozzle\')
ylabel(\'Density ratio\')
grid minor
subplot(4,1,3)
plot(x,p_nc)
xlabel(\'Length of nozzle\')
ylabel(\'Pressure ratio\')
grid minor
subplot(4,1,4)
plot(x,T_nc)
xlabel(\'Length of nozzle\')
ylabel(\'Temperature ratio\')
grid minor
sgtitle(\'Qualitative aspects of Quasi 1-D Supersonic Nozzle Flow\')
figure(5)
subplot(4,1,1)
hold on
plot(x,M_c)
xlabel(\'Length of nozzle\')
ylabel(\'Mach Number\')
grid minor
subplot(4,1,2)
plot(x,rho_c)
xlabel(\'Length of nozzle\')
ylabel(\'Density ratio\')
grid minor
subplot(4,1,3)
plot(x,p_c)
xlabel(\'Length of nozzle\')
ylabel(\'Pressure ratio\')
grid minor
subplot(4,1,4)
plot(x,T_c)
xlabel(\'Length of nozzle\')
ylabel(\'Temperature ratio\')
grid minor
sgtitle(\'Qualitative aspects of Quasi 1-D Supersonic Nozzle Flow\')
- CONSERVATIVE FORM
function [v,rho,T,m,M,p] = conservativeform(A,gamma,n,dx,nt,x,CFL)
m = zeros(1,n);
rho = zeros(1,n);
T = zeros(1,n);
%Initial conditions
rho(1:6) = 1;
T(1:6) = 1;
rho(7:16) = 1-0.366*(x(7:16)-0.5);
T(7:16) = 1-0.167*(x(7:16)-0.5);
rho(17:n) = 0.634-0.3879*(x(17:n)-1.5);
T(17:n) = 0.833-0.3507*(x(17:n)-1.5);
v = 0.59./(rho.*A); %Local velocity
%Solution vectors
u1 = rho.*A;
u2 = rho.*A.*v;
u3 = rho.*(((1/(gamma-1)).*T) + ((gamma/2)*v.^2)).*A;
u1_old = u1;
u2_old = u2;
u3_old = u3;
for k=1:nt
%Time domain
a = sqrt(T); % Local speed of sound
dt_calc = (CFL*dx)./(a + v); %(CFL X dx)/(a+V)
dt = min(dt_calc);
%Boundary conditions
%Inlet
u1(1) = A(1);
u2(1) = (2*u2(2))-u2(3);
u3(1) = (u1(1)*((T(1)/0.4))) + (0.7*v(1).^2*(A(1)));
%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);
%Flux variables for predictor
f1 = u2;
term1 = ((u2.^2)./u1);
term2 = (0.4/gamma).*(u3-(0.7.*(u2.^2./u1)));
f2 = term1 + term2;
f3 = ((gamma.*u2.*u3)./u1)-((0.7*0.4).*(u2.^3./u1.^2));
%%Predictor step
for i=2:n-1
%Source term
j2(i) = (1/gamma).*rho(i).*T(i).*((A(i+1)-A(i))./dx);
%Continuity equation
du1dt_p(i) = -(f1(i+1)-f1(i))/dx;
%Momentum equation
du2dt_p(i) = -((f2(i+1)-f2(i))/dx)+j2(i);
%Energy equation
du3dt_p(i) = -(f3(i+1)-f3(i))/dx;
%Step update
u1(i) = u1(i)+du1dt_p(i).*dt;
u2(i) = u2(i)+du2dt_p(i).*dt;
u3(i) = u3(i)+du3dt_p(i).*dt;
end
%Updating field variables
rho = u1./A;
v = u2./u1;
T = 0.4.*((u3./u1)-(0.7.*v.^2));
%%Corrector step
%Flux variables for corrector step
f1 = u2;
term1 = ((u2.^2)./u1);
term2 = (0.4/gamma).*(u3-(0.7.*(u2.^2./u1)));
f2 = term1 + term2;
f3 = ((gamma.*u2.*u3)./u1)-((0.7*0.4).*(u2.^3./u1.^2));
for i=2:n-1
j2(i) = (1/gamma).*rho(i).*T(i).*((A(i)-A(i-1))./dx);
du1dt_c(i) = -(f1(i)-f1(i-1))/dx;
du2dt_c(i) = -((f2(i)-f2(i-1))/dx)+j2(i);
du3dt_c(i) = -(f3(i)-f3(i-1))/dx;
end
%Compute average time derivative
du1dt = 0.5*(du1dt_p+du1dt_c);
du2dt = 0.5*(du2dt_p+du2dt_c);
du3dt = 0.5*(du3dt_p+du3dt_c);
%%Final Solution update
for i=2:n-1
u1(i) = u1_old(i) + du1dt(i)*dt;
u2(i) = u2_old(i) + du2dt(i)*dt;
u3(i) = u3_old(i) + du3dt(i).*dt;
end
%Updating field variables
rho = u1./A;
v = u2./u1;
T = 0.4.*((u3./u1)-(0.7.*v.^2));
m = rho.*A.*v;
p = rho.*T; %Local Pressure-ratio
M = v./sqrt(T); %Local Mach number
u1_old = u1;
u2_old = u2;
u3_old = u3;
%Plot of variation of MAF along the non dimensional distance at different
%time steps
figure(1)
hold on
grid minor
if k==10
plot(x,m_c,\'-\',\'linewidth\',2)
M1 = \"10dt\" ;
else if k==50
plot(x,m_c,\'g\',\'linewidth\',2)
M2 = \"50dt\";
else if k==100
plot(x,m_c,\'r\',\'linewidth\',2)
M3 = \"100dt\";
else if k==400
plot(x,m_c,\'b\',\'linewidth\',2)
M4 = \"500dt\";
else if k==1200
plot(x,m_c,\'black\',\'linewidth\',2)
M5 = \"1200dt\";
xlabel(\'Non-dimensional distance through nozzle\')
ylabel(\'Non-dimensional mass air flow rate\')
title(sprintf(\'Distribution of Mass air flow rate in the nozzle at different time instances\'))
legend([M1,M2,M3,M4,M5,M6]);
axis([0 3 0 1.2]);
end
end
end
end
end
end
end- NON-CONSERVATIVE FORM
function [v,rho,T,m,M,p] = Nonconservativeform(A,gamma,n,dx,nt,x,CFL)
%Initial conditions
T = 1-0.2314*x; %Temperature
v = (0.1 + 1.09*x).*T.^0.5; %Local velocity
rho = 1-0.3146*x; %Density
%Time domain
a = sqrt(T); % Local speed of sound
dt_calc = (CFL*dx)./(a + v); %(CFL X dx)/(a+V)
dt = min(dt_calc);
%%MacCormack method
tic
for k=1:nt %Time loop
rho_old = rho;
T_old = T;
v_old = v;
%%Predictor step
for i=2:n-1 %Forward differencing of spatial terms
drhodx_p(i) = (rho(i+1)-rho(i))/dx;
dvdx_p(i) = (v(i+1)-v(i))/dx;
dTdx_p(i) = (T(i+1)-T(i))/dx;
dlogAdx_p(i) = (log(A(i+1))-log(A(i)))/dx;
%Continuity equation
drhodt_p(i) = -rho(i)*dvdx_p(i)-rho(i)*v(i)*dlogAdx_p(i)-v(i)*drhodx_p(i);
%Momentum equation
dvdt_p(i) = -v(i)*dvdx_p(i)-(1/gamma)*(dTdx_p(i)+(T(i)/(rho(i)))*drhodx_p(i));
%Energy equation
dTdt_p(i) = -v(i)*dTdx_p(i)-(gamma-1)*T(i)*(dvdx_p(i)+v(i)*dlogAdx_p(i));
%Step update
rho(i) = rho(i)+(drhodt_p(i)*dt);
T(i) = T(i)+(dTdt_p(i)*dt);
v(i) = v(i)+(dvdt_p(i)*dt);
end
%%Corrector method
for i=2:n-1 %Backward differencing of spatial terms
dvdx_c(i) = (v(i)-v(i-1))/dx;
drhodx_c(i) = (rho(i)-rho(i-1))/dx;
dTdx_c(i) = (T(i)-T(i-1))/dx;
dlogAdx_c(i) = (log(A(i))-log(A(i-1)))/dx;
%Continuity equation
drhodt_c(i) = -rho(i)*dvdx_c(i)-rho(i)*v(i)*dlogAdx_c(i)-v(i)*drhodx_c(i);
%Momentum equation
dvdt_c(i) = -v(i)*dvdx_c(i)-(1/gamma)*(dTdx_c(i)+(T(i)/rho(i))*drhodx_c(i));
%Energy equation
dTdt_c(i) = -v(i)*dTdx_c(i)-(gamma-1)*T(i)*(dvdx_c(i) + v(i)*dlogAdx_c(i));
end
%Compute average time derivative
drhodt = 0.5*(drhodt_p + drhodt_c);
dvdt = 0.5*(dvdt_p + dvdt_c);
dTdt = 0.5*(dTdt_p + dTdt_c);
%Solution update
for m=1:n-1
rho(m) = rho_old(m)+drhodt(m)*dt;
v(m) = v_old(m)+dvdt(m)*dt;
T(m) = T_old(m)+dTdt(m)*dt;
end
%Boundary conditions
%Inlet BC
v(1) = 2.*v(2)-v(3);
%Outlet BC
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);
% %Pressure-ratio
p = rho.*T;
% %Mach number
M = v./sqrt(T);
% %Mass air flow rate(MAF)
m = rho.*v.*A;
% Plot of variation of MAF along the non dimensional distance at different
% time steps
figure(2)
hold on
grid minor
if k==10
plot(x,m_nc,\'-\',\'linewidth\',2)
M1 = \"10dt\" ;
else if k==50
plot(x,m_nc,\'g\',\'linewidth\',2)
M2 = \"50dt\";
else if k==100
plot(x,m_nc,\'r\',\'linewidth\',2)
M3 = \"100dt\";
else if k==400
plot(x,m_nc,\'b\',\'linewidth\',2)
M4 = \"400dt\";
else if k==800
plot(x,m_nc,\'yellow\',\'linewidth\',2)
M5 = \"800dt\";
else if k==1200
plot(x,m_nc,\'black\',\'linewidth\',2)
M6 = \"1200dt\";
xlabel(\'Non-dimensional distance through nozzle\')
ylabel(\'Non-dimensional mass air flow rate\')
title(sprintf(\'Distribution of Mass air flow rate in the nozzle at different time instances\'))
legend([M1,M2,M3,M4,M5,M6]);
axis([0 3 0 2]);
end
end
end
end
end
end
end
endCONSERVATIVE METHOD :
The conservation form, as the name suggests the one where mass is better conserverd over the geometry. This is particularly intuitive when looking at the conservation form of the N-S equations where in the continuity, momentum and the energy equations, the term `rhoAv` is present.
Variation of MAF across nozzle at different time steps:
The results have been viewed at different dts so as to allow convergence. At around 1200dt the MAF is a straight line which shows its convergence.
Apart from the mass flow rate we are also concerned with the variation of Field variables across the nozzle.
Variation of Field variables across nozzle:
The Field variables follow the expected pattern with an increase of Mach number and decrease of Temperature, Density and Pressure-ratio, shows the simulation has been modelled precisely.
NON-CONSERVATIVE METHOD :
Theoretically, the mass flow rate isnt conserved as much as it is in conservative method due to the fact that all variables vary across the length as opposed to the varying MAF alone acorss the nozzle.
Variation of MAF across nozzle at different time steps:
Variation of Field variables across nozzle:
The Field variables also follow the same pattern stating that MAF will provide as a better alternative in terms of methods of employment.
Variation of Field variables across nozzle:
The plot comparing the MAF of these two forms shows the MAF reducing for conservation till the throat and then increasing again. The conservation follows an almost linear path comparatively.
GRID-DEPENDENCY:
3 grids were chosen ( n = 31,61 & 91) . The results obtained were then studied. The values shown are at the throat section.
- 30 GRIDS:
rows Exact_Solution Conservative_method Non_Conservative_method
______________ ______________ ___________________ _______________________
\'Density ratio\' 0.634 0.65048 0.63866
\'Pressure ratio\' 0.833 0.84009 0.83645
\'Temperature ratio\' 0.528 0.54646 0.53421
\'Mach Number\' 1 0.98156 0.99939
The non conservational form gives much accurate results for 30 grid points than conservational form, goes on to show that when mass is conserved the field variables show deviations.
- 60 GRIDS
rows Exact_Solution Conservative_method Non_Conservative_method
_________________ ______________ ___________________ _______________________
\'Density ratio\' 0.634 0.6509 0.63766
\'Pressure ratio\' 0.833 0.84338 0.83544
\'Temperature ratio\' 0.528 0.54896 0.53273
\'Mach Number\' 1 0.9733 0.99978
- 90 GRIDS
rows Exact_Solution Conservative_method Non_Conservative_method
________________ ______________ ___________________ _______________________
\'Density ratio\' 0.634 0.63653 0.64044
\'Pressure ratio\' 0.833 0.83465 0.83681
\'Temperature ratio\' 0.528 0.53128 0.53592
\'Mach Number\' 1 1.0013 1.0026
From the table we can see that, increase in the grid spaces makes the FV approach exact solution for conservative method compared to non-conservative method. Hence we can conclude that Conservation form is grid dependant and Non-conservative form is grid independant. However increase in grid points increase the computational time for the simulation hence a suitable compromise needs to be selected in between accuracy and computaional cost.