Week 7 - Simulation of a 1D Supersonic nozzle flow simulation using Macormack Method
Objective :
To simulate the isentropic flow through a quasi 1D subsonic-supersonic nozzle by deriving both conservative and non conservative forms of governing equation and solve them by using Macormack Method in matlab.
Problem statement :
- Need to determine the steady-state temperature distribution for the flow-field variables and to investigate the difference between the two forms of governing equations by comparing their solutions.
- To plot steady-state distribution of primitive variables inside the nozzle.
- To plot time-wise variation of the primitive variables.
- To plot variation of Mass flow rate distribution inside the nozzle at different time steps during the time-marching process.
- To plot comparison of Normalized mass flow rate distributions of both forms.
Theory :
Mach number :-
It is a dimensionless quantity representing the ratio of flow velocity to the local speed of the sound.
`M = u/c`
Where,
M - Mach Number
u - Flow Velocity
c - Velocity Sound
Nozzle diagram :
The nozzle has three main sections
1. convergent.
2. Throat.
3. Divergent
- At convergent section the flow velocity is less than the sound velocity therefore there mach number is less than 1.
- At throat section the flow velocity is equal to the sound velocity therefore there mach number is equal to 1.
- At divergent section the flow velocity is greater than the sound velocity therefore there mach number is greater than 1.
So here, the flow at the inlet of the nozzle comes from a reservoir where pressure and temperature are denoted as by `rho_0` and `T_0` respectively. The cross sectional area of the reservoir is large (theoretically, A→∞">→∞) and hence the velocity is very small (V→">→0). Thus, `rho_0` and `T_0` are the stagnation values or total pressure and temperature respectively. The flow expands isentropically to supersonic speeds at the nozzle exits, where the exit pressure, temperature, velocity, and Mach number are denoted by `P_e`, `T_e`, `V_e` and `M_e` respectively. The flow is locally subsonic in the convergent section, sonic at the throat section and supersonic at the divergent section. The sonic flow (M=1) at the throat means that the local velocity at this location is equal to the local speed of sound. Using an asterisk to denote sonic flow values, we have at the throat V = `V^**` = `a^**` Similarly, the sonic flow values of pressure and temperature are denoted by `P^**` and `T^**`, respectively.
Although the area of nozzle changes as a function of distance along the nozzle x and therefore in reality the flow field is two dimensional (xy), we make an assumption that the flow properties vary only with x, this tendency to assume uniform flow properties across any given cross-section. Such flow is defined as Quasi-one-dimensional flow.
In order to solve quasi-one-dimensional problem, we need PDE or Integral form of equations for CFD analysis. These PDE or Integral form of the governing equations are in both conservative and non-conservative form.
Before, writing the governing equations in conservative and non-conservaticve form, we need to non-dimensionalize the primitive variables for removing the constraints of dimensions to analyse the numerical solution. It will certainly ease the method of analysis for a given problem.
Non Dimensional Quantities :
Also as part of the problem setup we represent quantities in terms of non-dimensional quantities
Temperature `T' = T/T_0`
Density `rho' = rho/rho_0`
Non - dimensional distance `X' = X/L`
Non dimensional velocity `V' = V/a_0` where as `a_0` is the local sonic velocity (`a_0 = sqrt (gamma RT_0)`)(a0=γRT0)">
Non-dimensional time `t' = t/((L/a_0))` where `L/a_0` is characteristic time scale.
Non dimensional area `A' = A/A^**`
Non Conservative form :
MacCormack method :
MacCormack's method is an explicit finite-difference technique which is second order accurate in both space and time. MacCormack method generally uses three step approach(Predictor, Corrector and averaging) to calculate a final numerical solution for an iteration. In predictor step, forward difference is taken for spatial derivatives whereas in corrector step, rearward (backward) difference is taken for spatial derivatives. Furthermore, the solutions obtained from these two steps are averaged and final numerical solution is calculated for a single iteration.
Initial and boundary conditions :
To start the time-marching calculations, we must stipulate initial conditions for ρ">ρρ, T and V as a function of x; that is, we must set up values of ρ">ρρ, T and V at time t = 0.
Non-conservative form :
Initial conditions :
`rho = 1 − 0.3146 x`
`T = 1 − 0.2314 x`
`V = (0.1 + 1.09 x) . T ^(1/2)`
Boundary conditions :
In our case, we are having a subsonic flow in the convergent section (M<1), choked flow at the throat section (M=1) and supersonic outflow in the divergent section (M>1). So by using method of characteristics for an unsteady, one-dimensional flow we can precisely allot inflow and outflow boundary conditions and can identify which quantities to state at the boundaries.
Fig (b) : Study of boundary conditions: subsonic inflow and supersonic outflow
What does these characteristics curve and streamline direction tells us?
- If any characteristic line is entering into the domain, implies that the value of one dependent flow field variable must be specified at that boundary.
- However, if any characteristic line is exiting the domain, implies that the value of the another dependent flow field variable must be allowed to float at that boundary i.e. it must be calculated in steps of time as function time wise solution of the flow field.
- Similarly, if streamline flows into the domain then specify third dependent flow field variable and if streamline flows out of the domain then third dependent flow field variable must be allowed to float.
In our case, at i = 1(inlet), we have primitive varibles like `rho`, A, V and T. So which variables must be taken for initial and boundary conditions? Here, we know the ρ">ρρ and T at node 1 are equal to the reservoir's ρ">ρρ and T which is assumed to be `rho = 1` and `T = 1`, so it becomes,
Inflow boundary condition :
ρ(1)=1"> `rho(1) = 1`and `T(1) = 1`(Since it is Dirichlet boundary condition)
However, velocity and area tends to infinity, so for velocity, we need to calculate it from the information obtained from inner nodes,
V(1)=2⋅V(2)-V(3)">V(1) = 2V(2) - V(3) (It is Neumann boundary condition)
Outflow boundary condition :
At the outflow boundary condition (Node =n), as we can see that, both the characteristics lines are out of the domain and streamline also flows out of the domain. This is because, at the outflow condition, the flow is supersonic (M>1), which is certainly greater than the speed of sound. So in this condition, information transport is very fast. Therefore, there are no flow field variables which require their values to be stipulated at the supersonic outflow boundary; all variables must be allowed to float at this boundary.
ρ"> `rho (n) = 2rho(n-1) - rho (n-2)`
`V (n) = 2V (n-1) - V (n-2)`
`T(n) = 2T (n-1) - T(n-2)`
The supersonic outflow boundary is purely Neumann boundary condition, which is extrapolating the boundary nodes using the information from inner nodes.
Conservative form :
In conservative form, we will initialize the nozzle profile rationally. As in conservative form, seen above from the governing equations, we have solution vectors, flux terms and a source term. So it is the prior matter of concern to state initial and boundary conditions in those terms rather than that of primitive variables.
Initial conditions :
These variations are slightly more realistic than those assumed in the nozzle shape and initial conditions in non-conservative form of governing equations discussed above. This is an anticipation that the stability behaviour of the finite-difference formulation using the conservation form of the governing equations might be slightly more sensitive, and therefore it is useful to start with more improved initial conditions than those discussed above in non-conservative form.
Assumptions for initial profile of the nozzle :
- For `0≤x'≤0.5≤x′≤0.5`
`rho' = 1`
`T' = 1`
- For `0≤x'≤0.5≤x′≤0.5`
`rho'= 1.0 - 0.366(x'-0.5)`
`T' = 1.0 - 0.167(x'-0.5)`
- For `0.5≤x'≤1.5≤x′≤1.5`
`rho' = 0.634 - 0.3879(x'-1.5)`
`T' = 0.833 - 0.3507(x'-1.5)`
Before stating the initial conditions in terms of solution vectors, we need to initialize velocity(V'). We have one of the dependent variable in our governing equation as `U_2`, which is physically local mass flow; that is `U_2 = rho' A' V'`. Therefore, for the initial conditions only; let us assume a constant mass flow through the nozzle and calculate V' as,
`V' = (U_2)/ (rho' A') = (0.59)/ (rho' A')`
Note: The value 0.59 is chosen for `U_2`, because it is close to exact analytical value of the steady-state mass flow (which is for this case is 0.579)
Therefore, the initial condition for V' as a function of x' is obtained by substituting the ρ′">`rho'` and A' variation. Finally, the initial conditions for `U_1`,`U_2` and `U_3` are obtained by substituting the above variations for `rho' T'` and `V'` into the definitions discusssed above.
`U_1 = rho' A'`
`U_2 = rho' A' V'`
`U_3 = rho ' ((e ')/ (gamma − 1) + gamma/2 (V ') ^2) A '`
where, e' = T'. Certainly, for the initial conditions described above, V' is calculated such that `U_2 = rho' A' V' = 0.59`.
Boundary conditions :
In conservative form, by referring fig (b) and absorbing the same method of characteristics approach we can allot boundary conditions for inflow and outflow boundary conditions. However, in this case, the boundary conditions will be in terms of solution vectors. Similarly as discussed above, at the subsonic inflow (M<1), we need to specify two quantities at the initial node, where we have `rho' ,A'` and `U_2 = rho' A' V' = 0.59`(Mass flow rate) specified and one quantity will be allowed to float by referring left hand characteristic line. However, at the supersonic outflow (M>1), the information is travelling exceeds the speed of sound, hence it is irrational to specify outflow boundary condition. Therefore, all the dependent variables are allowed to float and they will be calculated in steps of time as a function time-wise solution of the flow field.
Inflow boundary condition :
The inflow boundary condition can be given as,
`U_1 = rho' (1) A ' (1)` (Dirichlet boundary condition)
`U_2 = 2.U_2 (2) − U_2 (3)`
`V (1) = (U _2(1))/ (U_1 (1))`
Note: This velocity at the first node will be substituted below in calculating U3
`U_3 = U_1 (1) ' ((T ' (1))/ (gamma − 1) + gamma/2 (V ' (1))^2)` (Neumann boundary condition)
Outflow boundary condition :
All the dependent variables will be allowed to float. Hence, the outflow boundary can be given as,
`U_1(n) = 2.U_1(n-1) - U_1(n-2)` (Neumann boundary condition)
`U_2 (n) = 2.U_2(n-1) - U_2 (n-2)`(Neumann boundary condition)
`U_3 = 2.U_3(n-1) - U_3(n-2)` (Neumann boundary condition)
Time step calculation :
The governing system of equations is hyperbolic with respect to time. So there is stability constraint on this system which is known as Courant-Friedrich-Lewy (CFL) condition. The Courant number is a simple stability analysis of a linear hyperbolic equations that gives result CFL≤1">≤1 for an explicit numerical solution to be stable.
CFL condition can be given as,
CFL condition = (any quantity)*(time-step with which information would travel)/(element step size)
So, in our case, the quantity which is going to travel is velocity, or in other sense velocity is the quantity which will be travelling with certain time from one node to other. So here velocity of the fluid at a point in the flow is coupled with local speed of sound. The quantity can be expressed to be as a+v'
CFL = `("a + v'). (Delta t)/(Delta x)`
We will assume the Courant number as C = 0.5, accordingly will calculate the `Delta t`. In addition to it, we will calculate `(Delta t) _i ^t`
at all the grid points, i = 1 to i = n, and then we'll choose the minimum value of `Delta t` in our calculations.
After considering all these governing equations and their respective pre-requisites to solve the problem, then we will move towards to the MATLAB code for execution.
Main program :
%Quasi 1D Supersonic Nozzle Flow using Macormack method
clear all
close all
clc
n = 61; %Number of Nodes
x = linspace(0, 3, n); %Domain length
gamma = 1.4;
method = input('Enter (1) for Non-Conservative Form and (2) for Conservative Form: n');
if method == 1
[x, v, T, rho, M, P, mdot, M_thrt, rho_thrt, T_thrt, P_thrt, mdot_thrt, iteration_no, residual1, residual2, residual3] = NCF(n, x, gamma);
%Plotting
%figure(1)
%plot(x, mdot_thrt(different timesteps))
figure(2)
%plot(x, v, x, T, x, rho, x, M, x, P, LineWidth=2)
%legend('Velocity', 'Temperature Ratio', 'Density Ratio', 'Mach Number', 'Pressure Ratio')
subplot(5, 1, 1)
plot(x, M, LineWidth=2, Color='r')
ylabel('Mach Number')
grid on
title('Steady State Distribution of variables inside the nozzle - Non-Conservative Form')
subplot(5, 1, 2)
plot(x, P, LineWidth=2, Color='g')
ylabel('Pressure Ratio')
grid on
subplot(5, 1, 3)
plot(x, T, LineWidth=2, Color='b')
ylabel('Temperature Ratio')
grid on
subplot(5, 1, 4)
plot(x, rho, LineWidth=2, Color='y')
ylabel('Density Ratio')
grid on
subplot(5, 1, 5)
plot(x, mdot, LineWidth=2, Color='m')
xlabel('Nozzle Length')
ylabel('Mass Flow Rate')
grid on
figure(3)
plot(iteration_no, M_thrt, iteration_no, rho_thrt, iteration_no, T_thrt, iteration_no, P_thrt, iteration_no, mdot_thrt, LineWidth=2)
grid on
xlabel('Number of Timesteps')
legend('Mach Number', 'Density Ratio', 'Temperature Ratio', 'Pressure Ratio', 'Mass flow rate')
title('Time-wise variation of variable at nozzle throat - Non-Conservative Form')
figure(4)
plot(iteration_no, residual1, iteration_no, residual2, iteration_no, residual3, LineWidth=2)
grid on
xlabel('Number of Timesteps')
ylabel('Residuals')
legend('Density Ratio Residuals', 'Velocity Ratio Residuals', 'Temperature Ratio Residuals')
title('Residual Values - Non-Conservative Form')
elseif method == 2
[x, v, T, rho, M, P, mdot, M_thrt, rho_thrt, T_thrt, P_thrt, mdot_thrt, iteration_no, residual1, residual2, residual3] = CF(n, x, gamma);
%Plotting
%figure(1)
%plot(x, mdot_thrt(different timesteps))
figure(2)
%plot(x, v, x, T, x, rho, x, M, x, P, LineWidth=2)
%legend('Velocity', 'Temperature Ratio', 'Density Ratio', 'Mach Number', 'Pressure Ratio')
subplot(5, 1, 1)
plot(x, M, LineWidth=2, Color='r')
ylabel('Mach Number')
grid on
title('Steady State Distribution of variables inside the nozzle - Conservative Form')
subplot(5, 1, 2)
plot(x, P, LineWidth=2, Color='g')
ylabel('Pressure Ratio')
grid on
subplot(5, 1, 3)
plot(x, T, LineWidth=2, Color='b')
ylabel('Temperature Ratio')
grid on
subplot(5, 1, 4)
plot(x, rho, LineWidth=2, Color='y')
ylabel('Density Ratio')
grid on
subplot(5, 1, 5)
plot(x, mdot, LineWidth=2, Color='m')
xlabel('Nozzle Length')
ylabel('Mass Flow Rate')
grid on
figure(3)
plot(iteration_no, M_thrt, iteration_no, rho_thrt, iteration_no, T_thrt, iteration_no, P_thrt, iteration_no, mdot_thrt, LineWidth=2)
grid on
xlabel('Number of Timesteps')
legend('Mach Number', 'Density Ratio', 'Temperature Ratio', 'Pressure Ratio', 'Mass flow rate')
title('Time-wise variation of variable at nozzle throat - Conservative Form')
figure(4)
plot(iteration_no, residual1, iteration_no, residual2, iteration_no, residual3, LineWidth=2)
grid on
xlabel('Number of Timesteps')
ylabel('Residuals')
legend('U1 Residuals', 'U2 Residuals', 'U3 Residuals')
title('Residual Values - Conservative Form')
else
disp('Invalid Input')
return
end
%Non-Conservative Form
function [x, v, T, rho, M, P, mdot, M_thrt, rho_thrt, T_thrt, P_thrt, mdot_thrt, iteration_no, residual1, residual2, residual3] = NCF(n, x, gamma)
dx = x(2)-x(1);
%Initial Conditions
rho_i = 1 - (0.3146*x);
T_i = 1 - (0.2314*x);
v_i = (0.1+(1.09*x)).*(T_i.^0.5);
A = 1 + (2.2*((x-1.5).^2));
%Grid point for throat
nzzl_thrt = find(A==1);
%Field variables
rho = rho_i;
T = T_i;
v = v_i;
%Time Domain (iterations carried out=1400)
C = 0.5;
for i=1:n
dt_i(i)=C*(dx/(((T_i(i))^0.5)+(v_i(i))));
end
dt = min(dt_i);
nt = 1400; %Number of Iterations
%Macormack method
for j=1:nt
rho_old = rho;
v_old = v;
T_old = T;
for k = 2:n-1
%Defining spatial derivatives according to forward difference scheme
drho_dx = (rho(k+1) - rho(k))/dx;
dlnA_dx = (log(A(k+1)) - log(A(k)))/dx;
dv_dx = (v(k+1) - v(k))/dx;
dT_dx = (T(k+1) - T(k))/dx;
% Predictor step
drho_dt_p(k) = -(rho(k)*dv_dx) - ((rho(k)*v(k))*dlnA_dx) - (v(k)*drho_dx); %Continuity Equation
dv_dt_p(k) = -(v(k)*dv_dx) - ((1/gamma)*(dT_dx + ((T(k)/rho(k))*drho_dx))); %Momentum Equation
dT_dt_p(k) = -(v(k)*dT_dx) - ((gamma-1)*T(k)*(dv_dx + (v(k)*dlnA_dx))); %Energy Equation
% Solution update
rho(k) = rho_old(k)+(drho_dt_p(k)*dt);
v(k) = v_old(k)+(dv_dt_p(k)*dt);
T(k) = T_old(k)+(dT_dt_p(k)*dt);
end
for l=2:n-1
% Defining spatial derivatives according to backward difference scheme
drho_dx = (rho(l) - rho (l-1))/dx;
dlnA_dx = (log(A(l)) - log(A(l-1)))/dx;
dv_dx = (v(l) - v(l-1))/dx;
dT_dx = (T(l) - T(l-1))/dx;
% Corrector step
drho_dt_c(l) = -(rho(l) *dv_dx) - ((rho(l)*v(l))*dlnA_dx) - (v(l) *drho_dx);
dv_dt_c(l) = -(v(l) *dv_dx) - ((1/gamma)*(dT_dx + ((T(l)/rho(l))*drho_dx)));
dT_dt_c(l) = -(v(l)*dT_dx) - ((gamma-1)*T(l)*(dv_dx + (v(l)*dlnA_dx)));
end
for m = 2:n-1
%Average value for the particular time step
drho_dt(m) = (0.5*dt*(drho_dt_p(m) + drho_dt_c(m)));
dv_dt(m) = 0.5*dt*(dv_dt_p(m) + dv_dt_c(m));
dT_dt(m) = (0.5*dt*(dT_dt_p(m) + dT_dt_c(m)));
%Final Solution Update
rho(m) = rho_old(m) + drho_dt(m);
v(m) = v_old(m) + dv_dt(m);
T(m) = T_old(m) + dT_dt(m);
end
% Boundary Conditions
% Inlet
v(1) = (2*v(2)) - v(3);
% Outlet
rho(n) = (2*rho(n-1)) - rho(n-2);
v(n) = (2*v(n-1)) - v(n-2);
T(n) = (2*T(n-1)) - T(n-2);
mdot = rho.*v.*A; %Mass flow rate
M = v./(T.^0.5); %Mach number
P = rho.*T; %Pressure ratio
% Field variables at throat section
rho_thrt(j)=rho(nzzl_thrt);
T_thrt(j) = T(nzzl_thrt);
P_thrt(j) = P(nzzl_thrt);
M_thrt(j) = M(nzzl_thrt);
mdot_thrt(j) = rho(nzzl_thrt)*v(nzzl_thrt)*A(nzzl_thrt);
iteration_no(j) = j;
%Residuals at throat
residual1(j) = drho_dt(nzzl_thrt);
residual2(j) = dv_dt(nzzl_thrt);
residual3(j) = dT_dt(nzzl_thrt);
%Plotting
figure(1)
if j == 10
plot(x, mdot, Color='r', LineWidth=2)
hold on
grid on
elseif j == 100
plot(x, mdot, Color='g', LineWidth=2)
elseif j == 250
plot(x, mdot, Color='b', LineWidth=2)
elseif j == 500
plot(x, mdot, Color='y', LineWidth=2)
elseif j == 1000
plot(x, mdot, Color='m', LineWidth=2)
elseif j == 1400
plot(x, mdot, Color='c', LineWidth=2)
end
end
xlabel('Nozzle Length')
ylabel('Mass Flow Rate')
legend('dt=10', 'dt=100', 'dt=250', 'dt=500', 'dt=1000', 'dt=1400')
title('Mass Flow Rate at nozzle throat calculated at different timesteps - Non-Conservative Form')
end
% Conservative Form
function [x, v, T, rho, M, P, mdot, M_thrt, rho_thrt, T_thrt, P_thrt, mdot_thrt, iteration_no, residual1, residual2, residual3] = CF(n, x, gamma)
dx = x(2)-x(1);
%Initial Conditions for Conservative Form
for h = 1:n
if (x(h) <= 0.5)
rho(h) = 1.0;
T(h) = 1.0;
elseif (x(h) > 0.5) && (x(h) <= 1.5)
rho(h) = 1.0 - 0.366*(x(h) - 0.5);
T(h) = 1.0 - 0.167*(x(h) - 0.5);
elseif (x(h) > 1.5) && (x(h) <= 3.0)
rho(h) = 0.634 - 0.3879*(x(h) - 1.5);
T(h) = 0.833 - 0.3507*(x(h) - 1.5);
end
end
%Grid point for throat
A = 1 + (2.2*((x-1.5).^2));
nzzl_thrt = find(A==1);
v = 0.59./(rho.*A);
nt = 1400; %Number of Iterations
%Source Terms
U1 = rho.*A;
U2 = U1.*v;
U3 = U1.*((T/(gamma - 1)) + (gamma/2)*v.^2);
%Time Domain (iterations carried out=1400)
C = 0.5;
for i = 1:n
dt_i(i)=C*(dx/(((T(i)^0.5) + (v(i)))));
end
dt = min(dt_i);
%Macormack method
for j=1:nt
U1_old = U1;
U2_old = U2;
U3_old = U3;
%Flux Terms
F1 = U2;
F2 = (U2.^2)./U1 + ((gamma - 1)/gamma)*(U3 - (gamma/2)*((U2.^2)./U1));
F3 = gamma*((U3.*U2)./U1) - (gamma*(gamma - 1)/2)*((U2.^3)./(U1.^2));
for k = 2:n-1
%Predictor step
dU1_dt_p(k) = -(F1(k + 1) - F1(k))/dx; %Continuity Equation
J2 = (1/gamma)*(rho(k)*T(k)*(A(k + 1) - A(k))/dx);
dU2_dt_p(k) = -((F2(k + 1) - F2(k))/dx) + J2; %Momentum Equation
dU3_dt_p(k) = -(F3(k + 1) - F3(k))/dx; %Energy Equation
%Solution update
U1(k) = U1_old(k)+(dU1_dt_p(k)*dt);
U2(k) = U2_old(k)+(dU2_dt_p(k)*dt);
U3(k) = U3_old(k)+(dU3_dt_p(k)*dt);
end
rho = U1./A;
T = (gamma - 1)*((U3./U1) - (gamma/2)*((U2./U1).^2));
v = U2./U1;
F1 = U2;
F2 = (U2.^2)./U1 + ((gamma - 1)/gamma)*(U3 - (gamma/2)*((U2.^2)./U1));
F3 = gamma*((U3.*U2)./U1) - (gamma*(gamma - 1)/2)*((U2.^3)./(U1.^2));
for l=2:n-1
% Corrector step
dU1_dt_c(l) = -(F1(l) - F1(l - 1))/dx; %Continuity Equation
J2 = (1/gamma)*(rho(l)*T(l)*(A(l) - A(l - 1))/dx);
dU2_dt_c(l) = -((F2(l) - F2(l - 1))/dx) + J2; %Momentum Equation
dU3_dt_c(l) = -(F3(l) - F3(l - 1))/dx; %Energy Equation
end
for m = 2:n-1
%Average value for the particular time step
dU1_dt(m) = (0.5*dt*(dU1_dt_p(m) + dU1_dt_c(m)));
dU2_dt(m) = 0.5*dt*(dU2_dt_p(m) + dU2_dt_c(m));
dU3_dt(m) = (0.5*dt*(dU3_dt_p(m) + dU3_dt_c(m)));
%Final Solution Update
U1(m) = U1_old(m) + dU1_dt(m);
U2(m) = U2_old(m) + dU2_dt(m);
U3(m) = U3_old(m) + dU3_dt(m);
end
%Boundary Conditions
%Inlet
U1(1) = A(1);
U2(1) = 2*U2(2) - U2(3);
U3(1) = U1(1)*((T(1))/(gamma - 1) + (gamma/2)*(v(1)^2));
%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);
%Primitive Variables Update
rho = U1./A;
T = (gamma - 1)*((U3./U1) - (gamma/2)*((U2./U1).^2));
v = U2./U1;
P = rho.*T; %Pressure ratio
M = v./(T.^0.5); %Mach number
mdot = rho.*A.*v; %Mass flow rate
%Field variables at throat section
rho_thrt(j)=rho(nzzl_thrt);
T_thrt(j) = T(nzzl_thrt);
P_thrt(j) = P(nzzl_thrt);
M_thrt(j) = M(nzzl_thrt);
mdot_thrt(j) = rho(nzzl_thrt)*v(nzzl_thrt)*A(nzzl_thrt);
iteration_no(j) = j;
%Residuals at throat
residual1(j) = dU1_dt(nzzl_thrt);
residual2(j) = dU3_dt(nzzl_thrt);
residual3(j) = dU2_dt(nzzl_thrt);
%Plotting
figure(1)
if j == 10
plot(x, mdot, Color='r', LineWidth=2)
hold on
grid on
elseif j == 100
plot(x, mdot, Color='g', LineWidth=2)
elseif j == 250
plot(x, mdot, Color='b', LineWidth=2)
elseif j == 500
plot(x, mdot, Color='y', LineWidth=2)
elseif j == 1000
plot(x, mdot, Color='m', LineWidth=2)
elseif j == 1400
plot(x, mdot, Color='c', LineWidth=2)
end
end
xlabel('Nozzle Length')
ylabel('Mass Flow Rate')
legend('dt=10', 'dt=100', 'dt=250', 'dt=500', 'dt=1000', 'dt=1400')
title('Mass Flow Rate at nozzle throat calculated at different timesteps - Conservative Form')
endGRAPHS:
1) Steady State distribution of the primitive variables inside the nozzle
2) Time-wise variation of variable at nozzle throat
3) Mass Flow Rate at nozzle throat calculated at different timesteps
Residuals observed in the MacCormack Method
Grid Independence Test
number of nodes = 31
number of nodes = 61
CONCLUSION:
The isentropic flow simulation through a quasi 1D subsonic-supersonic nozzle has been executed as required. The required graphs are also plotted to show - the Steady State variation of the variables across the nozzle length, the time wise variation of these variables calculated at the nozzle throat against the number of iterations run in the simulation, the time wise variations of the mass flow rate calculated at the nozzle throat against the number iterations run in the simulation and the residual values observed throughout the simulation for both Non-Conservative and Conservative Forms respectively.
When the two forms of governing equations used here (Non-Conservative and Conservative) are compared we can conclude that they both result in a converging solution. However, while observing the Steady State Distribution of mass flow rate in the nozzle, the time wise variation of the variables within the nozzle throat and residual values, we can observe that in the Conservative form the variations in the initial few iterations result in a much smaller error than when compared to the ones observed in Non-Conservative form and hence also have smaller variations in the residual values with respect to the inital iterations. In CFD simulations the form of governing equations used has important implications on the solutions. if there are shocks, chemical reactions or any other sharp interfaces, then the Conservative form is used. However, there are various conditions, in many engineering situations Non-Conservative form is actually preferred while solving problems with shocks. For example, the Murman-Cole scheme, which uses a Non-Conservative form of governing equations, is highly preferred for transonic potential equations. It consists of a switch between central and upwind scheme.
Also a Grid Independence test has been performed on the simulations executed here. Three grids are considered having number of nodes = 31, 61 and 91. By comparing the graphical results in each of these cases, we can conclude that smaller the gridsize (i.e. greater the number of nodes) much lesser are the variations in the variables considered in the initial iterations i.e., the time wise variations are smaller when the gridsize is smaller. Also the residuals obtained in each form show smaller variations when smaller gridsizes are used. Hence, we can also say that the solution would converge faster if smaller gridsizes are used, however also to be noted is that the fact that this would slightly increase the computational requirement.