Simulation of a 1-D Super-sonic nozzle flow simulation using Macormack Method
Simulation of a 1D Super-sonic nozzle flow simulation using Macormack Method
Objectives:
- Simulation of 1D Super-sonic nozzle flow using the governing equation in both Conservative and Non-Conservative form.
- Demonstrating time-wise variation of density, temperature, velocity and Mach number at the throat.
- Demonstrating qualitative aspects of Quasi 1D nozzle flow.
- Compare normalized mass flow rate between the conservative and non-conservative form.
- Performing grid-dependency test.
- Comparison of the computation time of solutions by both conservative and non-conservative approaches.
------------------------------------------------------------------------------------------------------------
Here, we consider the steady, isentropic flow through a convergent-divergent nozzle. The flow at the inlet to the nozzle comes from a reservoir where the pressure and temperature are stagnation values as the cross-sectional area of the reservoir is large and hence the velocity is small. The flow expands isentropically to supersonic speeds at the nozzle exit. The flow is locally subsonic in the convergent section of the nozzle, sonic at the throat and supersonic at the divergent section.
In reality, the flow field is two dimensional but here we are making the assumption that the flow properties vary only in the x-direction. Such flow is defined as the quasi-one-dimensional flow.
Following are the governing equations in Non-Conservative form:
Continuity equation:
`(-drho\')/(dt\')= -rho\'(dv\')/(dx\')-rho\'v\'(d(lnA\'))/(dx\')-v\'(dp\')/(dx\')`
Momentum equation:
`(dv\')/(dt\')=-v\'(dv\')/(dx\')-1/gamma ((dT\')/(dx\')+(T\')/rho\'(dp\')/(dx\'))`
Energy Equation:
`(dT\')/(dt\')=-v\'(dT\')/(dx\')-(gamma-1)T\'((dv\')/(dx\')+v\'(d(lnA\'))/(dx\'))`
Following are the governing equations in Conservative form:
Contuinity Equation:
`(d(rho\'A\'))/(dt\')+(d(rho\'A\'v\'))/(dx\')=0`
Momentum Equation:
`(d(rho\'A\'v\'))/(dt\')+(d(rho\'A\'v\'^2+1/gammap\'A\'))/(dx\')=1/gammap\'(dA\')/(dx\')`
Energy Equation:
`(d(rho\'A\'((e\')/(gamma-1)+gamma/2v\'^2)))/(dt\')+(d(rho\'v\'A\'((e\')/(gamma-1)+gamma/2v\'^2)+p\'A\'v\'))/(dx\')=0`
The above equations in conservative and non-conservative form are expresed in non-dmensional terms. These terms are explained below:
Non-dimensional temperature: `T\'=T/T_0` where `T_0`= reservoir temperature.
Non-dimensional density: `rho\'=rho/rho_0` where `rho_0`= reservoir density.
Non-dimensional length: `x/L` where L = length of nozzle.
Speed of sound in the reservoir is `a_0` where `a_0=sqrt(gammaRT_0`
Non-dimensional velocity: `V\'=V/a_0`
`L/a_0` has the dimension of time so non-dimensional time is defined as: `t\'=t/(L/a_0)`
Non-dimensional area: `A\'=A/(A^#)` where `A^#` = area at throat.
The non-dimensional terms presented above vary from 0 to1 which is an aesthetic advantage when presenting the results.
The conservative form of governing equation will be needed to be written in the generic form so they can be solved numerically. For this, we introduce the following terms:
Solution Vector: U
Flux Vector: F
Source term: J
`U1=rho\'A\'`
`U_2=rho\'A\'V\'`
`U_3=rho\'((e\')/(gamma-1)+gamma/2V\'^2)A\'`
`F_1=rho\'A\'V\'`
`F_2=rho\'A\'V\'+1/gammap\'A\'`
`F_3= rho\'((e\')/(gamma-1)+gamma/2V\'^2)V\'A\'+p\'A\'V\'`
`J=1/gammap\'(dA\')/(dx\')`
When the flux vectors are written in primitive variables(`rho,T,V`) then numerical difficulties arise occasionally, hence the flux vectors and source terms are expressed using solution vectors.
`F_1=U_1`
`F_2=(U_2^2)/U_1+(gamma-1)/gamma(U_3-gamma/2(U_2)/U_1)`
`F_3=gamma(U_2U_3)/U_1-gamma(gamma-1)/2(U_2^3)/U_1^2`
`J=(gamma-1)/gamma(U-3-gamma/2(U_2^2)/U_1)(d(lnA))/(dx\')`
The conservative form of equations can be written using the above terms as follows:
Continuity Equation: `(dU_1)/(dt\')=-(dF_1)/(dx\')`
Momentum Equation: `(dU_2)/(dt\')=-(dF_2)/(dx\')+J`
Energy Equation: `(dU_3)/(dt\')=-(dF_3)/(dx\')`
All the other specifications were take from the Chapter 7 of the book: John D. Anderson Jr. Computational Fluid Dynamics: The Basics with Applications.
Matlab Programme:
Main Programme-
clear all
close all
clc
%% Solvig 1D nozzle flow using conservative form of govening equation
% Inputs
n = 31; % number of grid points
l = 3; % length of the nozzle
x = linspace(0,l,n); % grig
dx = x(2)-x(1); % spacing between grids
n_throat = (n-1)/2; %throat point
gamma = 1.4; % Ratio of Cp and Cv of air assuming ideal conditions
C = 0.5; % Courant Number
a = 1+2.2*(x-1.5).^2; % area profile for nozzle
nt = 1400;
%%------Solution for Non Conservative Form----------%%
tic
[rho,v,temp,p,M,m,rho_throat,temp_throat,p_throat,M_throat] = non_conservative_form(n,l,x,dx,gamma,C,a,n_throat);
timer = toc;
fprintf(\'Computational time = %d\',timer);
%%----Properties at throat----%%
figure(2);
subplot(4,1,1);
plot(rho_throat,\'b\',\'linewidth\',1.5);
title(\'Time-wise variations in quantity in Non-Conservative form\');
legend(\'Density\');
ylabel(\'rho/rho_{o}\');
axis([0 1400 0 1.3]);
subplot(4,1,2);
plot(temp_throat,\'r\',\'linewidth\',1.5);
legend(\'Temperature\');
ylabel(\'temp/temp_{o}\');
axis([0 1400 0.6 1]);
subplot(4,1,3);
plot(p_throat,\'g\',\'linewidth\',1.5);
legend(\'Pressure\');
ylabel(\'p/p_{o}\');
axis([0 1400 0.4 1.1]);
subplot(4,1,4);
plot(M_throat,\'m\',\'linewidth\',1.5);
legend(\'Mach Number\');
ylabel(\'M\');
axis([0 1400 0.6 1.4]);
xlabel(\'Number of iterations\');
%%-----Qualitative Aspects of 1D flow in Non-Conservative Form------%%
figure(3)
subplot(4,1,1);
plot(x,rho,\'m\',\'linewidth\',1.5);
title(\'Properties profiles in the 1D Supersonic Nozzle flow (Non-Conservative Form)\');
ylabel(\'rho/rho_{o}\');
legend(\'Density ratio\');
subplot(4,1,2);
plot(x,temp,\'r\',\'linewidth\',1.5);
ylabel(\'temp/temp_{o}\');
legend(\'Temperature ratio\');
subplot(4,1,3);
plot(x,p,\'g\',\'linewidth\',1.5);
ylabel(\'p/p_{o}\');
legend(\'Pressure ratio\');
subplot(4,1,4);
plot(x,M,\'b\',\'linewidth\',1.5);
ylabel(\'M\');
legend(\'Mach Number\',\'location\',\'northwest\');
%%---------Solving for Conservative form----------%%
tic
[ Mach_no, v_cf, temp_cf, p_cf, rho_cf, temp_throat_cf, rho_throat_cf, p_throat_cf,Mach_no_throat,mass_flowrate ] = conservative_form(n,l,x,dx,gamma,C,a,n_throat)
timer = toc;
fprintf(\'Computational time = %d\',timer);
%%----Properties at throat----%%
figure(5);
subplot(4,1,1);
plot(rho_throat_cf,\'b\',\'linewidth\',1.5);
title(\'Time-wise variations in quantity in Conservative form\');
legend(\'Density\');
ylabel(\'rho/rho_{o}\');
axis([0 1400 0 1.3]);
subplot(4,1,2);
plot(temp_throat_cf,\'r\',\'linewidth\',1.5);
legend(\'Temperature\');
ylabel(\'temp/temp_{o}\');
axis([0 1400 0.6 1]);
subplot(4,1,3);
plot(p_throat_cf,\'g\',\'linewidth\',1.5);
legend(\'Pressure\');
ylabel(\'p/p_{o}\');
axis([0 1400 0.4 1.1]);
subplot(4,1,4);
plot(Mach_no_throat,\'m\',\'linewidth\',1.5);
legend(\'Mach Number\');
ylabel(\'M\');
axis([0 1400 0.6 1.4]);
xlabel(\'Number of iterations\');
%%-----Qualitative Aspects of 1D flow in Non-Conservative Form------%%
figure(6)
subplot(4,1,1);
plot(x,rho_cf,\'m\',\'linewidth\',1.5);
title(\'Properties profiles in the 1D Supersonic Nozzle flow (Conservative Form)\');
ylabel(\'rho/rho_{o}\');
legend(\'Density ratio\');
subplot(4,1,2);
plot(x,temp_cf,\'r\',\'linewidth\',1.5);
ylabel(\'temp/temp_{o}\');
legend(\'Temperature ratio\');
subplot(4,1,3);
plot(x,p_cf,\'g\',\'linewidth\',1.5);
ylabel(\'p/p_{o}\');
legend(\'Pressure ratio\');
subplot(4,1,4);
plot(x,Mach_no,\'b\',\'linewidth\',1.5);
ylabel(\'M\');
legend(\'Mach Number\',\'location\',\'northwest\');
figure(7);
hold on;
plot(x, m, \'g\',\'LineWidth\', 1.75);
plot(x, mass_flowrate, \'m\',\'LineWidth\', 1.75);
legend(\'Conservative Form\', \'Non-Conservative Form\');
grid minor;
xlabel(\'x/L\');
ylabel(\'rhoVA/rho_{o}V_{o}A_{0}\');
title(\'Non-Dimensional Mass Flow comparison with Exact Solution\');Non-Conservative Form Solver:
function [rho,v,temp,p,M,m,rho_throat,temp_throat,p_throat,M_throat,drhodt_throat,dvdt_throat,k] = non_conservative_form(n,l,x,dx,gamma,C,a,n_throat,error,tol)
% Calculate the initial profiles
rho = 1-0.3146*x;
temp = 1-0.2314*x; % t=temperature
v = (0.1+1.09*x).*temp.^0.5;
k = 0; %Initializing time loop
% Time step
% Outer time loop
while error>tol;
k = k+1;
rho_old = rho;
v_old = v;
temp_old = temp;
m_old = rho_old.*v_old.*a;
% Adaptive time control step
timestep = C.*(dx./(sqrt(temp_old)+v_old));
dt = min(timestep);
%Predictor method
for j=2:n-1;
dvdx = (v(j+1)-v(j))/dx;
drhodx = (rho(j+1)-rho(j))/dx;
dtdx = (temp(j+1)-temp(j))/dx;
dlogadx = (log(a(j+1))-log(a(j)))/dx;
%Continuity equation
drhodt_p(j) = -rho(j)*dvdx-rho(j)*v(j)*dlogadx-v(j)*drhodx;
%Momentum equation
dvdt_p(j) = -v(j)*dvdx-(1/gamma)*(dtdx+(temp(j)/rho(j))*drhodx);
%Energy equation
dtdt_p(j) = -v(j)*dtdx-(gamma-1)*temp(j)*(dvdx+(v(j)*(dlogadx)));
% Solution Update
rho(j) = rho(j)+drhodt_p(j)*dt;
v(j) = v(j)+dvdt_p(j)*dt;
temp(j) = temp(j)+dtdt_p(j)*dt;
end
%Corrector Step
for j=2:n-1;
dvdx = (v(j)-v(j-1))/dx;
drhodx = (rho(j)-rho(j-1))/dx;
dtdx = (temp(j)-temp(j-1))/dx;
dlogadx = (log(a(j))-log(a(j-1)))/dx;
%Continuity equation
drhodt_c(j) = -rho(j)*dvdx-rho(j)*v(j)*dlogadx-v(j)*drhodx;
%Momentum equation
dvdt_c(j) = -v(j)*dvdx-(1/gamma)*(dtdx+(temp(j)/rho(j))*drhodx);
%Energy equation
dtdt_c(j) = -v(j)*dtdx-(gamma-1)*temp(j)*(dvdx+(v(j)*(dlogadx)));
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 Udpate
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
%Boundary Condition
%Inlet
v(1) = 2*v(2)-v(3);
%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);
%plot(x,v)
%title(sprintf(\'Time step = %d\',k));
%pause(0.03);
% 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(n_throat);
p_throat(k) = p(n_throat);
rho_throat(k) = rho(n_throat);
M_throat(k) = M(n_throat);
drhodt_throat(k) = drhodt(n_throat);
dvdt_throat(k) = dvdt(n_throat);
error = max(abs(m-m_old));
% 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,\'-.b\',\'linewidth\',1.75);
legend(\'50\\Deltat\',\'100\\Deltat\',\'150\\Deltat\',\'200\\Deltat\',\'700\\Deltat\');
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}\');
axis([0 3 0 2.5]);
grid on;
grid minor;
end
end
Conservative Form Solver:
%% Solution by Mac-Cormack method %%
function [ Mach_no, v_cf, temp_cf, p_cf, rho_cf, temp_throat_cf, rho_throat_cf, p_throat_cf,Mach_no_throat,mass_flowrate,k ] = conservative_form(n,l,x,dx,gamma,C,a,n_throat,error2,tol)
% Providing initial condition
% for density and temperature
for i = 1:length(x);
if x(i) >= 0 && x(i) < 0.5
rho_cf(i) = 1;
temp_cf(i) = 1;
else if x(i) >= 0.5 && x(i) < 1.5
rho_cf(i) = 1.0 -0.366*(x(i) - 0.5);
temp_cf(i) = 1.0 - 0.167*(x(i) - 0.5);
else if x(i) >= 1.5 && x(i) <= 3.0
rho_cf(i) = 0.634 - 0.3879*(x(i) - 1.5);
temp_cf(i) = 0.833 - 0.3507*(x(i) - 1.5);
end
end
end
end
% for velocity and pressure
v_cf = 0.59./(rho_cf.*a);
p_cf = rho_cf.*temp_cf;
% initializing solution vectors
U1 = rho_cf.*a;
U2 = rho_cf.*a.*v_cf;
U3 = rho_cf.*a.*((temp_cf./(gamma-1)) + ((gamma/2)*(v_cf.^2)));
%Adaptive time control Step
timestep = C*dx./(sqrt(temp_cf)+v_cf);
dt_cf = min(timestep);
k = 0; % Initializing time step
% Time loop
while error2>tol;
k = k+1;
mass_flowrate_old = rho_cf.*v_cf.*a;
% preserving old values
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);
%% Predictor Step %%
for j = 2:n-1
%Defining source term
J2(j)=(1/gamma).*rho_cf(j).*temp_cf(j).*((a(j+1)-a(j))/dx);
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;
%Solution Update
U1(j) =U1(j) + dU1dt_p(j)*dt_cf;
U2(j) =U2(j) + dU2dt_p(j)*dt_cf;
U3(j) =U3(j) + dU3dt_p(j)*dt_cf;
end
% Updating primitive variables and Flux vectors
rho_cf = U1./a;
v_cf = U2./U1;
temp_cf = (gamma -1)*(U3./U1 - (gamma/2)*v_cf.^2);
p_cf = rho_cf.*temp_cf;
%Flux vector
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);
%% Corrector Step %%
for j = 2:n-1
% Updating source
J2(j)=(1/gamma)*rho_cf(j).*temp_cf(j).*((a(j)-a(j-1))/dx);
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
% Computing average time derivative
dU1dt_avg = 0.5*(dU1dt_p+dU1dt_c);
dU2dt_avg = 0.5*(dU2dt_p+dU2dt_c);
dU3dt_avg = 0.5*(dU3dt_p+dU3dt_c);
% Final Solution update for solution vectors
for z = 2:n-1
U1(z) = U1_old(z) + dU1dt_avg(z)*dt_cf;
U2(z) = U2_old(z) + dU2dt_avg(z)*dt_cf;
U3(z) = U3_old(z) + dU3dt_avg(z)*dt_cf;
end
%% Applying boundary conditions
% At inlet
U1(1) = rho_cf(1)*a(1);
U2(1) = 2*U2(2)-U2(3);
U3(1) = U1(1)*(temp_cf(1)/(gamma-1) + (gamma/2)*v_cf(1)^2);
%At 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);
% Updating primitive variables
rho_cf=U1./a;
v_cf=U2./U1;
temp_cf=(gamma -1)*(U3./U1-(gamma/2)*v_cf.^2);
p_cf=rho_cf.*temp_cf;
% specifying mass flow rate and Mach number
mass_flowrate = rho_cf.*a.*v_cf;
Mach_no = v_cf./sqrt(temp_cf);
%computing values at throat
rho_throat_cf(k)=rho_cf(n_throat);
temp_throat_cf(k)=temp_cf(n_throat);
p_throat_cf(k)=p_cf(n_throat);
Mach_no_throat(k)=Mach_no(n_throat);
error2 = max(abs(mass_flowrate-mass_flowrate_old));
% Ploting mass flow rate at different time step
figure(5);
hold on;
if k == 100
plot(x,mass_flowrate,\'g\',\'linewidth\',1.25);
else if k == 200
plot(x,mass_flowrate,\'r\',\'linewidth\',1.25);
else if k == 700
plot(x,mass_flowrate,\'b\',\'linewidth\',1.25);
legend(\'100\\Deltat\',\'200\\Deltat\',\'700\\Deltat\');
end
end
end
title(\'Comparison of non-dimensional mass flow rates at different time-instants (Conservative Form)\',\'linewidth\',3);
xlabel(\'x/L (Non-dimensional distance through the nozzle)\');
ylabel(\'rhoVA/rho_{o}V_{o}A_{o}\');
axis([0 3 0 0.7]);
grid on;
grid minor;
end
end
Results:
For Grid points-31
1. Non-Dimensional mass flowrate distribution at various time-steps:
a) Non-conservative form:
b) Conservative Form:
The above plots are used to determine the minimum number of iterations needed for the convergence. The conservative form has less variation and the final mass-flow rate is very close to the analytical solution 0.579 when compared to the non-conservative form. But we also have to take into account that the better initial condition provided to the conservative form may be the reason for this better accuracy and we cannot simply drive the conclusion that the conservative form is better.
2. Time-wise variation of pressure, temperature, velocity and Mach number at the nozzle throat:
a) Non-Conservative Form:
b)Conservative Form:
The above plots depict the properties at throat during the solution. The non-dimensional values of the properties at the throat are well-known and that is why this section was chosen to check these calculations.
We can also check the stability of the two forms. Here the conservative form seems to be more stable. There are irregularities at the beginning in the non-conservative form but this dies down as the solution marches forward in time.
3. Qualitative aspects of Quasi 1D nozzle flow at steady state:
a) Non-Conservative form:
b) Conservative form:
The quantities in the above plots follow the trend that the nozzle flow problem should follow. The density, temperature, and pressure drops while moving towards the exit and the Mach number increases towards the end. These values match the analytical solution.
4. Comparison of the non-dimensional mass flowrates of both the forms with the exact solution:
The steady mass-flow results obtained with the conservative form are, on the whole, much closer to the exact solution(0.579). The conservation form does a better job of preserving mass throughout the flow-field, mainly because the mass-flow itself is one of the dependent variables in the equation, hence mass-flow is a primary result from these equations. In contrast, the dependent variables in the non-conservative form of the equation are primitive variables(density, temperature, velocity) and the mass-flow is obtained only as a secondary result. As the conservation form of the equation does a better job of conserving mass throughout the flow-field hence the name conservative form. But this does not establish the superiority of the conservative form over the non-conservative form which will be clear in the latter discussion.
Time and number of iterations taken by both the forms(for 31 grid points):
For Conservative form:
For Non-conservative form:
5.Grid Dependency Test:
A grid dependency test was ran by increasing the number of grid points from 31 to 61 and the results were checked for both forms of the equation and there was a little change in by doubling the number of grid points hence we have gird independence at 31 grid points.
Values are checked for primitive variables at the throat as we know values at the throat.
In the previous section where the comparison between the mass-flow rates of both forms are compared with the exact solution is made, a statement was made that though the mass is better conserved in the conservative form it does not establish the superiority of the conservative form over the non-conservative form and this is evident from the table above showing the values of primitive variables along with the Mach number at the throat. As it is clear that here the values are much closer to the exact solution for the Non-conservative form than for Conservative form and by doubling the number of grid points there is a very minor change in the quantities hence the grid independence.
Results for grid points - 61
1. Non-Dimensional mass flowrate distribution at various time-steps:
a) Non-conservative form:
b) Conservative Form:
2. Time-wise variation of pressure, temperature, velocity and Mach number at the nozzle throat:
a) Non-Conservative Form:
b)Conservative Form:
3. Qualitative aspects of Quasi 1D nozzle flow at steady state:
a) Non-Conservative form:
b)Conservative Form:
Time and number of iterations taken by both the forms(for 61 grid points):
Conservative form:
Non-conservative form: