Week 7 - Simulation of a 1D Supersonic nozzle flow simulation using Macormack Method
INTRODUCTION TO THE GIVEN PROBLEM
We are examining a continuous and steady isentropic flow through a converging-diverging nozzle, as illustrated in the diagram below.
At the nozzle's inlet, the flow originates from a reservoir with a theoretically large reservoir pressure (P∞ → ∞) and, consequently, a very low velocity (V → 0). As a result, we refer to these conditions as the stagnation values, denoted as P₀ and T₀, representing the total pressure and total temperature, respectively. The flow then undergoes an isentropic expansion, reaching supersonic velocities at the nozzle exit, characterized by exit pressure and temperature denoted with subscripts 'e.'
Within the nozzle, the flow experiences different local speeds: it is subsonic in the convergent section, reaches the sonic speed at the throat (where the cross-sectional area is at a minimum), and becomes supersonic in the divergent section. When we mention 'sonic flow' (Mach number M = 1) at the throat, it means that the local velocity matches the local speed of sound.
We can use our own notation, such as '*', to represent the throat values.
Several useful results have been derived from this analysis."
Mach number:
Mach number is a dimensionless quantity in fluid dynamics representing the ratio of flow velocity past a boundry to the local speed of sound
We can express the velocity of gas (u) as a function of the sonic velocity
M= u/c
- Where M is the Mach Number;
- When M < 1>
- M = 1 Sonic
- M > 1 Supersonic
- c = Sonic velocity (m/s)
- γ = Specific heat ratio
- R = Specific gas constant ( J / mol . K)
- T = Temperature of the fluid (K)
Set Up
MaCcormack mathod:
In the realm of computational fluid dynamics, the MacCormack method stands as a widely embraced discretization technique for numerically solving hyperbolic partial differential equations.
The MacCormack method achieves second-order accuracy without the need for second-order derivatives, thereby simplifying the computational intricacies. It is categorized as an explicit method.
Within the framework of the MacCormack method, the determination of a variable's value at the grid point (i, j) at the time (t + dt) follows this prescription:
Governing Equations:
The supersonic flow through a nozzle is governed by the Navier-Stokes equations, which can be derived for both conservative and non-conservative formulations.
The continuity equation, expressing the conservation of mass, can be stated as follows:
The symbols T₀ and ρ₀ represent the reservoir temperature and density, respectively. We establish the non-dimensional counterparts for temperature and density as follows
The shape of the nozzle, denoted as A = A(x), is specified and remains fixed over time. In this particular case, we adopt a parabolic area distribution:
A = 1 + 2.2 * (x - 1.5)^2 for 0 ≤ x ≤ 3.
For the density, temperature, and velocity, we have the following expressions:
ρ = 1 - 0.3146 * x
T = 1 - 0.2314 * x
V = (0.1 + 1.09 * x) * √T
Pressure (P) is determined by the product of density (ρ) and temperature (T), and the Mach number (M) is calculated as the ratio of velocity (u) to the square root of temperature (√T).
These initial conditions apply at t = 0.
Now, let's move on to the predictor step:
In this phase, we discretize the problem using forward difference equations.
To achieve this, we establish finite difference expressions for spatial derivatives. Additionally, for the sake of simplicity, we will omit the prime symbol to represent dimensionless variables."
And some of the other equations including for other variables like momentum and energy as well .
For grid generation we divide the length of the nozzle in number of points along the axis .The characteristics of the flow changes from subsonic to super sonic from right to left .The distance between the two grid points is dx.Smaller the mesh size the error also decreases.
In discretisation part we apply thr MaCcormach technique to solve the the conservative and non conservative forms of governing equations.
From these we get the predicted equation as follows: These are denoted by barred quantities.
In a similar approach, the conservative form of the governing equation can be solved using MacCormack's method with predictor and corrector.
SOLUTION:
Inorder for the code to be modular, the code is split into three parts.
Main program
Function for Conservative form
Function for Non-conservative form
Main Program:
close all
clear all
clc
%Length of nozzle
L=3;
%Number of grid points
n=[31 62 93 124 155];
%Number of time steps
nt=1:1400;
�L number
cfl = 0.5;
for i = 1:length(n)
x=linspace(0,L,n(i));
[mf_nc,mf_throat_nc]=non_conservative1(L,n(i),nt,cfl);
[mf_c,mf_throat_c]=conservative1(L,n(i),nt,cfl);
figure(7)
subplot(2,1,1)
plot(nt,mf_throat_c,'b',nt,mf_throat_nc,'r','linewidth',1)
xlabel('Time steps (nt)')
ylabel('Mass flow rate')
legend('Conservative form','Non-Conservative form')
title('Mass flow rate using "Non-Conservative form vs Conservative form" at throat point i=16')
subplot (2,1,2)
plot(x,mf_c,'b',x,mf_nc,'r','linewidth',1)
xlabel('Length of nozzle')
ylabel('Mass flow rate')
legend('Conservative form','Non-Conservative form')
title('Mass flow rate using "Non-Conservative form vs Conservative form" in the nozzle')
figure(8)
subplot (2,1,1)
plot(nt,mf_throat_c,'linewidth',1)
hold on
xlabel('Time steps (nt)')
ylabel('Mass flow rate')
legend('31 nodes','62 nodes','93 nodes')
title('Conservative form: Grid independency test at throat point i=16')
subplot (2,1,2)
plot(nt,mf_throat_nc,'linewidth',1)
hold on
xlabel('Time steps (nt)')
ylabel('Mass flow rate')
legend('31 nodes','62 nodes','93 nodes','124 nodes','155 nodes')
title('Non-Conservative form: Grid independency test at throat point i=16')
figure(9)
subplot (2,1,1)
plot(x,mf_c,'linewidth',1)
hold on
xlabel('Length of nozzle')
ylabel('Mass flow rate')
legend('31 nodes','62 nodes','93 nodes','124 nodes','155 nodes')
title('Conservative form: Grid independency test along the nozzle')
subplot (2,1,2)
plot(x,mf_nc,'linewidth',1)
hold on
xlabel('Length of nozzle')
ylabel('Mass flow rate')
legend('31 nodes','62 nodes','93 nodes','124 nodes','155 nodes')
title('Non-Conservative form: Grid independency test along the nozzle')
end
Explanation:
Input parameters were given for the variables x, dx, n, gamma, nt and cf1.
Non Conservative and Conservative forms are created as separate functions with the neccessary inputs.
Plotting is done for the following solutions
1)Timestep variation at throat area of both forms
2)Variation of Flow Rate Distribution of both forms
3)Comparison of Normalized Mass Flow Rate of both forms
Each plot is computed as subplots of four primitive variables such as Density, Temperature, Pressure and Mach Number at throat as well as Convergent/Divergent section.
Function for Non conservation form
% FUNCTION FOR NON-CONSERVATIVE FORM OF GOVERNING EQUATION
function [mf_nc,mf_throat_nc]=non_conservative1(L,n,nt,cfl)
L=3;
x=linspace(0,L,n);
dx=x(2)-x(1);
% initial conditions rho, a, v, T, mf, M, p
rho = 1-0.3146*x;
T= 1-0.2314*x;
v=(0.1+1.09*x).*T.^0.5;
a= 1+2.2*(x-1.5).^2;
M= v./(T.^0.5);
mf = rho.*a.*v;
p=rho.*T;
g = 1.4;
for z=1:31
u=(T(z))^(1/2);
dt(z) = cfl*dx/(u + v(z));
end
dt=min(dt);
for i = 1:length(nt)
%Saving a copy of old values
rho_old = rho;
T_old = T;
v_old = v;
%Predictor step
for j=2:n-1
drho_dx = (rho(j+1)-rho(j))/dx;
dv_dx = (v(j+1)-v(j))/dx;
dT_dx = (T(j+1)-T(j))/dx;
da_dx = (log(a(j+1))-log(a(j)))/dx;
%Continuty
drho_dt_p(j) = -rho(j)*dv_dx - rho(j)*v(j)*da_dx - v(j)*drho_dx;
%Momentum
dv_dt_p(j) = -v(j)*dv_dx - (1/g)*(dT_dx + (T(j)/rho(j))*drho_dx);
%Energy
dT_dt_p(j) = -v(j)*dT_dx - (g-1)*T(j)*(dv_dx + v(j)*da_dx);
%Updating primitives
v(j) = v_old(j) + dv_dt_p(j)*dt;
rho(j) = rho_old(j) + drho_dt_p(j)*dt;
T(j) = T_old(j) + dT_dt_p(j)*dt;
%Corrector step
drho_dx = (rho(j)-rho(j-1))/dx;
dv_dx = (v(j)-v(j-1))/dx;
dT_dx = (T(j)-T(j-1))/dx;
da_dx = (log(a(j))-log(a(j-1)))/dx;
%Continuity
drho_dt_c(j) = -rho(j)*dv_dx - rho(j)*v(j)*da_dx - v(j)*drho_dx;
%Momentum
dv_dt_c(j) = -v(j)*dv_dx - (1/g)*(dT_dx + (T(j)/rho(j))*drho_dx);
%Energy
dT_dt_c(j) = -v(j)*dT_dx - (g-1)*T(j)*(dv_dx + v(j)*da_dx);
end
%Average time derivatives
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 corected values
for k = 2:n-1
rho(k) = rho_old(k) + drho_dt_avg(k)*dt;
v(k) = v_old(k) + dv_dt_avg(k)*dt;
T(k) = T_old(k) + dT_dt_avg(k)*dt;
end
%Inlet boundary conditions
v(1) = 2*v(2)-v(3);
%Exit 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);
%Final update of primitives
p = rho.*T;
M= v./T.^0.5;
error(k) = abs(v(k)-v(k-1));
mf= rho.*a.*v;
%Finding values at throat point
T_throat(i) = T(16);
rho_throat(i) = rho(16);
v_throat(i) = v(16);
M_throat(i) = M(16);
p_throat(i) = p(16);
error_throat(i)=error(16);
mf_throat(i)=mf(16);
if n==31
for ii=[1:50:201 700]
if ii==i
figure(1)
plot(x,mf,'linewidth',1)
hold on
title ('Non-Conservative form - Variation of Mass flow rate distribution inside the nozzle at different time steps during the time-marching process')
xlabel('Length of nozzle')
ylabel('Mass flow rate')
elseif i==1400
plot(x,mf,'*r')
legend('1','51','101','151','201','700','1400')
end
end
end
end
if n==31
figure(2)
plot(x,v,'b',x,p,'r',x,T,'k',x,rho,'m',x,M,'g',x,mf,'c')
ylim([-0.1 3.7])
title (sprintf('Time step = %d',i))
xlabel('Length of nozzle')
ylabel('Non dimensional parameters')
legend('Velocity','Pressure','Temperature','Density','Mach number','Mass flow rate','location','northwest')
title('Non-Conservative form - Steady-state distribution of primitive variables inside the nozzle')
figure (3)
plot(nt,rho_throat)
hold on
plot(nt,v_throat,'linewidth',1)
plot(nt,p_throat,'linewidth',1)
plot(nt,T_throat,'linewidth',1)
plot(nt,M_throat,'linewidth',1)
plot(nt,mf_throat,'linewidth',1)
xlabel('Time steps (nt)')
ylabel('Non dimensional parameters')
legend('Density','Velocity','Pressure','Temperature','Mach number','Mass flow rate')
title('Non-Conservative form - Time-wise variation of the primitive variables at throat point i=16')
end
mf_nc=mf;
mf_throat_nc=mf_throat;
end
Explanation:
The purpose of the function is to calculate the values of
Mass flow rate, Pressure, Mach number, Density, velocity, Temperature, Density at throat, velocity at throat, Temperature at throat, Mass flow rate at throat, Pressure at throat, Mach number at throat by the known equations for all the grid points in the non conservative form of the equation
For this initial profiles for the Density, Temperature, Velocity and area of the nozzle has been defined.
Outer Time loop is executed, to calculate all the grid points for various formulas inside the loop.
MacCarmack technique is used to solve the numerical solution of governing equations such as Continuity, Momentum & Energy Equation.
Predictor, Corrector method is employed to solve the equations
Predictor & Corrector method are assigned within the for loop individually for calculating values.
The calculated values are averaged in the next step and updating the solution with old values.
Boundary Conditions are applied within the for loop and subsequently used for calculating the Mass Flow Rate, Mach Number & Pressure and also calculating the values throat section.
Finally plotting is done for Distance(x) & Mass Flow Rate for differrent Timesteps.
Function for Conservative form:
% FUNCTION FOR CONSERVATIVE FORM OF GOVERNING EQUATION
function [mf_c,mf_throat_c]=conservative1(L,n,nt,cfl)
L=3;
x=linspace(0,L,n);
dx=x(2)-x(1);
% initial conditions rho, a, v, T, M, mf, p
for i=1:length(x)
if (x(i)<=0.5)
rho(i) = 1;
T(i) = 1;
elseif (x(i)>0.5) & (x(i)<=1.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
end
a= 1+2.2*(x-1.5).^2; %area
g = 1.4; %gamma value
% solution vectors u1, u2, u3
u1 = rho.*a;
u2 = 0.59*ones(1,n);
v=u2./u1;
u3 = u1.*((T./(g-1))+(g/2).*v.^2);
M= v./(T.^0.5); %Mach number
mf = rho.*a.*v; %Mass flow rate
p=rho.*T; %Pressure
%Time step control
for z=1:31
u=(T(z))^(1/2);
dt(z) = cfl*dx/(u + v(z));
end
dt=min(dt);
for i = 1:length(nt)
%Saving a copy of old values
u1_old = u1;
u2_old = u2;
u3_old = u3;
%%Predictor Step
for j=2:n-1
�lculating flux vectors
f1 = u2;
f2 = ((u2.^2)./u1)+((g-1)/g)*(u3-((g/2)*(u2.^2)./u1));
f3 = ((g.*u2.*u3)./u1) - ((g*(g-1)/2)*(u2.^3)./(u1.^2));
%Source Term
j1 = (1/g)*rho(j)*T(j)*((a(j+1)-a(j))/dx);
%Continutiy Term
du1_dt_p(j)= -(f1(j+1)-f1(j))/dx;
%Momentum Term
du2_dt_p(j)= -(f2(j+1)-f2(j))/dx + j1;
%Energy Term
du3_dt_p(j)= -(f3(j+1)-f3(j))/dx;
%Intermediate solution update
u1(j) = u1_old(j) + du1_dt_p(j)*dt;
u2(j) = u2_old(j) + du2_dt_p(j)*dt;
u3(j) = u3_old(j) + du3_dt_p(j)*dt;
�coding the predicted values of primitive variable from solution vectors U1,U2,U3
rho=u1./a;
T=(g-1).*(u3./u1-(g/2).*(u2./u1).^2);
%Updating flux terms
f1=u2;
f2=u2.^2./u1+((g-1)/g).*(u3-(g/2)*u2.^2./u1);
f3=g.*u2.*u3./u1-0.5*g*(g-1).*u2.^3./u1.^2;
%Source term for corrector step
j2 = (1/g)*(rho(j)*T(j)*(a(j)-a(j-1))/dx);
%Corrector step continutiy Term
du1_dt_c(j) = -(f1(j)-f1(j-1))/dx;
%Corrector step momentum Term
du2_dt_c(j) = -(f2(j)-f2(j-1))/dx + j2;
%Corrector step energy Term
du3_dt_c(j) = -(f3(j)-f3(j-1))/dx;
end
%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 k = 2:n-1
u1(k) = u1_old(k) + du1_dt_avg(k)*dt;
u2(k) = u2_old(k) + du2_dt_avg(k)*dt;
u3(k) = u3_old(k) + du3_dt_avg(k)*dt;
end
%Inlet boundary conditions
u1(1) = a(1);
u2(1) = 2*u2(2)-u2(3);
u3(1) = u1(1)*(T(1)/(g-1) + (g*(u2(1)/u1(1))^2)/2);
%Exit boundary 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 update of primitives
rho = u1./a;
v = u2./u1;
T=(g-1).*(u3./u1-(g/2).*(u2./u1).^2);
p = rho.*T;
M= v./T.^0.5;
error(k) = abs(v(k)-v(k-1));
mf=rho.*a.*v;
%Finding values at throat point
T_throat(i) = T(16);
rho_throat(i) = rho(16);
v_throat(i) = v(16);
M_throat(i) = M(16);
p_throat(i) = p(16);
error_throat(i)=error(16);
mf_throat(i)=mf(16);
if n==31
for ii=[1:50:201 700]
if ii==i
figure(4)
plot(x,mf,'linewidth',1)
hold on
title ('Conservative form - Variation of Mass flow rate distribution inside the nozzle at different time steps during the time-marching process')
xlabel('Length of nozzle')
ylabel('Mass flow rate')
elseif i==1400
plot(x,mf,'*r')
legend('1','51','101','151','201','700','1400')
end
end
end
end
if n==31
figure(5)
plot(x,v,'b',x,p,'r',x,T,'k',x,rho,'m',x,M,'g',x,mf,'c')
ylim([-0.1 3.7])
title (sprintf('Time step = %d',i))
xlabel('Length of nozzle')
ylabel('Non dimensional parameters')
legend('Velocity','Pressure','Temperature','Density','Mach number','Mass flow rate','location','northwest')
title('Conservative form - Steady-state distribution of primitive variables inside the nozzle')
figure (6)
plot(nt,rho_throat)
hold on
plot(nt,v_throat,'linewidth',1)
plot(nt,p_throat,'linewidth',1)
plot(nt,T_throat,'linewidth',1)
plot(nt,M_throat,'linewidth',1)
plot(nt,mf_throat,'linewidth',1)
xlabel('Time steps (nt)')
ylabel('Non dimensional parameters')
legend('Density','Velocity','Pressure','Temperature','Mach number','Mass flow rate')
title('Conservative form - Time-wise variation of the primitive variables at throat point i=16')
end
mf_c=mf;
mf_throat_c=mf_throat;
end
Explanation:
The purpose of the function is to calculate the values of
Mass flow rate, Pressure, Mach number, Density, velocity, Temperature, Density at throat, velocity at throat, Temperature at throat, Mass flow rate at throat, Pressure at throat, Mach number at throat by the known equations for all the grid points in the conservative form of the equation
Initially, the profiles for the density, temperature has been defined in accordance with distance along the nozzle,piece wise values has been provided for the initial profiles
The Conservative form of the equation is split into three parts namely, solution vector (u), flux vector (f) and source term (j).
Outer Time loop is executed, to calculate all the grid points for various formulas inside the loop.
MacCarmack technique is used to solve the numerical solution of governing equations such as Continuity, Momentum & Energy Equation.
Predictor, Corrector method is employed to solve the equations
Predictor & Corrector method are assigned within the for loop individually for calculating values.
The calculated values are averaged in the next step and updating the solution with old values.
Boundary Conditions are applied within the for loop and subsequently used for calculating the Mass Flow Rate, Mach Number & Pressure and also calculating the values throat section.
Finally plotting is done for Distance(x) & Mass Flow Rate for differrent Timesteps
PLOTS:
(a) Non- Conservative form
Fig.3: Variation of mass flow rate distribution inside the nozzle at different time steps
Fig.4:Steady state distribution of primitive varibles inside the nozzle
Fig.5: Steady state variation of primitive variables at the throat of the nozzle
(b) Conservative form
Fig.6: Variation of mass flow rate distribution inside nozzle at different time steps
Fig.7: Steady state distribution of primitive variables inside the nozzle
Fig.8: Temporal variation of primitive variable at throat of the nozzle
(c) Conservative Vs Non-Conservative form
Fig.9: Comparison of mass flow rate for conservative and non conservative form in nozzle and at the throat of the nozzle
(d) Grid Independence study
Fig.10: Grid independence test at the throat of the nozzle
Grid independence test has been carried out for the both type of solution methodologies, i.e., conservative and non-conservative. Since the problem is 1 dimensional, only we are concerned with the number of grid points. The grid points has been increased from 31 to 155 and mass flow rate has been used as the variable to determine the grid independence. We could see in both cases as shown in figure 10&11 that, as the grid points are being increased, there is a very slight change in the mass flow rate.But taking the average of the mass flow rate along the length of the nozzle, we can see that, increasing the grid points further beyond 31 does not have a significant effect on the solution, which shows that the solution has became independent of the grid.Thus discretization using 31 grid points is optimum for the given quasi one-dimensional problem.
OBSERVATIONS:
For the non-conservative case, we can see that, during the initial time step, i.e., as the flow begins through the entry of the nozzle,there is a lot of fluctuations in the solution.But as the time progress, the solution becomes steady.But in contrast for the conservative case, the flow do not suffer much of the oscillation as observed in non-conservative case, this can be seen in figure 3&6.From the figure 9,we are ble to see that the solution reaches a steady state faster for the conservative case and also during the initial stage of the solution, the solution suffers oscillations.
The plot of the variation of the primitive variables along the length of the nozzle as in figure 4 &7 shows that, at the entry of the nozzle the velocity and mach number of the flow is low, indicating that the flow is subsonic, but this trend continues till the throat region of the nozzle, and there is a gradual increase in the mach number and velocity.But afterwards, the flow enters the supersonic regime and the magnitude of velocity and mach number increases.Subsequently,the density and pressure reduces in this regime.A similar trend is observed in the case of primitive varibles for both conservativve and non conservative form.
Primitive variables at the throat of the nozzle suffers large oscillations in the solution for non conservative case, but the oscillation is dampened out to an extend for the conservative case. The time taken for the solution to become steady is almost the same in both the cases.
CONCLUSIONS:
After examining the above plots and the results obtained, we can evidently conclude that, the rate of convergence is more in the case of conservative type of the solution and the solution is stable in every part of the domain.Whereas the non-conservative form experiences, oscillations to the solutions during the initial stages, but it dampens out gradually.Yet both the forms are dependable in calculating the solutions which is almost precise.