Quasi 1D flow through a nozzle
Main code where the functions are called:
%1D nozzle flow
clear all
close all
clc
%Grid spacing:
n = 31;
x = linspace(0,3,n);
dx = x(2) - x(1);
g = 1.4;
A = 1 + 2.2*((x-1.5).^2);
C = 0.5;
nt = 1000;
%Calculating location at throat
th = (n+1)/2;
tic
[vc,rho_c,Tc,mc,Mc,pc] = conservative(A, g, n, dx, nt, x, C);
t1 = toc
tic
[vnc,rho_nc,Tnc,mnc,Mnc,pnc] = non_conservative(A, g, n, dx, nt, x, C);
t2 = toc
figure(3)
plot(x, mnc)
hold on
plot(x,mc)
title(sprintf(\'Comparison of non-dimensional mass flow rates\'));
legend(\'Non-conservative form\', \'Conservative form\');
figure(4)
subplot(2,1,1)
plot(x, A/2, \'b\')
hold on
plot(x,-A/2, \'b\')
title(sprintf(\'Non-conservative form primitive variables at the %dth time step\', nt));
subplot(2,1,2)
hold on
plot(x,Mnc)
plot(x,rho_nc)
plot(x,pnc)
plot(x,Tnc)
legend(\'Mach number\', \'Non-dimensional Density\', \'Non-dimensional pressure\', \'Non-dimensional temperature\')
title(sprintf(\'Non-conservative form primitive variables at the %dth time step\', nt));
figure(5)
subplot(2,1,1)
plot(x, A/2, \'b\')
hold on
plot(x,-A/2, \'b\')
title(sprintf(\'Conservative form primitive variables at the %dth time step\', nt));
subplot(2,1,2)
hold on
plot(x,Mc)
plot(x,rho_c)
plot(x,pc)
plot(x,Tc)
hold off
legend(\'Mach number\', \'Non-dimensional Density\', \'Non-dimensional pressure\', \'Non-dimensional temperature\')
%Calculating the values of primitive variables at the throat
terms = {\'Non-dimensional density\'; \'Non-dimensional Temperature\'; \'Non-dimensional pressure\'; \'Non-dimensional Mach no.\'};
Exact_solution = [0.634;0.833;0.528;1];
Conservative_terms = [rho_c(th); Tc(th); pc(th); Mc(th)];
Non_conservative_terms = [rho_nc(th); Tnc(th); pnc(th); Mnc(th)];
table(terms, Exact_solution, Conservative_terms, Non_conservative_terms)The function to solve the conservation equations:
function [v,rho,T,U2,M,p] = conservative(A, g, n, dx, nt, x, C)
%Initial conditions:
for l = 1:n
if (0<=x(l))&&(x(l)<0.5)
rho(l) = 1;
T(l) = 1;
end
if (0.5<=x(l))&&(x(l)<1.5)
rho(l) = 1 - 0.366*(x(l)-0.5);
T(l) = 1 - 0.167*(x(l)-0.5);
end
if (1.5<=x(l))&&(x(l)<=3)
rho(l) = 0.634-0.3879*(x(l)-1.5);
T(l) = 0.833 - 0.3507*(x(l)-1.5);
end
end
v = 0.59./(rho.*A);
U1 = rho.*A;
U2 = rho.*A.*v;
U3 = U1.*((T/(g-1))+(g/2)*(v.^2));
%Non-conservative form:
for k = 1:nt
%Defining the time step
a = T.^0.5;
dt = C*(dx./(a+v));
dt = min(dt);
%Updating the old values
U1old = U1;
U2old = U2;
U3old = U3;
%Determining the flux variables
F1 = U2;
term1 = U2.^2./U1;
term2 = (U3 - (g/2)*term1);
F2 = term1+((1/g)*term2*(g-1));
term3 = (U2.*U3)./U1;
term4 = (U2.^3)./(U1.^2);
F3 = ((g*term3)-(((g-1)*g)/2)*term4);
for i = 2:n-1
%Predictor step
dU1dt_p(i) = -(F1(i+1) - F1(i))/dx;
dU2dt_p(i) = -((F2(i+1) - F2(i))/dx)+((1/g)*(rho(i)*T(i))*(A(i+1)-A(i)))/dx;
dU3dt_p(i) = -(F3(i+1) - F3(i))/dx;
%Updating the values:
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
%Determining the predicted primitive variables
rho = U1./A;
v = U2./U1;
T = (g-1)*((U3./U1) - ((g/2)*(v.^2)));
%Corrector step
F1 = U2;
term1 = U2.^2./U1;
term2 = (U3 - (g/2)*term1);
F2 = term1+((1/g)*term2*(g-1));
term3 = (U2.*U3)./U1;
term4 = (U2.^3)./(U1.^2);
F3 = ((g*term3)-(((g-1)*g)/2)*term4);
for i = 2:n-1
dU1dt_c(i) = -(F1(i) - F1(i-1))/dx;
dU2dt_c(i) = -((F2(i) - F2(i-1))/dx)+((1/g)*(rho(i)*T(i))*(A(i)-A(i-1)))/dx;
dU3dt_c(i) = -(F3(i) - F3(i-1))/dx;
end
dU1dt = 0.5*(dU1dt_p + dU1dt_c);
dU2dt = 0.5*(dU2dt_p + dU2dt_c);
dU3dt = 0.5*(dU3dt_p + dU3dt_c);
%Final solution:
for j = 2:n-1
U1(j) = U1old(j) + dU1dt(j)*dt;
U2(j) = U2old(j) + dU2dt(j)*dt;
U3(j) = U3old(j) + dU3dt(j)*dt;
end
%Inlet boundary conditions
U1(1) = A(1);
U2(1) = 2*(U2(2)) - U2(3);
U3(1) = U1(1) * ((T(1)/(g-1)) + (g/2)*(v(1)^2));
%Outlet 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);
%Determining the primitive variables
rho = U1./A;
v = U2./U1;
T = (g-1)*((U3./U1) - (g/2)*(v.^2));
M = v./(T.^0.5);
p = rho.*T;
%Plotting the variation of massflow rate and non dimensional distance
%with respect to the time step
figure(1)
hold on
if k == 100
plot(x,U2, \'g\', \'linewidth\', 2)
else if k == 200
plot(x,U2, \'m\', \'linewidth\', 2)
else if k == 700
plot(x,U2, \'y\', \'linewidth\', 2)
title(sprintf(\'Variation of non-dimensional mass flow rate across the nozzle for conservative form\'));
legend(\'100dt\', \'200dt\', \'700dt\')
end
end
end
end
return
end
Function to solve non-conservation equation:
function [v,rho,T,m,M,p] = non_conservative(A, g, n, dx, nt, x, C)
rho = 1 - 0.3146*x;
T = 1 - 0.2314*x;
v = (0.1 + 1.09*x).*(T.^0.5);
for k = 1:nt
%Determining the time parameters
a = T.^0.5;
dt = C*(dx./(a+v));
dt = min(dt);
%Updating the old values
rho_old = rho;
vold = v;
Told = T;
for i = 2:n-1
%Predictor step
dvdx_p = (v(i+1) - v(i))/dx;
dlogAdx_p = (log(A(i+1)) - log(A(i)))/dx;
drhodx_p = (rho(i+1) - rho(i))/dx;
dTdx_p = (T(i+1) - T(i))/dx;
%Continuity equation:
drhodt_p(i) = -rho(i)*dvdx_p - rho(i)*v(i)*dlogAdx_p - v(i)*drhodx_p;
%Momentum equation:
dvdt_p(i) = -v(i)*dvdx_p - (1/g)*(dTdx_p + (T(i)/rho(i))*drhodx_p);
%Energy equation:
dTdt_p(i) = -v(i)*dTdx_p - (g-1)*T(i)*(dvdx_p + (v(i)*dlogAdx_p));
%Updating the values:
rho(i) = rho(i) + drhodt_p(i)*dt;
v(i) = v(i) + dvdt_p(i)*dt;
T(i) = T(i) + dTdt_p(i)*dt;
end
%Corrector step
for i = 2:n-1
dvdx_c = (v(i) - v(i-1))/dx;
dlogAdx_c = (log(A(i)) - log(A(i-1)))/dx;
drhodx_c = (rho(i) - rho(i-1))/dx;
dTdx_c = (T(i) - T(i-1))/dx;
%Continuity equation:
drhodt_c(i) = -rho(i)*dvdx_c - rho(i)*v(i)*dlogAdx_c - v(i)*drhodx_c;
%Momentum equation:
dvdt_c(i) = -v(i)*dvdx_c - (1/g)*(dTdx_c + (T(i)/rho(i))*drhodx_c);
%Energy equation:
dTdt_c(i) = -v(i)*dTdx_c - (g-1)*T(i)*(dvdx_c + (v(i)*dlogAdx_c));
end
drhodt = 0.5*(drhodt_p + drhodt_c);
dvdt = 0.5*(dvdt_p + dvdt_c);
dTdt = 0.5*(dTdt_p + dTdt_c);
%Final solution:
for j = 2:n-1
rho(j) = rho_old(j) + drhodt(j)*dt;
v(j) = vold(j) + dvdt(j)*dt;
T(j) = Told(j) + dTdt(j)*dt;
end
%Inlet boundary conditions
v(1) = 2*v(2) - v(3);
%Outlet boundary conditions
rho(n) = 2*rho(n-1) - rho(n-2);
v(n) = 2*v(n-1) - v(n-2);
T(n) = 2*T(n-1) - T(n-2);
m = rho.*A.*v;
M = v./(T.^0.5);
p = rho.*T;
%Plotting the variation of massflow rate and non dimensional distance
%with respect to the time step
figure(2)
hold on
if k == 50
plot(x, m, \'r\', \'linewidth\', 2)
else if k == 100
plot(x,m, \'b\', \'linewidth\', 2)
else if k == 150
plot(x,m, \'g\', \'linewidth\', 2)
else if k == 200
plot(x,m, \'m\', \'linewidth\', 2)
else if k == 700
plot(x,m, \'y\', \'linewidth\', 2)
title(sprintf(\'Variation of non-dimensional mass flow rate across the nozzle for non-conservative form\'));
legend(\'50dt\', \'100dt\', \'150dt\', \'200dt\', \'700dt\')
axis([0 3 0 2]);
end
end
end
end
end
end
return
end
The aim of this project is to model the flow through a subsonic inlet - supersonic outlet nozzle. The quasi 1D flow through the nozzle has been modelled using MacCormack's method. MacCormack's method is an explicit downstream marching method. In the case of the 1D nozzle flow, we march along the length of the nozzle from the inlet to the outlet.
2 studies are conducted, using two methods of solving the Navier-Stokes equations. We use the conservation and non-conservation methods and subsequently study the effectiveness of each method.
The MacCormack method has been used to solve the equations in both methods giving us a direct relation of the effectiveness of either method. The MacCormack method involves the use of a predictor corrector approach to solve the unsteady equations in an explicit manner. Since this is an explicit solver, we are marching forwards in time. With respect to space, we march forwards in space in the predictor step and backwards in space in the corrector step.
The nozzle is defined with a non-dimensional area of `A = 1 + 2.2(x-1.5)^2`. The nozzle is axi-symmetric and is a converging diverging nozzle.
Since the code is written with each method as a function, the code is pasted above for ease of understanding and the logic used is explained in the below sections.
All the results obtained are compared with the data provided in the John D Anderson jr book, 'Computational Fluid Dynamics'. This allows us to effectively understand the fundamental differences between the conservation and non conservation methods qualitatively.
Based on the aforementioned book, the quantity we are measuring for both methods is the mass flow rate.
It is also to be noted that all values have been normalised and are non dimensional for simplicity of the solution. The primitive variables are all normalised with the use of their reservoir values, the velocity is normalised with the speed of sound at the reservoir and the area is normalised with the area at the throat.
Non-dimensional temperature: `T' = T/T_0`
Non-dimensional density: `rho' = rho/rho_0`
Non-dimensional distance through nozzle: `x' = x/L`
Non-dimensional speed of sound: `a_0 = sqrt(gamma R T_0)`
Non-dimensional velocity: `V'= V/a_0`
Non-dimensional time: `t'= t/(L/a_0)`
Non-dimensional area: `A'= A/A^**`
The conservation method:
The conservation method, as the name suggests is one where the mass is better conserved over the geometry. This is especially intuitive while looking at the conservative form of the N-S equations where in the continuty, momentum and energy equations, the term rho*A*v is present. This in simple terms is the mass flow rate appending the earlier comment that mass flow rate is conserved along the geometry.
Apart from the mass flow rate, we are also concerned with the variation of the primitive variables across the nozzle.
The primitive variables follow the expected pattern with an increase in Mach number across the nozzle and decrease in temperature, pressure and density. For a nozzle without shock, this is the expected result and shows that we have modelled the flow precisely.
Variation of mass flow rate across the nozzle:
The results have been viewed at 3 distinct time steps so as to allow for convergence. The most striking information in this plot is the yellow line. This value was taken at the 700th time step and shows the converged mass flow rate. More importantly, it shows the mass flow rate as nearly a straight line, showing conservation of the mass across the nozzle.
The non-conservation method:
Theoretically, the mass flow rate isn\'t conserved as well in this method as compared to the non-conservation method. This is simply due to the fact that all variables vary across the nozzle as opposed to the conservation method where just the mass was varied across the nozzle.
Continuity equation:
`(drho)/dt = -rho ((dV)/dx) - rho V (d(ln A))/dx - V(drho)/dx`
Momentum equation:
`(dV)/dt = -V(dV)/dx - 1/gamma ( (dT)/dx + T/rho(drho)/dx)`
The variation of the primitive variables across the nozzle:
The primitive variables follow the same pattern in both cases leading to a conclusion that the mass flow rate will provide the better basis for a comparison between the two methods.
Variation of mass flow rate across the nozzle:
Immediately, it looks like the yellow line (data at the 700th time step) is a straight line as well. However, the scale for both the plots is different. While the scale for the conservative method plot was from 0.5 to 0.7, the scale for the non-conservative method is from 0 to 2. This shows 2 things, the data varies a lot for the non-conservative method and that we need to look at the mass flow rates on the same plot.
Comparison of mass flow rates achieved after 1000 time steps for both methods:
The plot comparing the mass flow rates shows us the actual variation in the mass flow rate in the non-conservative method as compared to the conservative method. The curious case here is the mass flow rate reducing for the non-conservative method at the throat and then increasing again. The conservative method shows a more linear plot. This led to studying the values of the primitive variables at the throat.
Variation of primitive variables at the throat:
The table shows the values of the primitive variables at the throat. The above mentioned phenomena is clearly visible here where the non-conservation method values are lower at the throat. Much lower, to the degree that the values are much closer to the exact solution. While mass was conserved, the primitive variables showed a large error. To understand this better, a grid dependency test was set up.
Grid dependency:
3 grids were chosen, one with 30, one with 60 and another with 90 grid spaces. The results obtained were then studied. The results obtained for 30 grid spaces has been attached above.
Results for 60 grid spaces:
Results for 90 grid spaces:
The results for both the conservative and non-conservative methods seem to tend towards the exact values. For our purposes these results are acceptable. The results for the grids don\'t vary much for the non-conservative method and grid independency can be assumed. However, the results for the conservative method changes a lot from the grid with 30 spaces to the grid with 60 spaces. To absolutely claim grid independency, another grid spacing can be assumed and ran. For this purpose, a grid spacing of 120 was assumed and the results so obtained are as below:
These results show us that while the conservative form is grid independent, the non-conservation form is dependent on the grid again. This trend continues as the grid size is increased with the values going to unimaginable values. However, the conservative form holds the independency and provides a more accurate results with a more fine grid spacing.
The above table is for a very fine spacing of 120000 spaces. This trend continues as the grid becomes finer and can be explained as the primitive variables themselves being changed at each grid point in every iteration in the non-conservation method. This causes small errors in each value to be magnified in finer grids giving us an unimaginable solution. While the results for the conservative method holds at fine grid spacing, it isn\'t so for the non-conservative method.
The final test in the evaluation of conservation method vs non-conservation method, we look at the time taken for the simulation to run for 1000 time steps at different grid spacings.
The test was conducted for grid spacings of 30, 60, 90, 120 and 120000.
Convergence check:
Convergence was checked for both methods where 2 variables were checked for convergence to an order of 1e-4. The results yielded were similar for both methods and hence the time taken for completion of the simulation is calculated below.
Time taken for the simulation to run for 1000 time steps:
For 30 grid spaces:
t1 = 3.23
t2 = 2.30
For 60 grid spaces:
t1 = 3.9349
t2 = 2.9170
For 90 grid spaces:
t1 = 4.3848
t2 = 3.1531
For 120 grid spaces:
t1 = 4.1025
t2 = 3.1038
For 120000 grid spaces:
t1 = 62.3304
t2 = 49.6431
Where t1 and t2 are the times taken for the simulation to complete for conservation method and non-conservation method respectively. This clearly shows that the non-conservative method is always faster regardless of the grid used. However, for finer grids the results obtained using the non-conservation method are not satisfactory. This leads to the question, which method is more suitable. There are good and bad points to both and the answer depends on the user. A user who wants the simulation time and coding time to be less can choose the non-conservation form constraint to the size of the mesh and one who requires the mass be conserved can use the conservation method subject to more coding and an increase in the simulation time.