Solving Quasi 1D Supersonic Flow through a Convergent Divergent Nozzle using MATLAB
Abstract : Set up of a time marching, finite difference solution for the Quasi 1D subsonic supersonic nozzle flow. Non-dimensionalizing the Governing flow equations and setting them up with relevant initial and boundary conditions and solving using the Maccormacks technique. Simulating the solution for convergence and comparing the normalized mass flow rate between Non-conservative and Conservative forms and making conclusions. Performing grid dependence test.
I. Introduction :
We consider the steady state, isentropic flow through convergent-divergent nozzle as shown below. The flow at inlet to the nozzle comes from a reservoir. The cross-sectional area of the reservoir is large and hence the velocity is very small. Thus pressure and temperature at reservoir are stagnation values or total pressure and total temperatures, respectively. The flow expands isentropically to super-sonic speeds at nozzle exit. The flow is locally subsonic in the convergent section of the nozzle and sonic at the throat and supersonic at the divergent section. The sonic flow is M=1 , where M is mach no. means that local velocity at this location is equal to local speed of sound.
The area of nozzle changes as a function of distance along the nozzle (x), we make the assumption that the flow properties vary only with x, to assume uniform flow properties accross any given cross section. Such flow is defined as quasi-one-Dimensional.
II. Project Objectives :
- Write matlab code to Simulate and solve 1D subsonic supersonic nozzle flow equations using MacCormacks method.
- Implementing both Conservative and Non conservative forms of governing equations and making seperate functions both forms.
- To find out minimum number of cycles for which the simulation should be run to achieve convergence.
- Compare the normalized mass flow rate between Conservative and Non-conservative forms
- Timewise variations of density, temperature, velocity and mach no at nozzle throat.
- Distribution of Mass flow through nozzle at different times during time-marching approach to steady state.
- Perform grid dependence test.
III. Governing Flow equations (Non-Dimensional variables) :
Defining the Non-Dimensional variables as,
Temperature : `T' = T/T_0`
Density : `rho' = rho/rho_0`
Dimensionless length : `x' = x/L`
Speed of sound : `a_o = sqrt (lambdaRT_0)`
Velocity : `v' = v/a_0`
Time : `t' = t/(L/a_0)`
Area : `A' = A/A^(throat area)`
Introducing these Dimensionless variables in the governing flow equations, we get :
III. A) For Non-conservative form :
Continuity equation : `(delrho')/(delt') = -rho'((delV')/(delx')) - rho'V'((del(lnA'))/(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'))]`
III. B) For Conservative form :
Continuity equation : `(del(rho'A'))/(delt') + (del(rho'A'V'))/(delx')`
Momentum equation : `(del(rho'A'V'))/(delt') + (del[rho'A'V'^2 + (1/gamma)p'A'])/(delx') = 1/gamma p'(delA')/(delx')`
Energy equation : `(del[rho'((e')/(gamma-1) + lambda/2 V'^2)A'])/(delt') + (del[rho'((e')/(gamma-1) + lambda/2 V'^2)V'A' + p'A'V'])/(delx') = 0`
Now, the generic form of equation can be expressed as :
Solution vector : U
Flux term : F
Source term : J
Defining the elements of Solution vector U
`U1 = rho'A'`
`U2 = rho'A'V'`
`U3 = rho'((e')/(gamma-1) + gamma/2 V'^2)A'`
Pure form of Flux terms :
F1 = U1
`F_2 = U_2^2/(U_1) + (gamma-1)/gamma (U_3 - gamma/2 (U_2^2/(U_1)))`
`F_3 = gamma (U_2U_3)/U_1 - (gamma(gamma-1))/2 U_2^3/U_1^2`
Source term :
`J_2 = (gamma-1)/gamma (U_3 - gamma/2 U_2^2/U_1) (del(lnA'))/(delx')`
Primitive Variables :
We need primitive variables to decode the values after calculating the Solutions vector (U1,U2,U3) which is obtained at each time step
`rho' = U_1/A'`
`V' = U_2/U_1`
`T' = e' = gamma-1 (U_3/U_1 - gamma/2 V'^2)`
`p' = rho'T'`
IV. Now the Continuity, Momentum and Energy equations for Conservative form : (respectively) :
`(delU_1)/(delt') = - (delF_1)/(delx')`
`(delU_2)/(delt') = - (delF_2)/(delx') + J2`
`(delU_3)/(delt') = - (delF_3)/(delx')`
V. Time step :
A stability constraint exists in the system. The present governing equations are nonlinear PDE's . And an exact stability criterion an only be viewed as general guidance for present nonlinear problem. However it is a Good guidance.
The courant friedrichs lowry (CFL) condition for time step of 1 D flow , in Non dimensional quantities is given as:
`Deltat = min(C (Deltax)/(a + V))`
where, C = 0.5 (assumed in the problem)
a is local speed of sound
V is local flow velocity at a point in the flow.
`Deltat`is the increment in non dimensional time and `Deltax` is increment in nondimensional space.
- For full derivation refer : Computational fluid dynamics - Basics with applications by John D anderson
The problem Specification's has been taken from Chapter 7 of this book, and the code is ran for 1400 time steps as stated in there. The results obtained are directly compared and analyzed from the program created.
VI. Matlab Programs :
Main program :
This program calls the non conservaitive and conservative forms and plots all the results and graphs.
%% Simulation of a 1D Super-sonic nozzle flow simulation using Macormack Method
% First function called is of Non-conservative form and second one is of
% Conservative form.
% Ref - (CFD - basics with applications )John D anderson chap 7.
% code by AADIL SHAIKH
% date : 17/9/2019
clear all
close all
clc
%% inputs
n = 31;
x = linspace(0,3,n);
dx = x(2) - x(1);
throat = (n-1)/2 + 1;
gamma = 1.4;
c = 0.5;
% area
a = 1+2.2*(x-1.5).^2;
%time steps
nt = 1400;
%%
tic
[mass_flow_nc,rho_nc,v_nc,t_nc,M_nc,P_nc,rho_t_nc,t_t_nc,P_t_nc,M_t_nc] = non_conservative_form(n,x,dx,throat,gamma,c,a,nt);
Nc_time = toc;
tic
[u1,u2,u3,mass_flow_c,rho,v,t,M,P,rho_t,t_t,P_t,M_t] = Conservative_form(n,x,dx,throat,gamma,c,a,nt);
C_time = toc;
%%
for y = 1:n
analytical_massflow(y) = 0.579;
end
figure(7)
plot(x,mass_flow_nc,'linewidth',1.5)
hold on
plot(x,mass_flow_c,'linewidth',1.5)
plot (x, analytical_massflow,'linewidth',1.5)
xlabel('Non dimensional distance across nozzle')
ylabel('Mass flow variations')
legend('Non conservation form','Conservation form','Exact Solution','location','best')
saveas(gcf,'mass flow compare.png')
VI. A) Non-Conservative form :
function [mass_flow_nc,rho_nc,v_nc,t_nc,M_nc,P_nc,rho_t_nc,t_t_nc,P_t_nc,M_t_nc] = non_conservative_form(n,x,dx,throat,gamma,c,a,nt)
% calculate initial profiles
rho_nc = 1-0.3146*x;
t_nc = 1-0.2314*x; % t = temp
v_nc = (0.1+1.09*x).*t_nc.^0.5;
%% outer time loop
for k = 1:nt
%updating old values
rho_old = rho_nc;
v_old = v_nc;
t_old = t_nc;
% time step
A = (t_nc.^0.5);
dt = c.*(dx ./((A)+v_nc));
dt = min(dt);
% predictor method
for j = 2:n-1
dvdx = (v_nc(j+1)-v_nc(j))/dx;
drhodx = (rho_nc(j+1)-rho_nc(j))/dx;
dtdx = (t_nc(j+1)-t_nc(j))/dx;
dlogadx = (log(a(j+1))-log(a(j)))/dx;
% Continuity equation
drhodt_p(j) = -rho_nc(j)*dvdx - rho_nc(j)*v_nc(j)*dlogadx - v_nc(j)*drhodx;
% Momentum eqn
dvdt_p(j) = -v_nc(j)*dvdx - (1/gamma)*(dtdx + (t_nc(j)/rho_nc(j))*drhodx);
% Energy eqn
dtdt_p(j) = -v_nc(j)*dtdx - (gamma - 1)*t_nc(j)*(dvdx + v_nc(j)*dlogadx);
% Soln update
rho_nc(j) = rho_nc(j) + drhodt_p(j)*dt;
v_nc(j) = v_nc(j) + dvdt_p(j)*dt;
t_nc(j) = t_nc(j) + dtdt_p(j)*dt;
end
% Corrector method
for j = 2:n-1
dvdx = (v_nc(j)-v_nc(j-1))/dx;
drhodx = (rho_nc(j)-rho_nc(j-1))/dx;
dtdx = (t_nc(j)-t_nc(j-1))/dx;
dlogadx = (log(a(j))-log(a(j-1)))/dx;
% Continuity equation
drhodt_c(j) = -rho_nc(j)*dvdx - rho_nc(j)*v_nc(j)*dlogadx - v_nc(j)*drhodx;
% Momentum eqn
dvdt_c(j) = -v_nc(j)*dvdx - (1/gamma)*(dtdx + (t_nc(j)/rho_nc(j))*drhodx);
% Energy eqn
dtdt_c(j) = -v_nc(j)*dtdx - (gamma - 1)*t_nc(j)*(dvdx + v_nc(j)*dlogadx);
end
% Compute avg time derivative
drhodt = 0.5*(drhodt_p + drhodt_c);
dvdt = 0.5*(dvdt_p + dvdt_c);
dtdt = 0.5*(dtdt_p + dtdt_c);
% Final soln update
for i = 2:n-1
rho_nc(i) = rho_old(i) + drhodt(i)*dt;
v_nc(i) = v_old(i) + dvdt(i)*dt;
t_nc(i) = t_old(i) + dtdt(i)*dt;
end
%% Apply boundary conditions
% Inlet
v_nc(1) = 2*v_nc(2) - v_nc(3);
% outlet
v_nc(n) = 2*v_nc(n-1) - v_nc(n-2);
rho_nc(n) = 2*rho_nc(n-1) - rho_nc(n-2);
t_nc(n) = 2*t_nc(n-1) - t_nc(n-2);
M_nc = v_nc./(t_nc.^0.5); % Mach no
P_nc = rho_nc.*t_nc; % pressure
mass_flow_nc = rho_nc.*v_nc.*a;
rho_t_nc(k) = rho_nc(throat);
t_t_nc(k) = t_nc(throat);
P_t_nc(k) = P_nc(throat);
M_t_nc(k) = M_nc(throat);
%%
figure(1);
hold on
if k == 50
plot(x,mass_flow_nc,'color','g','linewidth',1.5)
elseif k == 100
plot(x,mass_flow_nc,'color','m','linewidth',1.5)
elseif k == 150
plot(x,mass_flow_nc,'color','b','linewidth',1.5)
elseif k == 200
plot(x,mass_flow_nc,'color','r','linewidth',1.5)
elseif k == 700
plot(x,mass_flow_nc,'color','c','linewidth',1.5)
elseif k == 1400
plot(x,mass_flow_nc,'o','color','bla','linewidth',1.5)
xlabel('Distance through nozzle(x)-(Nondimensional)')
ylabel('Mass Flow (Nondimensional)')
grid minor
legend('50 Deltat','100 Deltat','150 Deltat','200 Deltat','700 Deltat','1400 Deltat')
title({['Non Conservative form'];['Instantaneous distributions of Mass flow through nozzle at different times during time-marching approach to steady state']})
end
end
hold on
figure(2)
X = linspace(1,nt,nt);
subplot(2,1,1)
plot(X,rho_t_nc)
xlabel('no. of time steps')
ylabel('Density')
grid on
title({['Non Conservative form'];['Timewise variations of Density & Temperature at nozzle throat']})
subplot(2,1,2)
plot(X,t_t_nc)
xlabel('no. of time steps')
ylabel('Temperature')
grid on
figure(3)
subplot(2,1,1)
plot(X,P_t_nc)
xlabel('no. of time steps')
ylabel('Pressure')
grid on
title({['Non Conservative form'];['Timewise variations of Pressure & Mach no. at nozzle throat']})
subplot(2,1,2)
plot(X,M_t_nc)
xlabel('no. of time steps')
ylabel('Mach number')
grid on
endVI.B ) Conservative Form :
function [u1,u2,u3,mass_flow_c,rho,v,t,M,P,rho_t,t_t,P_t,M_t] = Conservative_form(n,x,dx,throat,gamma,c,a,nt)
% calculating initial profiles .
for i = 1:n
if (0 <= x(i)) && (x(i) <= 0.5)
rho(i) = 1;
t(i) = 1;
elseif (0.5 <= x(i)) && (x(i) <= 1.5)
rho(i) = 1 - 0.366*(x(i) - 0.5);
t(i) = 1 - 0.167*(x(i) - 0.5);
elseif (1.5 <= x(i)) && (x(i) <= 3.5)
rho(i) = 0.634 - 0.3879*(x(i) - 1.5);
t(i) = 0.833 - 0.3507*(x(i) - 1.5);
end
end
v = 0.59 ./(rho.*a);
u1 = rho.*a;
u2 = rho.*a.*v;
u3 = rho.*((t/(gamma - 1)) + (gamma/2).*v.^2).*a;
P = rho.*t;
for k = 1:nt
%updating old values
u1_old = u1;
u2_old = u2;
u3_old = u3;
% time step
A = (t.^0.5);
dt = c.*(dx ./((A)+v));
dt = min(dt);
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);
%predictor step
for j = 2:n-1
%j2 = ((gamma-1)/gamma)*(u3(j) - ((gamma/2)*(((u2(j)).^2)./u1(j))))*((log(a(j+1))-log(a(j)))/dx);
j2(j) = (1/gamma)*rho(j).*t(j).*((a(j+1) - a(j))/dx);
% continuity eqn
du1dt_p(j) = -((f1(j+1)-f1(j))/dx);
% momentum eqn
du2dt_p(j) = -((f2(j+1) - f2(j))/dx) + j2(j);
% energy eqn
du3dt_p(j) = -((f3(j+1) - f3(j))/dx);
% Solution update
u1(j) = u1(j) + du1dt_p(j)*dt ;
u2(j) = u2(j) + du2dt_p(j)*dt;
u3(j) = u3(j) + du3dt_p(j)*dt;
end
%updating primitive variables from new u1 u2 u3 values
rho = (u1)./(a);
v = u2./u1;
t = (gamma - 1)*(((u3)./(u1)) - (gamma/2)*((v).^2));
% updating predicted flux term valies from new updated u1 u2 u3 values
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);
% Corrector step
for j = 2:n-1
%j2 = ((gamma-1)/gamma)*(u3(j) - ((gamma/2)*(((u2(j)).^2)./u1(j))))*((log(a(j))-log(a(j-1)))/dx);
j2(j) = (1/gamma)*rho(j).*t(j).*((a(j) - a(j-1))/dx);
% continuity eqn
du1dt_c(j) = -((f1(j)-f1(j-1))/dx);
% momentum eqn
du2dt_c(j) = -((f2(j) - f2(j-1))/dx) + j2(j);
% energy eqn
du3dt_c(j) = -((f3(j) - f3(j-1))/dx);
end
% Compute avg time derivative
du1dt = 0.5*(du1dt_p + du1dt_c );
du2dt = 0.5*(du2dt_p + du2dt_c );
du3dt = 0.5*(du3dt_p + du3dt_c );
% final soln 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
% Applying boundary conditions
% inflow bc
u1(1) = rho(1)*a(1);
u2(1) = 2*u2(2) - u2(3);
u3(1) = u1(1)*((t(1)./(gamma-1)) + ((gamma/2)*v(1).^2));
%outflow bc
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);
% corrected values of primitive variables
rho = (u1)./(a);
v = u2./u1;
t = (gamma - 1)*(((u3)./(u1)) - (gamma/2)*((v).^2));
% solutions
M = v./((t).^0.5); % Mach no
P = rho.*t; % pressure
mass_flow_c = rho.*v.*a;
rho_t(k) = rho(throat);
t_t(k) = t(throat);
P_t(k) = P(throat);
M_t(k) = M(throat);
figure(4)
hold on
if k == 50
plot(x,mass_flow_c,'color','g','linewidth',1.5)
elseif k == 100
plot(x,mass_flow_c,'color','m','linewidth',1.5)
elseif k == 150
plot(x,mass_flow_c,'color','b','linewidth',1.5)
elseif k == 200
plot(x,mass_flow_c,'color','r','linewidth',1.5)
elseif k == 700
plot(x,mass_flow_c,'color','c','linewidth',1.5)
elseif k == 1400
plot(x,mass_flow_c,'o','color','bla','linewidth',1.5)
xlabel('Distance through nozzle(x)-(Nondimensional)')
ylabel('Mass Flow (Nondimensional)')
grid minor
legend('50 Deltat','100 Deltat','150 Deltat','200 Deltat','700 Deltat','1400 Deltat')
title({['Conservative form'];['Instantaneous distributions of Mass flow through nozzle at different times during time-marching approach to steady state']})
end
end
hold on
figure(5)
X = linspace(1,nt,nt);
subplot(2,1,1)
plot(X,rho_t)
xlabel('no. of time steps')
ylabel('Density')
grid on
title({['Conservative form'];['Timewise variations of Density & Temperature at nozzle throat']})
subplot(2,1,2)
plot(X,t_t)
xlabel('no. of time steps')
ylabel('Temperature')
grid on
figure(6)
subplot(2,1,1)
plot(X,P_t)
xlabel('no. of time steps')
ylabel('Pressure')
grid on
title({['Conservative form'];['Timewise variations of Pressure & Mach no. at nozzle throat']})
subplot(2,1,2)
plot(X,M_t)
xlabel('no. of time steps')
ylabel('Mach number')
grid on
endThe MacCormacks technique is implemented in both the Non-conservative and Conservative functons.
VII. Results :
VII. A) Timewise variations of density, temperature, velocity and Mach no. at nozzle throat (31 grid)
1. Non-conservative form :
2. Conservative form :
Initial disturbances occur in both the form's variatons but,
We can observe that the Timewise variations in the conservative form is less as compared to the Non conservative form.
VII. B) Instantaneous distributions of Mass flow through nozzle at different times during time marching approach to steady state (31 grid)
1. Non conservative form
2. Conservative form :
The mass flow distributions through the nozzle vary alot in both forms, however in Conservative solution the variation tends towards the steady state approach abit more than the non-conservative solution. It is a point to note that the initial conditions were different for both forms.
VII. C) Comparison of mass flows between Conservative and Non-conservative forms (31 grid):
Observing the plot above, mass flows of both solutions are very close to the analytical value 0.579 as given in the book. however, The conservation form gives a distribution much closer to being a constant. In contrast, the non conservation results have(on a magnified scale) a sizeable variation, with some spurious oscillations at both inflow and outflow boundaries.
This shows that the Conservation form does a better job of preserving the mass throughout the flow field. whereas in Nonconservation form the massflow is actually a secondary result due to equations used to obtain above results.
VII. D) Grid dependence test :
We will increase the grid size from 31 to 61 to see if the solution is dependent or independent of grid size.
All the values are fairly close to the Exact analytical solution however doubling the grid points did improve the numerical solution marginally. We notice that the non-conservative form values are a tad-bit more close to analytical solution than the conservation form.
The same can be said is true for all locations within the nozzle. in otherwords, the steady state numerical solutions are essentially same and we can therefore conclude that our original 31 grid points solution is essentially grid-independent for all practical purposes.
The grid independent solution does not agree exactly with the analytical results but it is close enough and have great accuracy.
VII. E) Minimum number of cycles for which the simulation should be run to achieve convergence and its Computational time :
Inspite the the fact that conservation form requires more tedious calculations and work it surprisingly converged faster than the Non-conservation form and Minimum cycles for which the simulation should be run displayed. These values are approximated from the graphs
VII. F) Similar Plots for 61 Grid points
1. Nonconservation form
2. Conservation form :
VII. G) Instantaneous distributions of Mass flow through nozzle at different times during time marching approach to steady state (61 grid)
1. Nonconservative form :
2. Conservative form :
VIII. Conclusions :
1. There is no Superiority of either form of governing equations over the other for the given flow problem in terms of accuracy of results. The non-conservation form seems to produce slightly more accurate results for primitive variables and conservation form seems to produce sightly more accurate results for flux variables.
2. The conservation form yields a better mass flow distribution and does a better job of conserving mass.
3. The non-conservation form leads to small residuals.
4. The Conservative method came out to be faster in this case.
5. The results in either case are certainly satisfactory.
6. The solution is grid independent .
IX : References
1. Computational fluid dynamics- basics with applications by John D anderson
The End
keywords - MATLAB, MATLAB-BASICS, NUMERICAL-ANALYSIS, CFD, MACCORMACK, NOZZLE FLOW