Flow simulation of a 1D Super-Sonic nozzle using the Mac-Cormack method
NUMERICAL SOLUTION OF 1D SUPERSONIC NOZZLE FLOW SIMULATION BY MAC-CORMACK METHOD
PROJECT OBJECTIVES:
i. Numerical solution of the governing equations in both conservative and non-conservative forms.
ii. Creating user defined functions for calculating the flow quantities in both Conservative and Non-Conservative forms.
iii. Demonstrating the time-wise variation of density, temperature, velocity and Mach number at the nozzle throat.
iv. Demonstrating the qualitative aspects of Quasi 1D nozzle flow at the steady state for both conservative and non-conservative forms.
v. Comparison of non-dimensional mass flow rates between the conservative and non-conservative forms.
vi. Performing a grid-dependency test for various flow quantites for (grid-sizes n=31 and n = 61)
vii. Comparison of the computation time of solutions by both conservative and non-conservative approaches.
______________________________________________________________________
(I) INTRODUCTION:
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 crosss sectional area of the nozzle is very large so that the velocity is almost `V~=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'`,`rho'`,`x'`,`V'`,`t'`and `A'` as:
Temperature: `T' = T/T_0`
Density: `rho' = rho/rho_0`
Non-dimensional Distance: `x'=x/L`
Non-dimensional Velocity: `V' = V/a_0` where `a_0` is the local sonic velocity (`a_0 = sqrt(gammaRT_0)`)
Non-dimensional Time: `t' = t/(L/a_0)` where `L/a_0` is the characteristic time-scale of the problem.
Non-dimensional Area: `A'=A/(A^star)`
GOVERNING EQUATIONS:
1. Non-Conservtive Form:
i. Continuity Equation:
`(delrho')/(delt') = rho'(delV')/(delx') - rho'V'(del(lnA))/(delx') - V'(delp')/(delx')`
ii. Momentum Equation:
`(delV')/(delt') = V'(delV')/(delx') - 1/gamma((delT')/(delx') + (T')/(rho')(delrho')/(delx'))`
iii. Energy Equation:
`(delT')/(delt') = V'(delT')/(delx') - (gamma-1)T'[(delV')/(delx') + V'(del(lnA'))/(delx')]`
2. Conservative Form:
i. Continuity Equation:
`(del(rho'A'))/(delt') + (del(rho'A'V'))/(delx') = 0`
ii. Momentum Equation:
`(del(rho'A'V'))/(delt') + (del(rho'A'V'^2+1/gammap'A'))/(delx')=1/gammap'(delA')/(delx')`
iii. Energy Equation:
`(del[rho'((T')/(gamma-1)+gamma/2V'^2)A'])/(delt') + (del[rho'((T')/(gamma-1) +gamma/2V'^2)V'A' + p'A'V'])/(delx') = 0`
The following governing 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 the following terms :
Solution Vector: `U`
Flux Vector: `F`
Source Terms: `J`
Hence we write the equations (continuity, momentum and energy) in these terms:
Solution Vectors:
`U1 = rho'A'`
`U2 = rho'A'V'`
`U3 = rho'((T')/(gamma-1)+gamma/xV'^2)A'`
Flux Vectors:
`F1 = rho'A'V'`
`F2 = rho'A'V'^2+1/gammap'A'`
`F3 = rho'((T')/(gamma-1) +gamma/2V'^2)V'A' + p'A'V'`
Source Terms:
`J2 = 1/gammap'(delA')/(delx')`
However if the Flux terms are written in terms of the primitive variables (`rho'`,`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) + (gamma-1)/gamma(U3 - gamma/2(U2)^2/(U1))`
`F3 = gamma(U2U3)/U1 - gamma(gamma-1)/2(U2)^3/(U1)^2`
Source Term: `J2 = (gamma-1)/gamma(U3 - gamma/2(U2)^2/(U1))(del(lnA'))/(delx')`
Decomposition to the primitive variables (`rho'`,`T'`,`V'`) we have:
Density: `rho = (U1)/(A')`
Velocity: `V' = (U2)/(U1)`
Temperature: `T' = (gamma-1)((U3)/(U1)-gamma/2V'^2)`
Pressure: `p'=rho'T'`
In the problem specification we have taken the number of timesteps as 1400 as mentioned in Chapter 7 of the book: John D. Anderson Jr. Computational Fluid Dynamics: The Basics with Applications. Other problem specifications were taken directly from the above book and the results obtained was used to compare and analyze the solutions obtained form the programs.
(II) MATLAB PROGRAMS:
1. MAIN SYNTAX: (DRIVER SCRIPT)
%% MATLAB Program to simulate the 1D Supersonic Nozzle flow using Mc-Cormack Method:
%% Analysis of Conservative and Non-conservative forms of governing equations and grid dependancy test:
clear all;
close all;
clc;
%% Defining the initial conditions:
L = 3; %Length of the nozzle
n = 61; %Manually entering the number of grid-points
x = linspace(0,L,n); % Initializing the 1D grid
dx = L/(n-1); % Grid spacings
gamma = 1.4; % Ratio of Cp and Cv for air assuming it to be ideal
nt = 1400; % number of timesteps
throat = (n-1)/2; % Throat point
a = 1 + 2.2*(x - 1.5).^2; % Area calculation
C = 0.5; % Courant Number (specified for the problem)
%% Solution by Non-Conservative form:
tic
[m, M, v, Temp, p, rho, Temp_throat, rho_throat, p_throat, M_throat] = non_conservative_form(x,dx,n,C,a,nt,gamma,throat);
time_elapsed = toc;
fprintf('The Computation time by Non-Conservative form: %0.4f',time_elapsed);
%% Time-wise variation of density, temperature, pressure and Mach number for the Non-Conservative form:
time_variation = figure(2);
subplot(4,1,1);
plot(rho_throat,'b');
title('Time-wise variations in quantity in Non-Conservative form');
legend('Density');
ylabel('\rho/\rho_{o}');
axis([0 1400 0 1.3]);
grid minor;
subplot(4,1,2);
plot(Temp_throat,'r');
legend('Temperature');
ylabel('T/T_{o}');
axis([0 1400 0.6 1]);
grid minor;
subplot(4,1,3);
plot(p_throat,'g');
legend('Pressure');
ylabel('p/p_{o}');
axis([0 1400 0.4 1.1]);
grid minor;
subplot(4,1,4);
plot(M_throat,'m');
legend('Mach Number');
ylabel('M');
axis([0 1400 0.6 1.4]);
grid minor;
xlabel('Number of iterations');
%% Qualitative Aspects of 1D flow in Non-Conservative Form:
qualitative_non_cons = figure(3);
subplot(4,1,1);
plot(x,M,'-b');
title('Qualitative Aspects of 1D Supersonic Nozzle flow: Non-Conservative form');
ylabel('M');
legend('Mach Number','location','southeast');
grid minor;
subplot(4,1,2);
plot(x,p,'-r');
ylabel('p/p_{o}');
legend('Pressure ratio');
grid minor;
subplot(4,1,3);
plot(x,rho,'-g');
ylabel('\rho/\rho_{o}');
legend('Density ratio','outside');
grid minor;
subplot(4,1,4);
plot(x,Temp,'-m');
ylabel('T/T_{o}');
legend('Temperature ratio','outside');
grid minor;
%% Solution by Conservative form:
tic
[m_c, M_c, v_c, Temp_c, p_c, rho_c, Temp_throat_c, rho_throat_c, p_throat_c, M_throat_c] = conservative_form(x,dx,n,C,a,nt,gamma,throat);
time_elapsed_cons = toc;
fprintf('The computation time by Conservative form: %0.4f',time_elapsed_cons);
%% Time-wise variation of density, temperature, pressure and Mach number for the Non-Conservative form:
time_variation_cons = figure(5);
subplot(4,1,1);
plot(rho_throat_c,'b');
title('Time-wise variations in quantity in Conservative Form');
legend('Density');
ylabel('\rho/\rho_{o}');
axis([0 1400 0 1.3]);
grid minor;
subplot(4,1,2);
plot(Temp_throat_c,'r');
legend('Temperature');
ylabel('T/T_{o}');
axis([0 1400 0.6 1]);
grid minor;
subplot(4,1,3);
plot(p_throat_c,'g');
legend('Pressure');
ylabel('p/p_{o}');
axis([0 1400 0.4 1.1]);
grid minor;
subplot(4,1,4);
plot(M_throat_c,'m');
legend('Mach Number');
ylabel('M');
axis([0 1400 0.6 1.4]);
grid minor;
xlabel('Number of iterations');
%% Qualitative Aspects of 1D flow in Conservative Form:
qualitative_cons = figure(6);
subplot(4,1,1);
plot(x,M_c,'-b');
title('Qualitative Aspects of 1D Supersonic Nozzle flow: Conservative form');
ylabel('M');
legend('Mach Number','location','southeast');
grid minor;
subplot(4,1,2);
plot(x,p_c,'-r');
ylabel('p/p_{o}');
legend('Pressure ratio');
grid minor;
subplot(4,1,3);
plot(x,rho_c,'-g');
ylabel('\rho/\rho_{o}');
legend('Density ratio','outside');
grid minor;
subplot(4,1,4);
plot(x,Temp_c,'-m');
ylabel('T/T_{o}');
legend('Temperature ratio','outside');
grid minor;
%% Comparison of Mass-Flow rates by Conservative and Non-Conservative methods and comparison with Exact Form:
mass_comparison = figure(7);
hold on;
plot(x, m, 'g','LineWidth', 1.75);
plot(x, m_c, 'm','LineWidth', 1.75);
line([0 3], [0.579 0.579],'LineWidth', 1.75);
legend('Conservative Form', 'Non-Conservative Form', 'Exact Solution');
grid minor;
xlabel('x/L');
ylabel('\rhoVA/\rho_{o}V_{o}A_{0}');
title('Non-Dimensional Mass Flow comparison with Exact Solution');
2. NON-CONSERVATIVE FORM OF THE GOVERNING EQUATION: (USER-DEFINED FUNCTION)
function [m, M, v, Temp, p, rho, Temp_throat, rho_throat, p_throat, M_throat] = non_conservative_form(x,dx,n,C,a,nt,gamma,throat)
% Initializing the variables:
rho = 1 - 0.3146*x;
Temp = 1- 0.2314*x;
v = (0.1 + 1.09*x).*Temp.^0.5;
%Outer Time Loop
for k = 1:nt
rho_old = rho;
Temp_old = Temp;
v_old = v;
%Adaptive time-step control
dt = min(C.*dx./(sqrt(Temp)+v));
%Predictor Part followed by Solution Update:
for j = 2: n-1
dvdx = (v(j+1)-v(j))/dx;
drhodx = (rho(j+1) - rho(j))/dx;
dTempdx = (Temp(j+1) - Temp(j))/dx;
dlog_a_dx = (log(a(j+1))-log(a(j)))/dx;
% Predictor Step:
%Continuity Equation:
drhodt_p(j) = -rho(j)*dvdx - rho(j)*v(j)*dlog_a_dx - v(j)*drhodx;
% Momentum Equation:
dvdt_p(j) = -v(j)*dvdx - (1/gamma)*(dTempdx + (Temp(j)/rho(j))*drhodx);
% Energy Equation:
dTdt_p(j) = -v(j)*dTempdx - (gamma-1)*Temp(j)*(dvdx + v(j)*dlog_a_dx);
% Solution Update step:
rho(j) = rho(j) + drhodt_p(j)*dt; %Density Update
v(j) = v(j) + dvdt_p(j)*dt; %Velocity Update
Temp(j) = Temp(j) + dTdt_p(j)*dt; %Temperature Update
end
% Corrector step
for j = 2: n-1
dvdx = (v(j)-v(j-1))/dx;
drhodx = (rho(j) - rho(j-1))/dx;
dTempdx = (Temp(j) - Temp(j-1))/dx;
dlog_a_dx = (log(a(j))-log(a(j-1)))/dx;
%Continuity Equation:
drhodt_c(j) = -rho(j)*dvdx - rho(j)*v(j)*dlog_a_dx - v(j)*drhodx;
% Momentum Equation:
dvdt_c(j) = -v(j)*dvdx - (1/gamma)*(dTempdx + (Temp(j)/rho(j))*drhodx);
% Energy Equation:
dTdt_c(j) = -v(j)*dTempdx - (gamma-1)*Temp(j)*(dvdx + v(j)*dlog_a_dx);
end
% Compute the average time derivative
drhodt = 0.5*(drhodt_p + drhodt_c);
dvdt = 0.5*(dvdt_p + dvdt_c);
dTdt = 0.5*(dTdt_p + dTdt_c);
% Final solution Update:
for i = 2:n-1
rho(i) = rho_old(i) + drhodt(i)*dt;
v(i) = v_old(i) + dvdt(i)*dt;
Temp(i) = Temp_old(i) + dTdt(i)*dt;
end
% Applying the boundary conditions
% Inlet
v(1) = 2*v(2) - v(3);
rho(1) = 1;
Temp(1) = 1;
% Outlet
v(n) = 2*v(n-1) - v(n-2);
rho(n) = 2*rho(n-1) - rho(n-2);
Temp(n) = 2*Temp(n-1) - Temp(n-2);
% Defining the values for pressure, mass flow rate and mach number:
p = rho.*Temp; %Non-Dimensional Pressure
m = rho.*a.*v; %Non-Dimensional Mass flow rate
M = v./sqrt(Temp); %Mach Number
%Calculation of Quantities at the Throat point of the nozzle:
Temp_throat(k) = Temp(throat);
p_throat(k) = p(throat);
rho_throat(k) = rho(throat);
M_throat(k) = M(throat);
% Comparison of non-dimensional mass-flow rates at different times:
mass = figure(1);
hold on;
if k == 50
plot(x,m,'b','linewidth',1.75);
else if k == 100
plot(x,m,'g','linewidth',1.75);
else if k == 150
plot(x,m,'r','linewidth',1.75);
else if k == 200
plot(x,m,'m','linewidth',1.75);
else if k == 700
plot(x,m,'-','linewidth',1.75);
plot(x(end),m(end),'b*');
end
end
end
end
end
end
title('Comparison of non-dimensional mass flow rates at different time-instants');
xlabel('x/L (Non-dimensional distance through the nozzle)');
ylabel('\rhoVA/\rho_{o}V_{o}A_{o}');
legend('50\Deltat','100\Deltat','150\Deltat','200\Deltat','700\Deltat');
axis([0 3 0 2.5]);
grid on;
grid minor;
end
3. CONSERVATIVE FORM OF THE GOVERNING EQUATIONS: (USER-DEFINED FUNCTION)
function [m_c, M_c, v_c, Temp_c, p_c, rho_c, Temp_throat_c, rho_throat_c, p_throat_c, M_throat_c] = conservative_form(x,dx,n,C,a,nt,gamma,throat)
%% Initialization step:
% Initializing density and temperature:
for i = 1:length(x)
if x(i) >= 0 && x(i) < 0.5
rho_c(i) = 1;
Temp_c(i) = 1;
else if x(i) >= 0.5 && x(i) < 1.5
rho_c(i) = 1.0 -0.366*(x(i) - 0.5);
Temp_c(i) = 1.0 - 0.167*(x(i) - 0.5);
else if x(i) >= 1.5 && x(i) <= 3.0
rho_c(i) = 0.634 - 0.3879*(x(i) - 1.5);
Temp_c(i) = 0.833 - 0.3507*(x(i) - 1.5);
end
end
end
end
% Initializing velocity and pressure:
v_c = 0.59./(rho_c.*a);
p_c = rho_c.*Temp_c;
% Initializing the solution vectors:
U1 = rho_c.*a;
U2 = rho_c.*a.*v_c;
U3 = rho_c.*a.*((Temp_c./(gamma - 1)) + (gamma/2)*v_c.^2);
% Time Iteration Loop
for k = 1:nt
% Preserving the old solution vectors:
U1_old = U1;
U2_old = U2;
U3_old = U3;
% Adaptive Time-step control:
dt = min(C.*dx./(sqrt(Temp_c) + v_c));
% Defining the 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);
%% Solution by Mc-Cormack Method:
for j = 2:n-1
% Defining the source term:
J2(j) = (1/gamma).*rho_c(j).*Temp_c(j).*((a(j+1)-a(j))/dx);
% Predictor step:
dU1dt_p(j) = -(F1(j+1) - F1(j))/dx;
dU2dt_p(j) = -((F2(j+1) - F2(j))/dx) + J2(j);
dU3dt_p(j) = -(F3(j+1) - F3(j))/dx;
% Updating the new values of solution vectors:
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 the primitive variables and flux vectors:
rho_c = U1./a;
v_c = U2./U1;
Temp_c = (gamma -1)*(U3./U1 - (gamma/2)*v_c.^2);
p_c = rho_c.*Temp_c;
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 j = 2:n-1
% Updating the source term:
J2(j) = (1/gamma)*rho_c(j).*Temp_c(j)*((a(j)-a(j-1))/dx);
% Corrector step:
dU1dt_c(j) = -(F1(j) - F1(j-1))/dx;
dU2dt_c(j) = -((F2(j) - F2(j-1))/dx) + J2(j);
dU3dt_c(j) = -(F3(j) - F3(j-1))/dx;
end
% Updating the solution vectors after corrector step:
% Averaging the predictor and corrector steps:
dU1dt = 0.5*(dU1dt_p + dU1dt_c);
dU2dt = 0.5*(dU2dt_p + dU2dt_c);
dU3dt = 0.5*(dU3dt_p + dU3dt_c);
% Calculating the new solution vectors:
for j = 2:n-1
U1(j) = U1_old(j) + dU1dt(j)*dt;
U2(j) = U2_old(j) + dU2dt(j)*dt;
U3(j) = U3_old(j) + dU3dt(j)*dt;
end
% Applying the boundary conditions:
% left side/ inlet conditions
% U1(1) = rho_c(1).*Temp_c(1);
U2(1) = 2*U2(2) - U2(3);
% U3(1) = U1(1)*(Temp_c(1)/(gamma-1)+(gamma/2)*v_c(1).^2);
% right side/outlet 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);
% Final updation of the primitive variables:
rho_c = U1./a;
v_c = U2./U1;
Temp_c = (gamma -1)*(U3./U1 - (gamma/2)*v_c.^2);
p_c = rho_c.*Temp_c;
%%
% Defining the Values of mach Number and mass flow rate:
M_c = v_c./sqrt(Temp_c);
m_c = rho_c.*a.*v_c;
% Calculating values of specific quantites at the throat:
M_throat_c(k) = M_c(throat); % Mach Number
p_throat_c(k) = p_c(throat); % Pressure Ratio
Temp_throat_c(k) = Temp_c(throat); % Non-dimensional Temperature
rho_throat_c(k) = rho_c(throat); % Density Ratio
% Comparison of non-dimensional mass-flow rates at different times:
mass2 = figure(4);
hold on;
if k == 100
plot(x,m_c,'b','linewidth',1.75);
else if k == 200
plot(x,m_c,'g','linewidth',1.75);
else if k == 700
plot(x,m_c,'r','linewidth',1.75);
end
end
end
end
title('Comparison of non-dimensional mass flow rates at different time-instants (Conservative Form)');
xlabel('x/L (Non-dimensional distance through the nozzle)');
ylabel('\rhoVA/\rho_{o}V_{o}A_{o}');
legend('100\Deltat','200\Deltat','700\Deltat');
axis([0 3 0.4 1.6]);
grid on;
grid minor;
end
(III) RESULTS:
1. Time-wise variation of pressure, temperature, velocity and Mach number at the nozzle throat:
(a) Non-Conservative Form: (Grid-size n=31)
(b) Conservative Form: (Grid-size n = 31)
From the above results we can say that the timewise variation of quantities at the nozzle throat (n=15) 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.
2. Qualitative aspects of Quasi 1D nozzle flow at steady state:
(a). Non-Conservative form: (Grid-size n = 31)
(b). Conservative Form: (Grid-size n = 31)
The profiles match exactly with the reference book.
3. Non-Dimensional mass flowrate distribution at various time-steps:
(a). Non-Conservative form: (Grid-size n = 31)
(b). Conservative Form: (Grid-size n = 31)
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.
4. Comparison of the non-dimensional mass flowrates of both the forms:
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.
5. Grid-dependancy test:
Here we would be examining the value of flow quantities by changing the value of grid-size. The results are summarized below:
From the above results we can see that for a grid-size of 30 points the Non-Conservative form gives a more accurate result compared to the Conservative form. However when we change the grid-size to 60 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.
Plot of solutions for the grid-size (n = 61)
Here we are comparing the plots of time-wise variations of non-dimensional mass flowrates and the flow quantities at the throat of the nozzle.
(a). Non-Conservative Form:
(b). Conservative Form:
6. Analysis of Computation time and number of timesteps required for convergence:
(IV) CONCLUSIONS AND REMARKS:
1. We used an adaptive time-step control for the following functions. it was calculated from the following formula: `Deltat = min(C(Deltax)/(a_i+V_i))`
Where `C` = Courant Number (CFL Criterion) = 0.5 specified in the problem
`a` = local sonic velocity `a = sqrt(T)`
`V` = local flow velocity
2. Mach Number was calculated as `M = V/sqrt(T)`
3. Conservation form produces slightly more accurate results for the flux quantities (Ex. `m'`) and the Non-conservation form produces slightly more accurate results for the primitive variables (`rho'`,`T'`,`V'`) but the non-conservative form is more accurate overall.
4. The non-conservative form leads to smaller residuals in computation.
5. 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.