Simulation of a 1D Super-sonic nozzle flow simulation using Macormack Method
In this challenge, the objective is to:
- To write the MATLAB code to solve the 1D flow across a converging diverging nozzle
- Solve the equation by non-conservative and conservative method separately
- Find out the number of cycles for convergence
- Compare the mass flow rate obtained by conservative and non-conservative method
- Determine which method is faster
We consider the steady, isentropic flow through a convergent-divergent nozzle. The flow at the inlet to the nozzle comes from a reservoir where pressure and temperatures are at stagnation. Also, it is assumed that the cross-sectional area of the nozzle is very large so that the velocity is almost equal to 0. The flow expands adiabatically to supersonic speeds at the nozzle exit. the flow is locally subsonic in the convergent section (M<1), sonic at the throat section (M = 1) and supersonic at the divergent section (M>1).
We define the necessary non-dimensional variables T', v', x', A, t', ρ' as:
Temperature: T' = T/T0
Density: ρ'=ρ/ρ0
Non-dimensional Distance: x' = x/L
Non-dimensional Velocity: v' = v/a0 where a0 is the sonic velocity (a0=√γRT0)
Non-dimensional Time: t' = t/La0 where La0 is the characteristic time-scale
Non-dimensional Area: A'=A/A*
The equations must be written in their generic form in order to be solvable by numerical techniques. To construct the generic form of the governing equations we define:
Solution Vector: U
Flux Vector: F
Source Terms: J
The equations (continuity, momentum and energy) are obtained from the book by John D Anderson and then converted into in these terms:
Solution Vectors:
U1=ρ'A'
U2=ρ'A'V'
U3=ρ'(T'/(γ−1)+(γ/2)V'2)A'
Flux Vectors:
F1=ρ'A'V'
F2=ρ'A'V'2+(1/γ)p'A'
F3=ρ'(T'/(γ−1)+(γ/2)V'2)V'A'+p'A'V'
Source Terms:
J2=(1/γ)p'(∂A'/∂x')
However if the Flux terms are written in terms of the primitive variables (ρ′, T', V' ) instabilities develop in the program and the solution becomes inconsistent. To avoid this error we would be defining the flux vectors F and the source terms J in terms of the solutions vectors U. Hence,
F1=U2
F2=((U2)2/U1)+ (γ−1/γ)(U3 − (γ/2)((U2)2/U1))
F3=(γ*U2*U3/U1)− (γ*(γ−1)/2)(U2)3/(U1)2
J2=((γ−1)/γ)(U3 − (γ/2)(U2)2/U1)*(∂(lnA')/∂x')
Decomposition to the primitive variables (ρ′, T', V') we have:
Density: ρ=U1/A'
Velocity: V'=U2/U1
Temperature: T'=(γ−1)*((U3/U1)−(γ/2)V'2)
Pressure: P'=ρ'T'
MAIN CODE:
n=31;
x=linspace(0,3,n);
m=0.579*ones(n);
tic;
% solving equations in the non-conservation form
[rho_nc,T_nc,p_nc,M_nc,m_nc,n_nc] = non_conservative_form;
time_nc=toc;
tic
% solving equations in the consevtion form
[rho_c,T_c,p_c,M_c,m_c,n_c] = conservative_form;
time_c=toc;
% Comparing the normalized mass flow rate between conservation & non-conservation forms
figure(7)
plot(x,m_nc,x,m_c,x,m,'g');
xlabel('Length');
ylabel('Mass flow rate');
title(sprintf('Comparison of mass flow variations obtained with conservation and non-conservation form'));
legend('Non-conservation form','Conservation form','Analytical solution');
% comparing the convergence rate
figure(8)
c = categorical({'Non-conservation form','Conservation form'});
b = bar(c, [n_nc n_c]);
b.FaceColor = 'flat';
b.CData(2,:) = [1 1 0];
ylabel('No. of iterations');
title(sprintf('Number of iterations taken to converge to a steady-state for a tolerance of 1e-4'));
% displaying the execution times on the output window
fprintf('Time taken to execute non-conservation form = %0.5g s', time_nc);
fprintf('\nTime taken to execute conservation form = %0.5g s', time_nc);
NON CONSERVATIVE METHOD:
function [rho,t,p,M,m,count]= non_conservative_form
% initialization
n=31;
nt=1400;
x = linspace(0,3,n);
dx = x(2)-x(1);
A = 1 + 2.2*(x-1.5).^2;
throat = (n+1)/2; % assuming equal no of grid points on either side of throat
nt_x = linspace(1,nt,nt);
gamma = 1.4;
error = 1;
count = 1;
m_old = 0;
C=0.5;
tol=1e-4;
% Initial conditions
rho = 1 - 0.3146*x;
t = 1 - 0.2314*x;
v = (0.1 + 1.09.*x).*t.^(1/2);
m= 0.579*ones(n);
m_initial=rho.*A.*v;
% loop to calculate dt
for k = 1:n
dt_x(k) = C*dx/(t(k)^(1/2)+v(k));
end
dt = min(dt_x);
%outer time loop
for j = 1:nt
% updating old values
rho_old = rho;
v_old = v;
t_old = t;
% predictor step
for i = 2:n-1
drho_dx = (rho(i+1) - rho(i))/dx;
dv_dx = (v(i+1) - v(i))/dx;
dt_dx = (t(i+1) - t(i))/dx;
dlogA_dx = (log(A(i+1)) - log(A(i)))/dx;
drho_dtp(i) = -rho(i)*dv_dx - rho(i)*v(i)*dlogA_dx - v(i)*drho_dx; % continuity equation
dv_dtp(i) = -v(i)*dv_dx - (1/gamma)*(dt_dx+(t(i)/rho(i))*drho_dx); % momentum equation
dt_dtp(i) = -v(i)*dt_dx - (gamma-1)*t(i)*(dv_dx + v(i)*dlogA_dx);% energy equation
% predictor step update
rho(i) = rho_old(i) + drho_dtp(i)*dt;
v(i) = v_old(i) + dv_dtp(i)*dt;
t(i) = t_old(i) + dt_dtp(i)*dt;
end
% corrector step
for i = 2:n-1
drho_dx = (rho(i) - rho(i-1))/dx;
dv_dx = (v(i) - v(i-1))/dx;
dt_dx = (t(i) - t(i-1))/dx;
dlogA_dx = (log(A(i))-log(A(i-1)))/dx;
drho_dtc(i) = - rho(i)*dv_dx - rho(i)*v(i)*dlogA_dx - v(i)*drho_dx; % continuity equation
dv_dtc(i) = - v(i)*dv_dx - (1/gamma)*(dt_dx + (t(i)/rho(i))*drho_dx); % momentum equation
dt_dtc(i) = - v(i)*dt_dx - (gamma-1)*t(i)*(dv_dx + v(i)*dlogA_dx); % energy equation
end
% average value of time derivatives
drho_dt = 0.5*(drho_dtp + drho_dtc);
dv_dt = 0.5*(dv_dtp + dv_dtc);
dt_dt = 0.5*(dt_dtp + dt_dtc);
% corrector step update
for i = 2:n-1
rho(i) = rho_old(i) + drho_dt(i)*dt;
v(i) = v_old(i) + dv_dt(i)*dt;
t(i) = t_old(i) + dt_dt(i)*dt;
end
% Inlet boundary condition
v(1) = 2*v(2) - v(3);
% Outlet boundary conditions
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);
% Calculating pressure, mass flow rate and Mach number
p = rho.*t;
M = v./t.^0.5;
m = rho.*v.*A;
m_t(j,:) = rho.*v.*A;
% Capturing flow field variables at the throat
rho_throat(j) = rho(throat);
T_throat(j) = t(throat);
p_throat(j) = p(throat);
M_throat(j) = M(throat);
% convergence
error = max(abs(m-m_old));
if error>tol
count = count + 1;
m_old = m;
end
end
% plotting the solutions
figure(1);
sgtitle(sprintf('Flow-field variables across nozzle in Non-conservation form \n No. of time steps = %g',nt));
subplot(4,1,1);
plot(x,M);
xlabel('Length');
ylabel('Mach no');
title(sprintf('Mach number'));
subplot(4,1,2);
plot(x,rho,'g');
xlabel('Length');
ylabel('Density');
title(sprintf('Non-dimensional Density'));
subplot(4,1,3);
plot(x,t,'m');
xlabel('Length');
ylabel('Temperature');
title(sprintf('Non-dimensional Temperature'));
subplot(4,1,4);
plot(x,p,'r');
xlabel('Length');
ylabel('Pressure');
title(sprintf('Non-dimensional Pressure'));
% Variarion of flow field variables at nozzle throat
figure(2)
plot(nt_x,rho_throat);
hold on
plot(nt_x,T_throat,'r');
plot(nt_x,p_throat,'g');
plot(nt_x,M_throat,'y');
title(sprintf(['Variation of Flow-field variables at nozzle throat Non-conservation form']));
xlabel('Number of time steps');
legend('Density','Temperature','Pressure','Mach number');
% Mass flow variation across nozzle at different time steps
figure(3)
plot(x,m_initial,x,m_t(50,:),x,m_t(100,:),x,m_t(150,:),x,m_t(200,:),x,m_t(700,:),x,m,'r');
xlabel('Non-dimensional distance through nozzle');
ylabel('Non-dimensional mass flow');
legend('0\Deltat','50\Deltat','100\Deltat','150\Deltat','200\Deltat','700\Deltat','Analytical solution');
title(sprintf('Non-dimensional mass flow rate across the nozzle Non-conservation form'));
end
CONSERVATIVE METHOD:
function [rho,t,p,M,m,count]= conservative_form
% Initialization
n=31;
nt=1400;
C=0.5;
tol=1e-4;
x = linspace(0,3,n);
dx = x(2)-x(1);
A = 1 + 2.2*(x-1.5).^2;
throat = (n+1)/2; % assuming equal no of grid points on either side of throat
nt_x = linspace(1,nt,nt);
gamma = 1.4;
error = 1;
count = 1;
m_old = 0;
m = 0.579*ones(n);
% Initial conditions
u2_i = 0.59; % initial mass flow rate, at t=0;
for i= 1:n
if ((x(i)>=0) && (x(i)<=0.5))
rho(i) = 1;
t(i) = 1;
elseif ((x(i)<1.5) && (x(i)>0.5))
rho(i) = 1 - 0.366*(x(i)-0.5);
t(i) = 1 - 0.167*(x(i)-0.5);
elseif ((x(i)>=1.5) && (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
v(i) = u2_i/(rho(i)*A(i));
u1(i) = rho(i)*A(i);
u2(i) = rho(i)*A(i)*v(i);
u3(i) = rho(i)*A(i)*(t(i)/(gamma-1) + (gamma/2)*v(i)^2);
end
% copy of initial values
u1_initial = u1;
u2_initial = u2;
u3_initial = u3;
% loop to calculate dt
for k = 1:n
dt_x(k) = C*dx/(t(k)^(1/2)+v(k));
end
dt = min(dt_x);
% time loop
for j = 1:nt
u1_old = u1;
u2_old = u2;
u3_old = u3;
% predictor step
for i = 2:n-1
% computing flux and source terms
f1_i = u2(i);
f1_ip1 = u2(i+1);
f2_i = u2(i)^2/u1(i) + ((gamma-1)/gamma)*(u3(i) - (gamma/2)*u2(i)^2/u1(i));
f2_ip1 = u2(i+1)^2/u1(i+1) + ((gamma-1)/gamma)*(u3(i+1) - (gamma/2)*u2(i+1)^2/u1(i+1));
j2 = (1/gamma)*rho(i)*t(i)*((A(i+1) - A(i))/dx);
f3_i = gamma*u2(i)*u3(i)/u1(i) - (gamma*(gamma-1)/2)*u2(i)^3/u1(i)^2;
f3_ip1 = gamma*u2(i+1)*u3(i+1)/u1(i+1) - (gamma*(gamma-1)/2)*u2(i+1)^3/u1(i+1)^2;
du1_dtp(i) = -(f1_ip1 - f1_i)/dx; % continuity eq
du2_dtp(i) = j2 - (f2_ip1 - f2_i)/dx; % momentum eq
du3_dtp(i) = -(f3_ip1 - f3_i)/dx; % energy eq
% updating the predicted values
u1(i) = u1(i) + du1_dtp(i)*dt;
u2(i) = u2(i) + du2_dtp(i)*dt;
u3(i) = u3(i) +du3_dtp(i)*dt;
end
% Updating primitive variables after predictor step
rho = u1./A;
v = u2./u1;
t = (gamma-1)*((u3./u1) - (gamma/2)*v.^2);
% corrector step
for i = 2:n-1
f1_i = u2(i);
f1_ic1 = u2(i-1);
f2_i = u2(i)^2/u1(i) + ((gamma-1)/gamma)*(u3(i) - (gamma/2)*u2(i)^2/u1(i));
f2_ic1 = u2(i-1)^2/u1(i-1) + ((gamma-1)/gamma)*(u3(i-1) - (gamma/2)*u2(i-1)^2/u1(i-1));
j2 = (1/gamma)*rho(i)*t(i)*((A(i) - A(i-1))/dx);
f3_i = gamma*u2(i)*u3(i)/u1(i) - (gamma*(gamma-1)/2)*u2(i)^3/u1(i)^2;
f3_ic1 = gamma*u2(i-1)*u3(i-1)/u1(i-1) - (gamma*(gamma-1)/2)*u2(i-1)^3/u1(i-1)^2;
du1_dtc(i) = -(f1_i - f1_ic1)/dx;
du2_dtc(i) = j2 - (f2_i - f2_ic1)/dx;
du3_dtc(i) = -(f3_i - f3_ic1)/dx;
end
%average value
du1_dt = 0.5*(du1_dtp + du1_dtc);
du2_dt = 0.5*(du2_dtp + du2_dtc);
du3_dt = 0.5*(du3_dtp + du3_dtc);
for i = 2:n-1
% computing corrected values
u1(i) = u1_old(i) + du1_dt(i)*dt;
u2(i) = u2_old(i) + du2_dt(i)*dt;
u3(i) = u3_old(i) + du3_dt(i)*dt;
end
% Inlet BC
u1(1) = rho(1)*A(1);
u2(1) = 2*u2(2) - u2(3);
v(1) = u2(1)/(rho(1)*A(1));
u3(1) = u1(1)*(t(1)/(gamma-1) + (gamma/2)*v(1)^2);
% Outlet 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);
% computing the values of the primitive variables
for i= 1:n
rho(i) = u1(i)/A(i);
v(i) = u2(i)/u1(i);
t(i) = (gamma-1)*(u3(i)/u1(i) - gamma/2*v(i)^2);
end
% Pressure, mass flow rate and Mach number
p = rho.*t;
M = v./t.^0.5;
m = rho.*v.*A;
m_t(j,:) = rho.*v.*A;
% Flow field variables at the throat
rho_throat(j) = rho(throat);
T_throat(j) = t(throat);
p_throat(j) = p(throat);
M_throat(j) = M(throat);
% Convergence
error = max(abs(m-m_old));
if error>tol
count = count + 1;
m_old = m;
end
end
% plotting steady state solutions of flow-field variables
figure(4);
sgtitle(sprintf('Flow-field variables across nozzle Conservation form \n No. of time steps = %g',nt));
subplot(4,1,1);
plot(x,M);
xlabel('Length');
ylabel('Mach number');
title(sprintf('Mach number'));
subplot(4,1,2);
plot(x,rho,'r');
xlabel('Length');
ylabel('Density');
title(sprintf('Non-dimensional Density'));
subplot(4,1,3);
plot(x,t,'g');
xlabel('Length');
ylabel('Temperature');
title(sprintf('Non-dimensional Temperature'));
subplot(4,1,4);
plot(x,p,'m');
xlabel('Length');
ylabel('Pressure');
title(sprintf('Non-dimensional Pressure'));
%Variarion of flow field variables at nozzle throat
figure(5)
plot(nt_x,rho_throat);
hold on
plot(nt_x,T_throat);
plot(nt_x,p_throat);
plot(nt_x,M_throat);
title(sprintf('Variation of flowfield variables at nozzle throat - Conservation form'));
xlabel('Number of time steps');
legend('Density','Temperature','Pressure','Mach number');
% Mass flow variation across nozzle at different time steps
figure(6)
plot(x,u2_initial,'--',x,m_t(50,:),x,m_t(100,:),x,m_t(150,:),x,m_t(200,:),x,m_t(700,:),x,m,':');
xlabel('Non-dimensional distance through nozzle');
ylabel('Non-dimensional mass flow');
legend('0\Deltat','50\Deltat','100\Deltat','150\Deltat','200\Deltat','700\Deltat','Analytical solution');
title(sprintf('Mass flow rate across the nozzle Conservation form'));
end
IMAGES
RESULTS
From the above results we can say that the timewise variation of quantities at the nozzle throat (n=16) is much less for the Conservative form than the Non-Conservative form and although there are slight disturbances during the initial steps of iteration, it soon becomes stable. The results closely match with those given in reference. The profiles match exactly with the reference book for qualitative analysis of flow field properties.
Although the non-conservation form and the conservation form have different initial conditions we can see that the variations in the mass flowrates are quite less in the conservative form and it is quite close to the exact analytical value of 0.579 at the steady state. As evident from the plot the conservative form gives a mass flowrate variation that is close to a constant value and simultaneously close to the exact value of 0.579 at the steady state. The non-conservative form shows some oscillations around the inflow and outflow boundaries and is less satisfactory compared to the conservation form.
For the grid dependency test, we can say that for a grid-size of 31 points the Non-Conservative form gives a more accurate result compared to the Conservative form. However when we change the grid-size to 61 points the values of the flow quantites are quite close when calculated by either of the methods. Additionally both the methods give results which are clsoer to the exact value in this case.
Conservation form produces slightly more accurate results for the flux quantities m′">
and the Non-conservation form produces slightly more accurate results for the primitive variables but the non-conservative form is more accurate overall.
The non-conservative form leads to smaller residuals in computation.
The non-conservative form is faster in comparison to the conservative form when it comes to the computation times. Here conservative form requires marginally more time due to increased number of calculations.
Time for conservative form is 4.0757 s
Time for non-conservative form is also 4.0757 s