SImulation of 1D Superosnic nozzle flow using Maccormack's method
Objectives:
- Implement both the Conservative and Non-conservative form of governing equations.
- Display the timewise variation of the field variables.
- Display the steady state variation of the flow field variables for the conservative and Non conservative form of the Quasi 1D flow.
- Comparison of mass flow rates for Conservative and Non conservative forms
- Performing a grid dependancy test.
In this problem, we study the case of isentropic flow through a convergent-divergent nozzle. The flow is subsonic at the inlet (M<1), sonic at the throat (M=1) and supersonic at the exit (M>1). We assume that the flow comes to the inlet from the reservoir. The area tends to infinity and hence the velocity is almost zero. But this is not the case, since at the inlet, the area has a finite value and hence M>0. Here, we assume the field variables vary only with x, hence its said to be quasi-one dimensional flow.
The governing eqautions are obtained in the PDE form using non-dimensionalized variables T',ρ', p' and A' ( Temperature, density, pressure and area respectively). The nozzle flow is a non-linear hyperbolic PDE problem.
Governing Equations:
1) Non Conservative form
a) Continuity equation
`(delrho)/(delt) = -rho*(delV)/(delx) - rho*V*(del(lnA))/(delx) - V*(delrho)/(delx)`
b) Momentum Equation
`(delV)/(delt)= -V*(delV)/(delx) - (1/gamma)*((delT)/(delx) + T/rho*(delrho)/(delx))`
c) Energy Equation
`(delT)/(delt)= -V*(delT)/(delx) - (gamma - 1)*T*((delV)/(delx) + V*(del(lnA))/(delx))`
2) Conservative form (In generic form)
Solution vectors:
`U1 = rho*A`
`U2 = rho*A*v`
`U3 = rho*A*((T/(gamma - 1)) + ((gamma/2)*(V^2)))`
Flux vectorss:
`F1 = U2`
`F2 = ((U2)^2)/(U1) + ((gamma-1)/gamma)*(U3 - (gamma*(U2)^2)/(2*U1))`
`F3 = ((gamma*U2*U3)/(U1)) - ((gamma*(gamma-1)/2)*((U2)^3/(U1)^2))`
Source term:
`J2 = 1/gamma*p*((delA)/(delx))`
The problem is solved using Maccormack method which is a predictor-corrector type method. The number of grid points 'n' chosen is 31. Different initial conditions for 'rho','v','T' and 'A' are assumed for the non conservative and conservative forms. A Courant number of 0.5 is chosen. The time step 'dt' is calculated based on the Courant number, speed of sound and initial velocity.
Following is the code for the main program :
Main Program
close all
clear all
clc
% Main program to simulate 1D Subsonic-Supersonic isentropic nozzle flow using Maccormacks's technique
% Non conservative and Conservative forms are used
% Input conditions
L = 3; % Length of nozzle
n = 31; % Number of grid points
x = linspace(0,L,n); % Grid
dx = x(2) - x(1); % Distance between each node
gamma = 1.4; % Specific heat ratio
nt = 1400; % Number of time steps
C = 0.5; % Assumption of Courant number
% Initial conditions (Non dimensionalized)
rho = 1 - 0.3146*x; % Density
T = 1 - 0.2314*x; % Temperature
v = (0.1 + (1.09*x)).*T.^0.5; % Velocity
A = 1 + 2.2*(x - 1.5).^2; % Area
% Throat position % A = A* at throat
n_throat = (n - 1)/2; % Throat is at 1.5m
a = sqrt(T); % Non dimensionalized Speed of sound
% Time step using CFL criteria
dt = min((C*dx)./(a + v)); % Non conservative
% --Non Conservative Solver--
tic
[mass_flow rho T p M] = Non_conservative_1D_supersonic_nozzle(n,x,dx,gamma,nt,rho,T,v,A,n_throat,a,dt);
t1 = toc;
fprintf('Computation time (Non conservative) = %0.3f \n',t1);
% --Conservative Solver--
tic
[mass_flow_c rho_c T_c p_c M_c] = Conservative_1D_supersonic_nozzle(n,x,dx,gamma,nt,C,A,n_throat);
t2 = toc;
fprintf('Computation time (Conservative) = %0.3f \n',t2);
% Steady state variations of mass flow rate of Conservative and Non conservative forms
figure(7)
plot(x,mass_flow,'g','linewidth',1.25); % Non conservative
hold on
plot(x,mass_flow_c,'r','linewidth',1.25); % Conservative
y = 0.579; % Exact analytical solution
line([0 3],[y y],'linewidth',1.25);
title('Steady state mass flow variations');
legend('Non conservative','Conservative','Exact analytical solution','Location','north')
xlabel('x/L')
ylabel('Non dimensional mass flow (\rhoVA)/(\rho_{o}a_{o}A^{*})')
grid minor
Following is the function for Non conservative form :
% Non conservative solver for 1D Subsonic-Supersonic nozzle flow using Maccormack's method
function [mass_flow rho T p M] = Non_conservative_1D_supersonic_nozzle(n,x,dx,gamma,nt,rho,T,v,A,n_throat,a,dt)
for k = 1:nt % Time loop
% Creating a copy of density, velocity and temperature
rho_old = rho;
v_old = v;
T_old = T;
for i = 2:n-1 % Loop for inner nodes
% Predictor step (using forward differencing)
dv_dx = (v(i+1) - v(i))/dx;
drho_dx = (rho(i+1) - rho(i))/dx;
dT_dx = (T(i+1) - T(i))/dx;
dlogA_dx = (log(A(i+1)) - log(A(i)))/dx;
% Continuity equation
drho_dt_p(i) = (-rho(i)*dv_dx) - (rho(i)*v(i)*dlogA_dx) - (v(i)*drho_dx);
% Momentum equation
dv_dt_p(i) = (-v(i)*dv_dx) - ((1/gamma)*(dT_dx + (T(i)/rho(i))*drho_dx));
% Energy equation
dT_dt_p(i) = (-v(i)*dT_dx) - ((gamma - 1)*T(i)*(dv_dx + v(i)*dlogA_dx));
% First solution update
rho(i) = rho(i) + drho_dt_p(i)*dt;
v(i) = v(i) + dv_dt_p(i)*dt;
T(i) = T(i) + dT_dt_p(i)*dt;
end
% Corrector step (using backward differencing)
for i = 2:n-1 % Loop for inner nodes
dv_dx = (v(i) - v(i-1))/dx;
drho_dx = (rho(i) - rho(i-1))/dx;
dT_dx = (T(i) - T(i-1))/dx;
dlogA_dx = (log(A(i)) - log(A(i-1)))/dx;
% Continuity equation
drho_dt_c(i) = (-rho(i)*dv_dx) - (rho(i)*v(i)*dlogA_dx) - (v(i)*drho_dx);
% Momentum equation
dv_dt_c(i) = (-v(i)*dv_dx) - ((1/gamma)*(dT_dx + (T(i)/rho(i))*drho_dx));
% Energy equation
dT_dt_c(i) = (-v(i)*dT_dx) - ((gamma - 1)*T(i)*(dv_dx + v(i)*dlogA_dx));
end
% Computing Average time derivative
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 solution update
for j = 2:n-1
rho(j) = rho_old(j) + drho_dt_avg(j)*dt;
v(j) = v_old(j) + dv_dt_avg(j)*dt;
T(j) = T_old(j) + dT_dt_avg(j)*dt;
end
% Applying Boundary conditions
% Inlet
rho(1) = 1;
v(1) = 2*v(2) - v(3);
T(1) = 1;
% Outlet
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 = v./sqrt(T); % Mach number
mass_flow = rho.*A.*v; % Mass flow rate (Non dimensional)
p = rho.*T; % Pressure (Non dimensional)
% Nozzle throat flow field variables
p_throat(k) = p(n_throat); % Pressure
rho_throat(k) = rho(n_throat); % Density
T_throat(k) = T(n_throat); % Temperature
M_throat(k) = M(n_throat); % Mach number
% Instantaneous distribution of mass flow rate as a function of distance at different time steps
figure(3)
if (k == 1)
plot(x,mass_flow,'--b')
hold on
elseif (k == 50)
plot(x,mass_flow,'g','linewidth',1.25)
elseif (k == 100)
plot(x,mass_flow,'m','linewidth',1.25)
elseif (k == 150)
plot(x,mass_flow,'r','linewidth',1.25)
elseif (k == 200)
plot(x,mass_flow,'c','linewidth',1.25)
elseif (k == 700)
plot(x,mass_flow,'k','linewidth',1.25)
title('Instantaneous distribution of mass flow rate as a function of distance at different time steps (Non conservative')
xlabel('Non dimensional distance(x)');
ylabel('Non dimensional mass flow (\rhoVA)/(\rho_{o}a_{o}A^{*})');
legend('0\Deltat','50\Deltat','100\Deltat','150\Deltat','200\Deltat','700\Deltat');
grid minor
end
end
% Plotting results for conditions at nozzle throat
figure(1)
subplot(4,1,1);
suptitle('Timewise variation of flow field at nozzle throat - Time step vs Field variables (Non conservative)')
plot(p_throat,'r')
title('Pressure');
ylabel('p/p_{o}');
axis([0 1400 0.4 1.2])
grid minor
subplot(4,1,2);
plot(rho_throat,'b')
title('Density');
ylabel('\rho/\rho_{o}');
axis([0 1400 0.6 1.1])
grid minor
subplot(4,1,3);
plot(T_throat,'g')
title('Temperature');
ylabel('T/T_{o}');
axis([0 1400 0.6 1.1])
grid minor
subplot(4,1,4);
plot(M_throat,'m')
title('Mach number');
xlabel('Number of time steps(k)');
ylabel('Mach number');
axis([0 1400 0.5 1.5])
grid minor
% Plotting steady state values
figure(2)
subplot(4,1,1);
suptitle('Steady state profile of properties of Quasi-1D flow (Non conservative)')
plot(x,p,'r')
ylabel('p/p_{o}');
axis([0 3 0 1.2]);
legend('Pressure ratio')
grid minor
subplot(4,1,2);
plot(x,rho,'b')
ylabel('\rho/\rho_{o}');
axis([0 3 0 1.2]);
legend('Density ratio')
grid minor
subplot(4,1,3);
plot(x,T,'g')
ylabel('T/T_{o}');
axis([0 3 0 1.2]);
legend('Temperature ratio')
grid minor
subplot(4,1,4);
plot(x,M,'m')
ylabel('Mach number');
legend('Mach number','Location','best')
axis([0 3 0 4]);
grid minorFollowing is the function for Conservative form :
% Conservative solver for 1D Subsonic-Supersonic nozzle flow using Maccormack's method
function [mass_flow_c rho_c T_c p_c M_c] = Conservative_1D_supersonic_nozzle(n,x,dx,gamma,nt,C,A,n_throat)
% Initial conditions
for i = 1:length(x)
if (x(i) >= 0 && x(i) < 0.5)
rho_c(i) = 1;
T_c(i) = 1;
elseif (x(i) >= 0.5 && x(i) < 1.5)
rho_c(i) = 1 - 0.366*(x(i) - 0.5);
T_c(i) = 1 - 0.167*(x(i) - 0.5);
elseif (x(i) >= 1.5 && x(i) <= 3)
rho_c(i) = 0.634 - 0.3879*(x(i) - 1.5);
T_c(i) = 0.833 - 0.3507*(x(i) - 1.5);
end
end
rho_c_old = rho_c;
v_c = 0.59./(rho_c.*A); % Analytical mass flow rate = 0.59
p_c = rho_c.*T_c; % Non dimenisonal Pressure
mass_flow_c_initial = rho_c.*A.*v_c;
% Initializing solution vectors
U1 = rho_c.*A;
U2 = rho_c.*A.*v_c;
U3 = rho_c.*A.*((T_c/(gamma - 1)) + ((gamma/2)*(v_c.^2)));
a_c = sqrt(T_c); % Non dimensionalized speed of sound
dt_c = min((C.*dx)./(a_c + v_c)); % Time step using CFL criteria
for k = 1:nt % time loop
% Creating a copy of solution vectors
U1_old = U1;
U2_old = U2;
U3_old = U3;
% Defining flux terms
F1 = U2;
F2 = (U2.^2./U1) + ((gamma-1)/gamma)*(U3 - (gamma*U2.^2)./(2*U1));
F3 = ((gamma*U2.*U3)./U1) - ((gamma*(gamma-1)/2)*(U2.^3./U1.^2));
for i = 2:n-1 % Loop for inner nodes
J2(i) = (1/gamma).*rho_c(i).*T_c(i).*((A(i+1)-A(i))/dx);
dF1_dx = (F1(i+1) - F1(i))/dx;
dF2_dx = (F2(i+1) - F2(i))/dx;
dF3_dx = (F3(i+1) - F3(i))/dx;
% Predictor step
dU1_dt_p(i) = - dF1_dx;
dU2_dt_p(i) = - dF2_dx + J2(i);
dU3_dt_p(i) = - dF3_dx;
% First solution update
U1(i) = U1(i) + dU1_dt_p(i)*dt_c;
U2(i) = U2(i) + dU2_dt_p(i)*dt_c;
U3(i) = U3(i) + dU3_dt_p(i)*dt_c;
end
% Updating flux and primitive variables
F1 = U2;
F2 = (U2.^2./U1) + ((gamma-1)/gamma)*(U3 - (gamma*U2.^2)./(2*U1));
F3 = ((gamma*U2.*U3)./U1) - ((gamma*(gamma-1)/2)*(U2.^3./U1.^2));
rho_c = U1./A;
v_c = U2./U1;
T_c = (gamma -1)*(U3./U1 - (gamma/2)*v_c.^2);
for i = 2:n-1 % Corrector step
J2(i) = (1/gamma).*rho_c(i).*T_c(i).*((A(i)-A(i-1))/dx);
dF1_dx = (F1(i) - F1(i-1))/dx;
dF2_dx = (F2(i) - F2(i-1))/dx;
dF3_dx = (F3(i) - F3(i-1))/dx;
dU1_dt_c(i) = - dF1_dx;
dU2_dt_c(i) = - dF2_dx + J2(i);
dU3_dt_c(i) = - dF3_dx;
end
% Computing 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 solution update
for j = 2:n-1
U1(j) = U1_old(j) + dU1_dt_avg(j)*dt_c;
U2(j) = U2_old(j) + dU2_dt_avg(j)*dt_c;
U3(j) = U3_old(j) + dU3_dt_avg(j)*dt_c;
end
% Applying boundary conditions
% Inlet
U1(1) = rho_c(1)*A(1); % rho = 1 at i = 1
U2(1) = 2*U2(2) - U2(3);
U3(1) = U1(1).*(((T_c(1)/(gamma-1)) + ((gamma/2)*v_c(1).^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 primitive variables
rho_c = U1./A;
v_c = U2./U1;
T_c = (gamma - 1)*((U3./U1) - ((gamma/2)*v_c.^2));
p_c = rho_c.*T_c;
M_c = v_c./sqrt(T_c); % Mach number
mass_flow_c = rho_c.*A.*v_c; % Mass flow rate (Non dimensional)
% Nozzle throat flow field variables
p_throat_c(k) = p_c(n_throat); % Pressure
rho_throat_c(k) = rho_c(n_throat); % Density
T_throat_c(k) = T_c(n_throat); % Temperature
M_throat_c(k) = M_c(n_throat); % Mach number
figure(4)
if (k == 1)
plot(x,mass_flow_c_initial,'--b')
hold on
elseif (k == 100)
plot(x,mass_flow_c,'m','linewidth',1.25)
elseif (k == 200)
plot(x,mass_flow_c,'c','linewidth',1.25)
elseif (k == 700)
plot(x,mass_flow_c,'k','linewidth',1.25)
title('Instantaneous distribution of mass flow rate as a function of distance at different time steps (Conservative)')
xlabel('Non dimensional distance(x)');
ylabel('Non dimensional mass flow (\rhoVA)/(\rho_{o}a_{o}A^{*})');
legend('0\Deltat','100\Deltat','200\Deltat','700\Deltat');
ylim([0.4 0.7]);
grid minor
end
end
figure(5)
subplot(4,1,1);
suptitle('Steady state profile of properties of Quasi-1D flow (Conservative)')
plot(x,p_c,'r')
ylabel('p/p_{o}');
axis([0 3 0 1.2]);
legend('Pressure ratio')
grid minor
subplot(4,1,2);
plot(x,rho_c,'b')
ylabel('\rho/\rho_{o}');
axis([0 3 0 1.2]);
legend('Density ratio')
grid minor
subplot(4,1,3);
plot(x,T_c,'g')
ylabel('T/T_{o}');
axis([0 3 0 1.2]);
legend('Temperature ratio')
grid minor
subplot(4,1,4);
plot(x,M_c,'m')
ylabel('Mach number');
legend('Mach number','Location','best')
axis([0 3 0 4]);
grid minor
figure(6)
subplot(4,1,1);
suptitle('Timewise variation of flow field at nozzle throat - Time step vs Field variables (Conservative)')
plot(p_throat_c,'r')
title('Pressure');
ylabel('p/p_{o}');
axis([0 1400 0.4 1.2])
grid minor
subplot(4,1,2);
plot(rho_throat_c,'b')
title('Density');
ylabel('\rho/\rho_{o}');
axis([0 1400 0.6 1.1])
grid minor
subplot(4,1,3);
plot(T_throat_c,'g')
title('Temperature');
ylabel('T/T_{o}');
axis([0 1400 0.6 1.1])
grid minor
subplot(4,1,4);
plot(M_throat_c,'m')
title('Mach number');
xlabel('Number of time steps(k)');
ylabel('Mach number');
axis([0 1400 0.5 1.5])
grid minorResults:
1) Timewise variation of flow field variables at nozzle throat
i) Non conservative
ii) Conservative
The pressure, density, temperature and Mach number at the nozzle throat are plotted as a function of number of time steps. We see that there are oscillations that take place initially in the Non conservative form which dies out as the number of time steps increases. This is called as the 'rubber band effect' , where initially it is stretched and as steady state is reached, the field variables are less stretched or 'relaxed'. Whereas in the Conservative form, the oscillations are of small amplitude and die down quickly and becomes stable.
2) Minimum number of iterations to converge:
We see that the minimum number of iterations to converge for the Conserative form is more. This is because we are not directly solving the primtive variables. Rather the solution vectors are solved after which the flux vectors are obtained which are then used to obtain primitve variables. In the Conservative form, the flux vector is the dependant quantity. Whereas in the Non conservative form, the primitive variables are the dependant quantity.
3) Steady state values of properties
i) Non conservative
ii) Conservative
From both the steady state plots, we see that for the pressure, density and temperature drops at the outlet at the nozzle as expected. The Mach number increases as the flow reaches the outlet as it reaches supersonic speed.
4) Distribution of Mass flow rates as a function of time steps (Non conservative)
In this graph, the non dimensionalized mass flow rate is plotted as a function of distance along the nozzle. The dotted line represents the initial mass flow rate at '0deltat'. At '50deltat', the mass flow rate has changed considerably. From '100deltat' to '200deltat', we see that the mass flow keeps varying in the nozzle due to transient variations. Later, it starts becoming stable. At '700deltat' the mass flow has converged to a steady state value. It is close to the analytical result of 0.579 (mass flow rate).
5) Distribution of Mass flow rates as a function of time steps (Conservative)
Compared to the previous plot for Non conservative, we see that there is less transient variations in the mass flow at different time steps and converges in a smoother way. The variations in mass flow is much less severe in the Conservative approach.
6) Steady state mass flow variations (Conservative and Non conservative)
Here the steady state mass flow predicted by the Conservative form is more stable compared to the Non conservative form. The Non conservative has oscillations at the inlet and outlet of the nozzle. The steady state mass flow obtained by Conservative form is more closer to the analytical solution (0.579) compared to the Non conservative form. The mass flow is well conserved in the Conservative form as it is the dependant variable (along with the other flux variables) whereas the primtive variables are the dependant quantity in the Non conservative form.
7) Grid Dependance test
Here the effect of Grid dependance is employed to the problem. We see in both the Conservative and Non conservative form, when we increase the number of grid points to 61, the value of primitive variables are closer to the analytical solution. But the values of primitive variables at 31 grid points is also close to the values at 'n=61'. Hence we can say that at 'n=31', we achieved grid independance.
We see that the values of Non conservative form is closer to the analytical solution. By increasing the number of grid points to 61 for the Conservative form, the primitve variable values are closer to the analytical solution but less accurate compared to the Non conservative form with 31 grid points. The residuals are smaller in case of the Non conservative approach.
In conclusion, the superiority of one method over the other depends on the conservation,accuracy, residuals and computational efficiency. But we cannot establish the superiority of Conservative or Non conservative in this given case.