Week 7 - Simulation of a 1D Supersonic nozzle flow simulation using Macormack Method
AIM
To numerically simulate the isentropic flow through a quasi-1D subsonic-supersonic nozzle using the MacCormack technique.
1. Deriving both the Conservation and Non-Conservation forms of the governing equations.
2. Solving both the Conservation and Non-Conservation forms using the MacCormack technique.
3. Compare the numerical analysis in MATLAB of both forms.
THEORY
Ramjets, Scramjets, and Rockets all use nozzles to accelerate hot exhaust to produce thrust, as described by Newton's third law of motion. The amount of thrust produced by the engine depends on the mass flow rate through the engine, the exit velocity of the flow, and the pressure at the exit of the engine. The values of these three flow variables are all determined by the nozzle design.
The nozzle is a device with a design that allows the hot gases to flow at high velocity to generate the much-needed thrust for motion. The configuration of the nozzle is such that the inflow of gases occurs through the converging part of the nozzle, otherwise called the Convergent section, and the outflow of gases occurs through the diverging part of the nozzle, otherwise called the Divergent section.
This nozzle configuration is called a convergent-divergent, or CD, nozzle. In a CD nozzle, the hot exhaust leaves the combustion chamber and converges down to the minimum area, or throat, of the nozzle. The throat size is chosen to choke the flow and set the mass flow rate through the system. The flow in the throat is sonic which means the Mach number is equal to one in the throat. Downstream of the throat, the geometry diverges and the flow is isentropically expanded to a supersonic Mach number that depends on the area ratio of the exit to the throat. The expansion of a supersonic flow causes the static pressure and temperature to decrease from the throat to the exit, so the amount of the expansion also determines the exit pressure and temperature. The exit temperature determines the exit speed of sound, which determines the exit velocity. The exit velocity, pressure, and mass flow through the nozzle determine the amount of thrust produced by the nozzle.
In the following problem, we consider the flow to be steady, inviscid and isentropic, flowing through a C-D nozzle as shown in the figure. The flow at the inlet comes through the reservoir. The cross-sectional area is considered to be infinity at the opening, hence the velocity is 0.
The primitive variables at the reservoir are marked as
`P_o` = Total Pressure
`T_o` = Total Temperature
`\rho_o` = Density
`V` = Velocity = 0
`A = infty`
The flow expands isentropically to supersonic speeds at the exit of the nozzle, where the exit pressure, temperature, velocity, and Mach number is `p_e`, `T_e`, `V_e`, and `M_e` respectively. The flow is locally subsonic at the convergent section, sonic at the throat (minimum area), and supersonic at the divergent section. The area of the stream tube changes as a function of x, that is, `A =A(x)`. Strictly speaking, this flow is three-dimensional; the flow-field variables are functions of x, y, and z. In particular, the velocity at the boundary of the stream tube must be tangent to the boundary, and hence it has components in the y and z directions as well as the axial x direction. However, if the area variation is moderate, the components in the y and z directions are small in comparison with the component in the x direction. In such a case, the flow-field variables can be assumed to vary with x only (i.e., the flow can be assumed to be uniform across any cross-section at a given x station). Such a flow, where `A = A(x)`, `p = p(x)`, `ρ = ρ(x)`, `V = V (x)` is defined as quasi-one-dimensional flow.
The governing equations for quasi-one-dimensional, steady, isentropic flow can be expressed as
Continuity- `\rho_1 V_1A_1 = \rho_2 V_2 A_2`
Momentum- `p_1 A_1 + \rho_1 V_1^2 A_1 + int_(A_1)^(A_2) p dA = p_2 A_2 + \rho_2 V_2^2 A_2`
Energy- `h_1 + V_1^2/2 = h_2 + V_2^2/2`
For calorically perfect gas we have
`p = \rhoRT`
`h = C_pT`
Before moving on to a solution to the governing equations, let us examine some physical characteristics of a quasi-one-dimensional flow. To help this examination, we first obtain some differential expressions for the governing equations, in contrast to the algebraic equations obtained above. The differential form of the Euler's equation are as follows-
`d(\rhoVA) = 0`
`dp = -\rhoVdV`
`dh + VdV =0`
On combining these three equations above we have
`(d\rho)/ \rho = (dV)/V = (dA)/A`
From the momentum differential equation we have
`dp = -\rhoVdV`
`(dp)/(d\rho) (d\rho)/\rho = -VdV`
From the equation of speed of sound, we have
`(dp)/(d\rho) = a^2
Substituting we get `(d\rho)/(\rho) =- V/a^2 dV`
`(d\rho)/(\rho) = -V^2/a^2 (dV)/V = -M^2 (dV)/V`
Hence comparing with the earlier relation of `\rho`, `V`, and `A`, we get
`(dA)/A = (M^2-1) (dV)/V`
1. For 0 ≤ M < 1>
2. For M > 1 (supersonic flow), the quantity in parentheses in the equation above is positive. Hence, an increase in velocity (positive dV) is associated with an increase in area (positive dA). Likewise, a decrease in velocity (negative dV) is associated with a decrease in the area (negative dA). For a supersonic flow, to increase the velocity, we must have a divergent duct, and to decrease the velocity, we must have a convergent duct. They are the direct opposite of the trends for subsonic flow.
3. For M = 1 (sonic flow), the equation shows that dA = 0 even though a finite dV exists. Mathematically, this corresponds to a local maximum or minimum in the area distribution. Physically, it corresponds to a minimum area.
Following are the Governing Equations for the Quasi 1-D Invisid, Unsteady, Isentropic Subsonic-Supersonic Nozzle Flow -
Non-dimensional Governing Equations
1. Non-Conservation form of Governing Equations
Continuity Equation - `(∂ ρ ' )/(∂ t ' )= − ρ ' (∂ V ')/( ∂ x ') − ρ ' V ' (∂ ( ln A ' ))/( ∂ x ') − V ' (∂ ρ ' )/(∂ x ')`
Momentum Equation - `(∂ V ')/(∂ t ') = − V ' (∂ V ')/( ∂ x ') − 1/ γ ( (∂ T ')/( ∂ x ') + (T ' )/(ρ ') ⋅ (∂ ρ ')/( ∂ x ') )`
Energy Equation - `(∂ T ')/( ∂ t ') = − V ' (∂ T ')/( ∂ x ') − ( γ − 1 ) . T ' ( (∂ V ')/( ∂ x ') + V ' (∂ ( ln A ' ))/( ∂ x ') )`
2. Conservation form of Governing Equations
Continuity Equation - `(∂ ( ρ ' A ' ))/ (∂ t ') + (∂ ( ρ ' A ' V ' ))/(∂ x ') = 0`
Momentum Equation - `(∂ ( ρ ' A ' V ' ))/( ∂ t ') + (∂ ( ρ ' A ' V '^ 2 + ( 1 /γ ) p ' A ' ))/ (∂ x ') = 1/ γ ⋅ p ' (∂ A ')/(∂ x ')`
Energy Equation - `(∂ ( ρ ' A ' ( (e ')/(γ − 1) + γ /2 V '^ 2 ) ) )/(∂ t ') + (∂ ( ρ ' A ' V ' ( (e ')/( γ − 1) + γ/ 2 V '^ 2 ) + p ' A ' V ' ) )/(∂ x ') = 0`
Upon simplification we get
Continuity Equation - `(∂ U_1) /(∂ t ') = − (∂ F_1)/(∂ x ')`
Momentum Equation - `(∂ U_2) /(∂ t ') = − (∂ F_2 )/(∂ x ') + J_2`
Energy Equation - `(∂ U_3) /(∂ t ') = − (∂ F_3 )/(∂ x ')`
Discretization
To implement a finite-difference solution, we divide the Nozzle length into a number of discrete grid points along the X - direction, as shown in the following. In this problem, we make the assumption that flow variables only vary along the X - direction. The variation along the Y - direction is neglected as it is very small.
MacCormack Method
MacCormack method is a second-order explicit finite difference technique, a widely used discretization scheme for the numerical solution of hyperbolic partial differential equations.
Let f(x) be a function. Then the following steps are followed
i. Initial profile
Before the calculation, the initial profile for all the flow field variables which are velocity, density and temperature has to be determined. These are estimated values in the boundary.
ii. Predictor step
`((df)/dt)_p = (f(i+1) -f(i))/(trianglex)`
iii. Solution update
`f_p = f_ i + ((df)/dt)_p. trianglet`
`f_p`= Flow field variables calculated after the predictor step
`f_i` = Flow field variable in the initial step
`trianglet` = time step size
iv. Corrector step
`((df_p)/dt)_c = (f(i) - f(i-1))/ (trianglex)`(∂fdt)avg=0.5((∂fdt)p+(∂fpdt)c)">
`((df)/dt)_(avg) = 0.5(((df)/dt)_p + ((df_p)/dt)_c)`
v. Final solution update
`f_(sol) = f_ i + ((df)/dt)_(avg). trianglet`
`f_(sol)` = Final solution of the flow field variables for the first time step
Boundary Conditions
Initial Conditions
Nozzle Shape
The nozzle shape, `A = A ( x )`, is specified and held fixed, independent of time. For this case the parabolic area distribution given by
`A = 1 + 2.2 ( x − 1.5 )^ 2` `0 ≤ x ≤ 3`
Note that `x = 1.5` is the throat of the nozzle, that the convergent section occurs for `x < 1> and that the divergent section occurs for `x > 1.5` as shown in the figure.
Two initial conditions are considered based on the Non-Conservative or Conservative approach taken
1. Non-Conservative approach
2. Conservative approach
MATLAB Program
1. Main Program
clear all
close all
clc
%SOLVING QUASI 1D SUPERSONIC NOZZLE FLOW EQUATIONS
%INPUT PARAMETERS
%No. of grid points
n = input('Enter the No. of Nodes or Grid Size [Recommended 31 or 61]) : ');
gamma = 1.4;
x = linspace(0,3,n);
dx = x(2) - x(1);
%Area & Throat profiles
a = 1 + 2.2*(x - 1.5).^2;
b = a.*(-1);
%No. of time step
nt = 1400;
%Courant Number
C = 0.5;
ct = 1;
%FUNCTION FOR ANALYTICAL FORM
tic;
[Mass_flow_rate_A, Pressure_A, M_A, rho_A, T_A, i, a, b] = Analytical_form(x, dx, n, gamma, nt, C)
Time_elapsed_A = toc;
sprintf('Time Elasped for Analytical form: %0.4f', Time_elapsed_A)
% Plotting the Exact Analytical Solution
figure(7)
Supt_A = suptitle('Numerical Solution by Analytical Approach');
set(Supt_A,'FontSize',10,'FontWeight','normal');
subplot(5,1,1)
plot(x,a,'color',[0.8 0.1 0.1],'linewidth',1.5);
title('Supersonic Nozzle Profile along X-Direction')
hold on;
plot(x,b,'color',[0.9 0.1 0.1],'linewidth',1.5);
plot([x(round(n/2)) x(round(n/2))],[-5 5],':','color','k');
plot(x,0,'--ko','MarkerSize',2,'MarkerEdgeColor','k','MarkerFaceColor',[0.7,0.7,0.7]);
plot([x(1) x(n)],[0 0],':','color','k')
subplot(5,1,2)
plot(x,M_A,'-o','color','r','linewidth',1.5,'MarkerSize',2,'MarkerEdgeColor','k','MarkerFaceColor','k');
title('Mach Number vs x / L');
grid on;
xlabel('x / L');
ylabel('M');
subplot(5,1,3)
plot(x,Pressure_A,'-o','color','r','linewidth',1.5,'MarkerSize',2,'MarkerEdgeColor','k','MarkerFaceColor','k');
title('Pressure Ratio vs x / L');
grid on;
xlabel('x / L');
ylabel('P/Po');
subplot(5,1,4)
plot(x,rho_A,'-o','color','r','linewidth',1.5,'MarkerSize',2,'MarkerEdgeColor','k','MarkerFaceColor','k');
title('Density Ratio vs x / L');
grid on;
xlabel('x / L');
ylabel('\rho/\rhoo');
subplot(5,1,5)
plot(x,T_A,'-o','color','r','linewidth',1.5,'MarkerSize',2,'MarkerEdgeColor','k','MarkerFaceColor','k');
title('Temperature Ratio vs x / L');
grid on;
xlabel('x / L');
ylabel('T/To');
%FUNCTION FOR NON CONSERVATIVE FORM
tic;
[Mass_flow_rate_NC, Numerical_Mass_flow_rate_NC, Pressure_NC, Mach_number_NC, rho_NC, v_NC, T_NC, rho_throat_NC, v_throat_NC, T_throat_NC, Mass_flow_rate_throat_NC, Pressure_throat_NC, Mach_number_throat_NC, k] = Non_conservative_form(x, dx, n, gamma, nt, C);
Time_elapsed_NC = toc;
sprintf('Time Elasped for Non Conservative form: %0.4f', Time_elapsed_NC);
%PLOTTING
%Comparison of non-dimensional Mass flow rates at different time steps
figure(2)
for q = 1:6
plot(x,Numerical_Mass_flow_rate_NC(q,:),'linewidth',1.5);
hold on;
end
%Labelling
title('Mass flow rates at different Timesteps - Non-Conservative form')
xlabel('Distance')
ylabel('Mass flow rate')
plot(x,Mass_flow_rate_A,'-o','color',[0.85 0.3 0.1],'linewidth',1.5,'MarkerSize',2,'MarkerEdgeColor','k','MarkerFaceColor','k');
legend({'Timestep = 50','Timestep = 100','Timestep = 200','Timestep = 400','Timestep = 700', 'Timestep = 1400', 'Analytical Value'})
axis([0 3 0 2])
grid on
hold on
%Steady state simulation for primitve values for Non Conservative form
figure(3);
Supt_NC = suptitle('Numerical Solution by Non-Conservation Approach');
set(Supt_NC,'FontSize',10,'FontWeight','normal');
subplot(5,1,1)
plot(x,a,'color',[0.8 0.1 0.1],'linewidth',1.5);
hold on;
plot(x,b,'color',[0.9 0.1 0.1],'linewidth',1.5);
title('Supersonic Nozzle Profile along X-Direction')
hold on;
plot([x(round(n/2)) x(round(n/2))],[-5 5],':','color','k');
plot(x,0,'--ko','MarkerSize',2,'MarkerEdgeColor','k','MarkerFaceColor',[0.7,0.7,0.7]);
plot([x(1) x(n)],[0 0],':','color','k')
�nsity
subplot(5,1,2);
plot(x,rho_NC,'color','b','linewidth',1)
ylabel('\rho/\rhoo')
axis([0 3 0 1])
grid on;
title('Variation of Flow Rate Distribution - Non-Conservative')
%Temperature
subplot(5,1,3);
plot(x,T_NC,'color','b','linewidth',1)
ylabel('T/To')
axis([0 3 0 1])
grid on;
%Pressure
subplot(5,1,4);
plot(x,Pressure_NC,'color','b','linewidth',1)
ylabel('P/Po')
axis([0 3 0 1])
grid on;
%Mach number
subplot(5,1,5);
plot(x,Mach_number_NC,'color','b','linewidth',1)
xlabel('Nozzle Length')
ylabel('Mach number')
axis([0 3 0 4])
grid on;
%FUNCTION FOR CONSERVATIVE FORM
tic;
[Mass_flow_rate_C, Numerical_Mass_flow_rate_C, Pressure_C, Mach_number_C, rho_C, v_C, T_C, rho_throat_C, v_throat_C, T_throat_C, Mass_flow_rate_throat_C, Pressure_throat_C, Mach_number_throat_C, i] = Conservative_form(x, dx, n, gamma, nt, C);
Time_elapsed_C = toc;
sprintf('Time Elasped for Conservative form: %0.4f', Time_elapsed_C);
%PLOTTING
%Comparison of non-dimensional Mass flow rates at different time steps
figure(5)
for q = 1:6
plot(x,Numerical_Mass_flow_rate_C(q,:),'linewidth',1.5);
hold on;
end
%Labelling
title('Mass flow rates at different Timesteps - Conservative Form)')
xlabel('Distance')
ylabel('Mass flow rate')
plot(x,Mass_flow_rate_A,'-o','color',[0.85 0.3 0.1],'linewidth',1.5,'MarkerSize',2,'MarkerEdgeColor',[0.5 0.2 0.2],'MarkerFaceColor',[0.85 0.3 0.1]);
axis([0 3 0 2])
grid on
hold on
legend({'Timestep = 50','Timestep = 100','Timestep = 200','Timestep = 400','Timestep = 700', 'Timestep = 1400', 'Analytical Value'})
%Steady state simulation for primitve values for Conservative form
figure(6);
Supt_C = suptitle('Numerical Solution by Conservative Approach');
set(Supt_C,'FontSize',10,'FontWeight','normal');
subplot(5,1,1)
plot(x,a,'color',[0.8 0.1 0.1],'linewidth',1.5);
title('Supersonic Nozzle Profile along X-Direction')
hold on;
plot(x,b,'color',[0.9 0.1 0.1],'linewidth',1.5);
plot([x(round(n/2)) x(round(n/2))],[-5 5],':','color','k');
plot(x,0,'--ko','MarkerSize',2,'MarkerEdgeColor','k','MarkerFaceColor',[0.7,0.7,0.7]);
plot([x(1) x(n)],[0 0],':','color','k')
�nsity
subplot(5,1,2);
plot(x,rho_C,'color','b','linewidth',1)
ylabel('\rho/\rhoo')
axis([0 3 0 1])
grid on;
title('Variation of Flow Rate Distribution - Conservative')
%Temperature
subplot(5,1,3);
plot(x,T_C,'color','b','linewidth',1)
ylabel('T/To');
axis([0 3 0 1]);
grid on;
%Pressure
subplot(5,1,4);
plot(x,Pressure_C,'color','b','linewidth',1)
ylabel('P/Po')
axis([0 3 0 1])
grid on;
%Mach number
subplot(5,1,5);
plot(x,Mach_number_C,'color','b','linewidth',1)
xlabel('Nozzle Length')
ylabel('Mach number')
axis([0 3 0 4])
grid on;
%COMPARISON PLOT of Normalized Mass Flow Rate of both forms
figure(8)
hold on
plot(x,Mass_flow_rate_NC,'r','linewidth',1.25)
hold on;
plot(x,Mass_flow_rate_C,'b','linewidth',1.25)
hold on;
grid on
legend('Non Conservative','Conservative');
xlabel('Length of the Domain')
ylabel('Mass Flow Rate')
title('Comparison of Normalized Mass Flow Rate of both forms')2. Function for Analytical Form
% Analytical Calculation
function [Mass_flow_rate_A, Pressure_A, M_A, rho_A, T_A, i, a, b] = Analytical_form(x, dx, n, gamma, nt, C)
a = 1+2.2*(x-1.5).^2; % Area
b = a.*(-1);
M = 0;
% Mach Number
for i = 1:round(n/2) % From Inlet to Throat
resi = 1;
while resi > 1e-9
M = M + 1e-4;
resi = ((1/M^2)*((2/(gamma+1))*(1+((gamma-1)/2)*M^2))^((gamma+1)/(gamma-1)))-(a(i)^2);
end
M_A(i) = M;
end
for i = round(n/2):n % From Throat to Outlet
resi = 1;
while resi > 1e-9
M = M + 1e-4;
resi = ((1/M^2)*((2/(gamma+1))*(1+((gamma-1)/2)*M^2))^((gamma+1)/(gamma-1)))-(a(i)^2);
resi = resi*(-1);
end
M_A(i) = M;
end
rho_A = (1+(((gamma-1)/2).*M_A.^2)).^(-1/(gamma-1));
T_A = (1+(((gamma-1)/2).*(M_A.^2))).^(-1);
Pressure_A = rho_A.*T_A;
Mass_flow_rate_A = ones(n).*0.579;3. Function for Non-Conservative Form
%Function for Non-Conservative form
function [Mass_flow_rate_NC, Numerical_Mass_flow_rate_NC, Pressure_NC, Mach_number_NC, rho_NC, v_NC, T_NC, rho_throat_NC, v_throat_NC, T_throat_NC, Mass_flow_rate_throat_NC, Pressure_throat_NC, Mach_number_throat_NC, dt, k] = Non_conservative_form(x, dx, n, gamma, nt, C)
%Intial boundary conditions
rho_NC = 1 - 0.3146*x;
T_NC = 1 - 0.2341*x;
v_NC = (0.1 + 1.09*x).*T_NC.^0.5;
%Area & Throat profiles
a = 1 + 2.2*(x - 1.5).^2;
b = a.*(-1);
throat = find(a==1);
%Outer Time Loop
for k = 1:nt
%Time-step control
dt = min(C.*dx./(sqrt(T_NC) + v_NC));
%Saving a copy of Old Values
rho_old = rho_NC;
T_old = T_NC;
v_old = v_NC;
%Predictor Method
for j = 2:n-1
%Simplifying the Equation's spatial terms
dv_dx = (v_NC(j + 1) - v_NC(j))/dx;
drho_dx = (rho_NC(j + 1) - rho_NC(j))/dx;
dT_dx = (T_NC(j + 1) - T_NC(j))/dx;
dlog_a_dx = (log(a(j + 1)) - log(a(j)))/dx;
%Continuity Equation
drho_dt_P(j) = -rho_NC(j)*dv_dx - rho_NC(j)*v_NC(j)*dlog_a_dx - v_NC(j)*drho_dx;
%Momentum Equation
dv_dt_P(j) = -v_NC(j)*dv_dx - (1/gamma)*(dT_dx + (T_NC(j)*drho_dx/rho_NC(j)));
%Energy Equation
dT_dt_P(j) = -v_NC(j)*dT_dx - (gamma - 1)*T_NC(j)*(dv_dx + v_NC(j)*dlog_a_dx);
%Solution Update
rho_NC(j) = rho_NC(j) + drho_dt_P(j)*dt;
v_NC(j) = v_NC(j) + dv_dt_P(j)*dt;
T_NC(j) = T_NC(j) + dT_dt_P(j)*dt;
end
%Corrector Method
for j = 2:n-1
%Simplifying the Equation's spatial terms
dv_dx = (v_NC(j) - v_NC(j - 1))/dx;
drho_dx = (rho_NC(j) - rho_NC(j - 1))/dx;
dT_dx = (T_NC(j) - T_NC(j - 1))/dx;
dlog_a_dx = (log(a(j)) - log(a(j - 1)))/dx;
%Continuity Equation
drho_dt_C(j) = - rho_NC(j)*dv_dx - rho_NC(j)*v_NC(j)*dlog_a_dx - v_NC(j)*drho_dx;
%Momentum Equation
dv_dt_C(j) = - v_NC(j)*dv_dx - (1/gamma)*(dT_dx + (T_NC(j)*drho_dx/rho_NC(j)));
%Energy Equation
dT_dt_C(j) = - v_NC(j)*dT_dx - (gamma - 1)*T_NC(j)*(dv_dx + v_NC(j)*dlog_a_dx);
end
%Compute Average Time derivative
drho_dt = 0.5.*(drho_dt_P + drho_dt_C);
dv_dt = 0.5.*(dv_dt_P + dv_dt_C);
dT_dt = 0.5.*(dT_dt_P + dT_dt_C);
%Final Solution Update
for i = 2:n-1
rho_NC(i) = rho_old(i) + drho_dt(i)*dt;
v_NC(i) = v_old(i) + dv_dt(i)*dt;
T_NC(i) = T_old(i) + dT_dt(i)*dt;
end
%Applying Boundary Conditions
%At Inlet
rho_NC(1) = 1;
v_NC(1) = 2*v_NC(2) - v_NC(3);
T_NC(1) = 1;
%At Outlet
rho_NC(n) = 2*rho_NC(n-1) - rho_NC(n - 2);
v_NC(n) = 2*v_NC(n - 1) - v_NC(n - 2);
T_NC(n) = 2*T_NC(n - 1) - T_NC(n - 2);
�fining Mass flow rate, Pressure and Mach number
Mass_flow_rate_NC = rho_NC.*a.*v_NC;
Pressure_NC = rho_NC.*T_NC;
Mach_number_NC = (v_NC./sqrt(T_NC));
�lculating Variables at throat section
rho_throat_NC(1,k) = rho_NC(throat);
T_throat_NC(1,k) = T_NC(throat);
v_throat_NC(1,k) = v_NC(throat);
Mass_flow_rate_throat_NC(1,k) = Mass_flow_rate_NC(throat);
Pressure_throat_NC(1,k) = Pressure_NC(throat);
Mach_number_throat_NC(1,k) = Mach_number_NC(throat);
if k==50
Numerical_Mass_flow_rate_NC(1,:) = rho_NC.*v_NC.*a;
elseif k==100
Numerical_Mass_flow_rate_NC(2,:) = rho_NC.*v_NC.*a;
elseif k==200
Numerical_Mass_flow_rate_NC(3,:) = rho_NC.*v_NC.*a;
elseif k==400
Numerical_Mass_flow_rate_NC(4,:) = rho_NC.*v_NC.*a;
elseif k==700
Numerical_Mass_flow_rate_NC(5,:) = rho_NC.*v_NC.*a;
elseif k==1400
Numerical_Mass_flow_rate_NC(6,:) = rho_NC.*v_NC.*a;
end
%PLOTTING
%Timestep variation of properties at throat area in Non Conservative form
ts = 1:k;
figure(1)
if rem(k,50)==0
�nsity
subplot(4,1,1)
plot(ts,rho_throat_NC,'color','k','linewidth',1.5)
ylabel('\rho/\rhoo');
title('Timestep variation at throat area - Non-Conservative form')
grid on;
hold on;
%Temperature
subplot(4,1,2)
plot(ts,T_throat_NC,'color','k','linewidth',1.5)
ylabel('T/To');
grid on;
hold on;
%Pressure
subplot(4,1,3)
plot(ts,Pressure_throat_NC,'color','k','linewidth',1.5)
ylabel('P/Po');
grid on;
hold on;
%Mach Number
subplot(4,1,4)
plot(ts,Mach_number_throat_NC,'color','k','linewidth',1.5)
xlabel('Timesteps')
ylabel('Mach Number');
grid on;
hold on;
pause(0.01)
end
end
end
4. Function for Conservative Form
%Function for Conservative form
function [Mass_flow_rate_C, Numerical_Mass_flow_rate_C, Pressure_C, Mach_number_C, rho_C, v_C, T_C, rho_throat_C, v_throat_C, T_throat_C, Mass_flow_rate_throat_C, Pressure_throat_C, Mach_number_throat_C, i] = Conservative_form(x, dx, n, gamma, nt, C)
�termining Initial profile of Density and Temperature
for i = 1:n
if (x(i) >= 0 && x(i) <= 0.5)
rho_C(i) = 1;
T_C(i) = 1;
elseif (x(i) >= 0.5 && x(i) <= 1.5)
rho_C(i) = 1 - 0.366*(x(i) - 0.5);
T_C(i) = 1 - 0.167*(x(i) - 0.5);
elseif (x(i) >= 1.5 && x(i) <= 3)
rho_C(i) = 0.634 - 0.3879*(x(i) - 1.5);
T_C(i) = 0.833 - 0.3507*(x(i) - 1.5);
end
end
%Area & Throat profiles
a = 1 + 2.2*(x - 1.5).^2;
b = a.*(-1);
throat = find(a==1);
%Initial profile fot Velocity and Pressure
v_C = 0.59./(rho_C.*a);
Pressure_C = rho_C.*T_C;
�lculate the Initial Solution Vectors(U)
U1 = rho_C.*a;
U2 = rho_C.*a.*v_C;
U3 = rho_C.*a.*(T_C./(gamma - 1)) + ((gamma/2).*v_C.^2);
%Outer Time Loop
for i =1:nt
%Time-step control
dt = min(C.*dx./(sqrt(T_C) + v_C));
%Saving a copy of Old Values
U1_old = U1;
U2_old = U2;
U3_old = U3;
�lculating the Flux vector(F)
F1 = U2;
F2 = ((U2.^2)./U1) + ((gamma - 1)/gamma)*(U3 - (gamma/2)*((U2.^2)./U1));
F3 = ((gamma*U2.*U3)./U1) - (gamma*(gamma - 1)*((U2.^3)./(2*U1.^2)));
%Predictor Method
for j = 2:n-1
�lculating the Source term(J)
J2(j) = (1/gamma)*rho_C(j)*T_C(j)*((a(j + 1) - a(j))/dx);
%Continuity Equation
dU1_dt_P(j) = -((F1(j + 1) - F1(j))/dx);
%Momentum Equation
dU2_dt_P(j) = -((F2(j + 1) - F2(j))/dx) + J2(j);
%Energy Equation
dU3_dt_P(j) = -((F3(j + 1) - F3(j))/dx);
%Updating New values
U1(j) = U1(j) + dU1_dt_P(j)*dt;
U2(j) = U2(j) + dU2_dt_P(j)*dt;
U3(j) = U3(j) + dU3_dt_P(j)*dt;
end
%Updating Flux vectors
F1 = U2;
F2 = ((U2.^2)./U1) + ((gamma - 1)/gamma)*(U3 - ((gamma*0.5*(U2.^2))./(U1)));
F3 = ((gamma*U2.*U3)./U1) - ((gamma*0.5*(gamma - 1)*(U2.^3))./(U1.^2));
%Corrector Method
for j = 2:n-1
%Updating the Source term(J)
J2(j) = (1/gamma)*rho_C(j)*T_C(j)*((a(j) - a(j - 1))/dx);
%Simplifying the Equation's spatial terms
dU1_dt_C(j) = -((F1(j) - F1(j - 1))/dx);
dU2_dt_C(j) = -((F2(j) - F2(j - 1))/dx) + J2(j);
dU3_dt_C(j) = -((F3(j) - F3(j - 1))/dx);
end
%Compute Average Time derivative
dU1_dt = 0.5*(dU1_dt_P + dU1_dt_C);
dU2_dt = 0.5*(dU2_dt_P + dU2_dt_C);
dU3_dt = 0.5*(dU3_dt_P + dU3_dt_C);
%Final Solution Update
for k = 2:n-1
U1(k) = U1_old(k) + dU1_dt(k)*dt;
U2(k) = U2_old(k) + dU2_dt(k)*dt;
U3(k) = U3_old(k) + dU3_dt(k)*dt;
end
%Applying Boundary Conditions
%At Inlet
U1(1) = rho_C(1)*a(1);
U2(1) = 2*U2(2) - U2(3);
v_C(1) = U2(1)./U1(1);
U3(1) = U1(1)*((T_C(1)/(gamma - 1)) + ((gamma/2)*(v_C(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);
%Final Primitive value update
rho_C = U1./a;
v_C = U2./U1;
T_C = (gamma - 1)*((U3./U1 - (gamma/2)*(v_C).^2));
�fining Mass flow rate, Pressure and Mach number
Mass_flow_rate_C = rho_C.*a.*v_C;
Pressure_C = rho_C.*T_C;
Mach_number_C = (v_C./sqrt(T_C));
�lculating Variables at throat section
rho_throat_C(i) = rho_C(throat);
T_throat_C(i) = T_C(throat);
v_throat_C(i) = v_C(throat);
Mass_flow_rate_throat_C(i) = Mass_flow_rate_C(throat);
Pressure_throat_C(i) = Pressure_C(throat);
Mach_number_throat_C(i) = Mach_number_C(throat);
figure(4)
if i==50
Numerical_Mass_flow_rate_C(1,:) = rho_C.*v_C.*a;
elseif i==100
Numerical_Mass_flow_rate_C(2,:) = rho_C.*v_C.*a;
elseif i==200
Numerical_Mass_flow_rate_C(3,:) = rho_C.*v_C.*a;
elseif i==400
Numerical_Mass_flow_rate_C(4,:) = rho_C.*v_C.*a;
elseif i==700
Numerical_Mass_flow_rate_C(5,:) = rho_C.*v_C.*a;
elseif i==1400
Numerical_Mass_flow_rate_C(6,:) = rho_C.*v_C.*a;
end
%PLOTTING
%Timestep variation of properties at throat area in Conservative form
ts = 1:i;
figure(4)
if rem(i,50)==0
�nsity
subplot(4,1,1)
plot(ts,rho_throat_C,'color','k','linewidth',1.5)
ylabel('\rho/\rhoo');
title('Timestep variation at throat area - Conservative Form')
grid on;
hold on;
%Temperature
subplot(4,1,2)
plot(ts,T_throat_C,'color','k','linewidth',1.5)
ylabel('T/To');
grid on;
hold on;
%Pressure
subplot(4,1,3)
plot(ts,Pressure_throat_C,'color','k','linewidth',1.5)
ylabel('P/Po');
grid on;
hold on;
%Mach Number
subplot(4,1,4)
plot(ts,Mach_number_throat_C,'color','k','linewidth',1.5)
xlabel('Timesteps')
ylabel('Mach Number');
grid on;
hold on;
pause(0.01)
end
end
endRESULT AND DISCUSSION
1. Analytic Solution
a. Characteristics Plot of Flow-Field Variables across the Nozzle at different Nodes along the X-Direction For the Analytical Solution.
2. Non-Conservative form solution
a. Characteristics Plot of Flow Variables across Nozzle at different Nodes along the X-Direction for the Non-Conservative Method.
b. Time variation plot of non-dimension Characteristic Flow Variables at the nozzle throat for the Non-Conservative Method.
c. Mass Flow rate at Different No. of Time Steps For Non-Conservative Method.
3. Conservative form solution
a. Characteristics Plot of Flow Variables across Nozzle at different Nodes along the X-Direction for the Conservative Method.
b. Time variation plot of non-dimension Characteristic Flow Variables at the nozzle throat for the Conservative Method.
c. Mass Flow rate at Different No. of Time Steps for Conservative Method.
Comparison of Steady solution Non-Dimensional Mass Flow Rate between Non-Conservative and Conservative forms.
Although the final solution of both forms is almost similar, the Non-Conservative form oscillates more than the Conservative form. The latter also gives an almost stable solution along the length. Hence, we can conclude that the Conservative form attains stability in a lesser number of iterations with lesser time.
GRID INDEPENDENCE TEST
Following is the sum up of the results we got after the simulation. (The table compares the values of the flow variables at the Nozzle Throat)
In order for the grid independency test, we simulated twice with a different number of nodes. First with 31 Grid Points and second with 61 Grid Points. It is observed that the 61 grid points give results closer to the analytical solution and therefore are considered to be more accurate. Increasing the number of grids would give an even more accurate result, however compromising the computation time.
From the time-variation graph, it can be seen that the for Non-Conservative form, the flow variables reach a steady state in about 500-600 time steps. However, for the Conservative form, the flow variables reach the steady state only after 700 time steps. This shows that the Non-Conservative form converges faster than the Conservative form.
-----------------------------------------------------------------------------------------------------------------------------