Simulation of a 1D Super-sonic Nozzle Flow using Macormack Method
Objective: -
In this challenge, we are going to write a code solve the Quasi 1D supersonic nozzle flow equations using the Macormack Method implementing both the conservative and non-conservative forms of the governing equations and also perform a grid dependence test.
Solution: -
For the current challenge, we are going to write a code to solve Quasi 1D supersonic nozzle flow implementing both conservative and non-conservative forms using the McCormack method.
The above figure shows the generic Isentropic nozzle flow schematic which we are going to assume for this simulation (Figure taken from CFD Basics and Applications by J. D. Anderson).
MacCormack method is a widely used discretization scheme for the numerical solution of hyperbolic partial differential equations. It is an explcit scheme which is second order accurate and comprises of two steps namely predictor and corrector which use forward and backward differencing schemes respectively.
The formulae for the Conservative and the Non-Conservative forms are given below:
Conservative Form:
Continuity Equation:
`(delU_1)/(delt) = -(delF_1)/(delx)`
Momentum Equation:
`(delU_2)/(delt) = -(delF_2)/(delx) + J_2`
Energy Equation:
`(delU_3)/(delt) = -(delF_3)/(delx)`
where
`U_1 = rho*A,`
`U_2 = rho*A*v and`
`U_3 = rho*A*((T)/(gamma-1) + ((gamma)/(2)*v^2))` are the Solution vectors
`F_1 = rho*A*v,`
`F_2 = rho*A*v + (1/gamma)*rho*A and`
`F_3 = rho*(e/(gamma-1)+(gamma/2)*v^2)*v*A +rho*v*A` are the flux vectors
`J_2 = (1/gamma)*rho*((delA)/(delx))` is known as the source term
Non-Conservative Form:
Continuity Equation:
`(delrho)/(delt) = -rho*(delv)/(delx) -rho*v*(delln(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*(delln(A))/(delx))`
Note:
All of the above equations are non-dimensional.
Code:
Main Program
clear all
close all
clc
% Inputs
n = 31;
nt = 1500;
x = linspace(0,3,n);
dx = x(2) - x(1);
gamma = 1.4;
c = 0.5;
throat = (n-1)/2;
a = 1 + 2.2*(x - 1.5).^2; % Area
[ sim_time_c, th_mach_c, th_press_c, th_temp_c, th_rho_c, mach_no_c, t_c, rho_c, p_c,m_dot_c ] = conservative( n, nt, x, dx, c, a, gamma, throat )
[ sim_time_nc, th_mach_nc, th_press_nc, th_temp_nc, th_rho_nc, mach_no_nc, t_nc, rho_nc, p_nc,m_dot_nc ] = non_conservative( n, nt, x, dx, c, a, gamma, throat )
fprintf('Simulation time for Conservative Form is %0.3g seconds', sim_time_c)
fprintf('\nSimulation time for Non-Conservative Form is %0.3g seconds', sim_time_nc)
% Plots
% Timewise Variation at the Nozzle Throat
figure(3)
subplot(411)
hold on
plot(th_mach_c, 'b')
plot(th_mach_nc, 'r')
legend('Conservative Form', 'Non-Conservative Form');
ylabel('Mach Number')
title('Timewise Variation at the Nozzle Throat')
grid minor
subplot(412)
hold on
plot(th_press_c, 'b')
plot(th_press_nc, 'r')
legend('Conservative Form', 'Non-Conservative Form');
ylabel('Pressure Ratio')
grid minor
subplot(413)
hold on
plot(th_rho_c, 'b')
plot(th_rho_nc, 'r')
legend('Conservative Form', 'Non-Conservative Form');
ylabel('Density Ratio')
grid minor
subplot(414)
hold on
plot(th_temp_c, 'b')
plot(th_temp_nc, 'r')
legend('Conservative Form', 'Non-Conservative Form');
xlabel('Number of Iterations')
ylabel('Temperature Ratio')
grid minor
% Qualitative Aspects of Quasi 1-D Nozzle Flow
figure(4)
subplot(411)
hold on
plot(x, mach_no_c, '-b+')
plot(x, mach_no_nc, 'r')
leg5 = legend('Conservative Form', 'Non-Conservative Form');
set(leg5, 'Location', 'northeastoutside')
ylabel('Mach Number')
title('Qualitative Aspects of Quasi 1-D Nozzle Flow')
grid minor
subplot(412)
hold on
plot(x, p_c, '-b+')
plot(x, p_nc, 'r')
leg6 = legend('Conservative Form', 'Non-Conservative Form');
set(leg6, 'Location', 'northeastoutside')
ylabel('Pressure Ratio')
grid minor
subplot(413)
hold on
plot(x, rho_c, '-b+')
plot(x, rho_nc, 'r')
leg7 = legend('Conservative Form', 'Non-Conservative Form');
set(leg7, 'Location', 'northeastoutside')
ylabel('Density Ratio')
grid minor
subplot(414)
hold on
plot(x, t_c, '-b+')
plot(x, t_nc, 'r')
leg8 = legend('Conservative Form', 'Non-Conservative Form');
set(leg8, 'Location', 'northeastoutside')
xlabel('Non-Dimensional Length of Nozzle')
ylabel('Temperature Ratio')
grid minor
figure(5)
hold on
plot(x, m_dot_c, 'b')
plot(x, m_dot_nc, 'r')
line([0 3], [0.579 0.579], 'color', 'g')
ylim([0.575 0.597])
leg8 = legend('Conservative Form', 'Non-Conservative Form', 'Exact Solution');
set(leg8, 'Location', 'northeastoutside')
title('Non-Dimensional Mass Flow Distribution');
xlabel('Non-Dimensional Length of Nozzle (x/l)')
ylabel('Mass Flow Rate Ratio')
grid minorFunction for Conservative Form
function [ sim_time_c, th_mach_c, th_press_c, th_temp_c, th_rho_c, mach_no_c, t_c, rho_c, p_c, m_dot_c ] = conservative( n, nt, x, dx, c, a, gamma, throat )
% Conservative Form
% Calculating Initial Profiles
for i = 1:n
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.5)
rho_c(i) = 0.634 - 0.3879*(x(i) - 1.5);
t_c(i) = 0.833 - 0.3507*(x(i) - 1.5);
end
end
v_c = 0.59./(rho_c.*a);
p_c = rho_c.*t_c;
% time steps
dt = min(c*dx./(sqrt(t_c) + v_c));
% Defining 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));
% Outer Time Loop
tic
for k = 1:nt
% Saving old values
rho_old = rho_c;
v_old = v_c;
t_old = t_c;
U1_old = U1;
U2_old = U2;
U3_old = U3;
% Flux vectors
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);
ctr = 1;
% Predictor
for i = 2:n-1
% Source term
J2 = (1/gamma)*rho_c(i).*t_c(i).*((a(i + 1) - a(i))/dx);
dU1dt_p(i) = -(F1(i+1) - F1(i))/dx;
dU2dt_p(i) = (-(F2(i+1) - F2(i))/dx) + J2;
dU3dt_p(i) = -(F3(i+1) - F3(i))/dx;
%Solution 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;
rho_c(i) = U1(i)./a(i);
v_c(i) = U2(i)./U1(i);
t_c(i) = (gamma - 1)*((U3(i)./U1(i)) - (gamma/2)*(v_c(i)^2));
end
for i = 2:n-1
F1(i) = U2(i);
F2(i) = (U2(i).^2./U1(i)) + ((gamma - 1)/gamma)*(U3(i) - (gamma/2)*((U2(i).^2)./U1(i)));
F3(i) = gamma*(U2(i).*U3(i)./U1(i)) - ((gamma*(gamma - 1))/2)*(U2(i).^3./U1(i).^2);
ctr = ctr + 1;
end
% Corrector
for i = 2:n-1
% Source term
J2 = (1/gamma).*rho_c(i).*t_c(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;
dU3dt_c(i) = -(F3(i) - F3(i-1))/dx;
end
% 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 Soution 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
% Apply Boundary Conditions
% Inlet
U2(1) = 2*U2(2) - U2(3);
%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);
% Flow field variables update
rho_c = U1./a;
v_c = U2./U1;
t_c = (gamma - 1)*((U3./U1) - (gamma/2)*v_c.^2);
% Calculating Mass Flow Rate
p_c = rho_old.*t_old;
mach_no_c = v_old./(t_old).^0.5;
m_dot_c = rho_old.*a.*v_old;
% Calculating values at throat
th_mach_c(k) = mach_no_c(throat);
th_press_c(k) = p_c(throat);
th_rho_c(k) = rho_c(throat);
th_temp_c(k) = t_c(throat);
% Plotting Mass Flow at different time steps
if n == 31
mass_flow_c = figure(1);
hold on
if k == 100
plot(x, m_dot_c, 'r')
elseif k == 200
plot(x, m_dot_c, 'b')
elseif k == 300
plot(x, m_dot_c, 'g')
elseif k == 550
plot(x, m_dot_c, 'm')
title('Non-Dimensional Mass Flow Distribution at Different Time Steps for Conservative Form')
legend('100 Deltat', '200 Deltat', '300 Deltat', '550 Deltat')
ylim([0.5 0.7])
xlabel('Non-Dimensional Length of Nozzle (x/l)')
ylabel('Mass Flow Rate Ratio')
grid minor
end
end
end
sim_time_c = toc
end
Function for Non-Conservative Form
function [ sim_time_nc, th_mach_nc, th_press_nc, th_temp_nc, th_rho_nc, mach_no_nc, t_nc, rho_nc, p_nc, m_dot_nc ] = non_conservative( n, nt, x, dx, c, a, gamma, throat )
% Non-Conservative Form
% Calculate Initial profiles
rho_nc = 1 - 0.3146*x; % Density
t_nc = 1 - 0.2314*x; % Temperature
v = (0.1 + 1.09*x).*t_nc.^0.5; % Velocity
% time steps
dt = min(c*dx./(t_nc.^0.5+v));
% Outer Time Loop
tic
for k = 1:nt
rho_old = rho_nc;
v_old = v;
t_old = t_nc;
% Predictor
for j = 2:n-1
dvdx = (v(j+1)-v(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 eq
drhodt_p(j) = -rho_nc(j)*dvdx - rho_nc(j)*v(j)*dlogadx - v(j)*drhodx;
% Momentum eq
dvdt_p(j) = -v(j)*dvdx - (1/gamma)*(dtdx + (t_nc(j)/rho_nc(j)*drhodx));
% Energy eq
dtdt_p(j) = -v(j)*dtdx - (gamma-1)*t_nc(j)*(dvdx + v(j)*dlogadx)
% Solution Update
v(j) = v(j) + dvdt_p(j)*dt;
t_nc(j) = t_nc(j) + dtdt_p(j)*dt;
rho_nc(j) = rho_nc(j) + drhodt_p(j)*dt;
end
% Corrector
for j = 2:n-1
dvdx = (v(j)-v(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 eq
drhodt_c(j) = -rho_nc(j)*dvdx - rho_nc(j)*v(j)*dlogadx - v(j)*drhodx;
% Momentum eq
dvdt_c(j) = -v(j)*dvdx - (1/gamma)*(dtdx + (t_nc(j)/rho_nc(j)*drhodx));
% Energy eq
dtdt_c(j) = -v(j)*dtdx - (gamma-1)*t_nc(j)*(dvdx + v(j)*dlogadx);
end
% Average Time Derivative
dvdt = 0.5*(dvdt_p + dvdt_c);
drhodt = 0.5*(drhodt_p + drhodt_c);
dtdt = 0.5*(dtdt_p + dtdt_c);
% Final Solution Update
for i = 2:n-1
v(i) = v_old(i) + dvdt(i)*dt;
rho_nc(i) = rho_old(i) + drhodt(i)*dt;
t_nc(i) = t_old(i) + dtdt(i)*dt;
end
% Apply boundary conditions
% Inlet
v(1) = 2*v(2) - v(3);
% Outlet
v(n) = 2*v(n-1) - v(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);
% Calculating Mass Flow Rate
p_nc = rho_old.*t_old;
mach_no_nc = v_old./(t_old).^0.5;
m_dot_nc = rho_old.*a.*v_old;
% Calculating values at throat
th_mach_nc(k) = mach_no_nc(throat);
th_press_nc(k) = p_nc(throat);
th_rho_nc(k) = rho_nc(throat);
th_temp_nc(k) = t_nc(throat);
% Plotting Mass Flow at different time steps
if n == 31
mass_flow_nc = figure(2);
hold on
if k == 50
plot(x, m_dot_nc, 'r')
elseif k == 200
plot(x, m_dot_nc, 'b')
elseif k == 300
plot(x, m_dot_nc, 'g')
elseif k == 500
plot(x, m_dot_nc, 'm')
title('Non-Dimensional Mass Flow Distribution at Different Time Steps for Non-Conservative Form')
legend('50 Deltat', '200 Deltat', '300 Deltat', '500 Deltat')
xlabel('Non-Dimensional Length of Nozzle (x/l)')
ylabel('Mass Flow Rate Ratio')
grid minor
end
end
end
sim_time_nc = toc
end
Results and Observation: -
Non-Dimensional Length of Nozzle for 30 Points
Non-Dimensional Length of Nozzle for 60 Points
From the figures 4 and 5, we can say that the conservation form is much closer to the exact solution when compared to the Non-conservation form. The non-conservation form has a lot of variations consisting of discrepancies or irregular oscillations at both tthe inflow and the outflow boundaries.
Comparison of Steady state results: Conservation vs Non-Conservation Form
Inference: -
The conservation forms shows a better mass flow distribution showing that it's mass conservation is better than the non-conservative form.
The non-conservation form has slight irregularities which depict the quality of the algorithm used. In this census, the non-conservation form does a better job of conserving mass.
It is very important to understand that both conservation and non-conservation forms are the same. Although the non-conservation form produces accurate results for the primitive variables, the conservation form produces accurate results for the flux variables. In either way, both of the results satisfy the case under scrutiny.
As we have established the similarity of the forms, it is also important to note that the non-conservation form has an advantage on the basis of calculational effort as the conservation form requires a higher amount of work to solve the primitive variables from the flux variables.