Simulation of a 1D super-sonic nozzle flow using MacCormack Method
Objective
Simulation of a 1D super-sonic nozzle flow using the following parameters
The flow of information in this report will be as follows
- Understanding of the problem
- Main Program
- Output
- Inference
1. Understanding of the problem
In this report, we will be solving a 1D super-sonic nozzle flow with a steady and isentropic flow characteristic through a convergent-divergent nozzle shown in Fig(a).
Fig. (a)
The following consideration was taken into account before stepping into the process.
(a) The flow to the nozzle from the inlet comes from a reservoir where the pressure and temperature are denoted by `ρ_0` (Total Pressure) and `T_0` (Total Temperature), respectively.
(b) The cross-sectional area of the reservoir is considered to be very large and hence the velocity at the inlet from the reservoir is very low.
(c) The flow is subsonic in the convergent section, sonic at the nozzle throat and supersonic in the divergent section.
(d) The flow is considered to be isentropically expanded to supersonic speed at the nozzle exit where the pressure, temperature Mach number, and velocity are denoted by `p_e, T_e, M_e` and `V_e`, respectively.
Governing Equation
To solve the above-described problem, we need to formulate the governing equations. Well, good thing we have these at our disposal. Let’s have a look at them.
Continuity:
`∂/(∂t) ∭_∀ρd∀+∬_SρV⋅dS=0`
Momentum:
`∭_∀(∂ρV)/(∂t)∀+∬_SρV(V⋅dS)=∭_∀ρFd∀-∬_SpdS`
Energy:
`∂/(∂t) ∭_∀ρ(e+V^2/2)d∀+∬_S(e+V^2/2)(ρV⋅dS)+∬_SpV.dS+w ̇ _(visc)+w ̇_(body) =∭_∀q ̇∀ρd∀`
where ∀ = control volume
ρ = density
V = velocity
S = surface on the control volume
`w ̇ _(visc)` = Viscous forces
`w ̇ _(body)` = Body forces
`e` = internal energy
`q ̇` = heat added per unit mass per unit time.
Now that we have our governing equations, why wait? Well, even though the above equations were written in the general sense and that would mean just writing it down for one-dimensional domain would manufacture us with the solution, such is not the case for the problem we are trying to solve. The best way to define the reason behind such a phenomenon is in the way we twisted the scope of our problem; we assume that the flow properties will remain uniform across any chosen cross-section in a respective dimension while in reality, the nozzle flow occurs in two-dimension. Therefore, as we have made such an assumption, calling it quasi-one-dimensional flow, the above equations aren’t necessarily appropriate for us. For quasi-one-dimensional flow, we will consider `q ̇` = 0, `w ̇ _(visc)`= 0 and `w ̇ _(body)` = 0.
Now, before manipulation, we shall be using normalized primitive terms for easier computation and representation of data, the following formulation was used for normalizing the primitive terms.
`ρ=ρ/ρ_0`
`T=T/T_0`
`V=V/a_0 , a_0` = speed of sound in the reservoir `[ a_0 = √(T_0 ) ]`
`t= (t/L)/a_0, L`= domain length
`x=x/L`
`A=A/A^** , A^**` = sonic throat area.
Now that we have normalized the primitive terms, we can proceed to the manipulation of the above generalized integral form of governing equations for quasi-one-dimensional flow. Upon manipulation, we get the following formulation for conservation and non-conservation form.
Conservation Form
Continuity:
`(∂(ρA))/(∂t)=-(∂(ρAV))/(∂x)`
Momentum:
`∂/(∂t)(ρAV)=-∂/(∂x)[ρAV^2+(1/γ)pA]+( 1/γ )p (∂A)/(∂x)`
Energy:
`∂/(∂t)[ρ(e/(γ-1)+γ/2 V^2)A]=-∂/(∂x)[ρ(e/(γ-1)+γ/2 V^2 )VA+pAV]`
Upon further simplification, the above equation can be written as
Continuity:
`(∂U_1)/(∂t)=-(∂F_1)/(∂x)` (1)
Momentum:
`(∂U_2)/(∂t)=-(∂F_2)/(∂x)+J_2` (2)
Energy:
`(∂U_3)/(∂t)=-(∂F_3)/(∂x)` (3)
where `U_1= ρA`
`U_2= ρAV`
`U_3= [ρ(e/(γ-1)+γ/2 V^2)A]`
`F_1= ρAV`
`F_2= [ρAV^2+(1/γ)pA]`
`F_3= [ρ(e/(γ-1)+γ/2 V^2 )VA+pAV]`
`J_2= (1/γ )p (∂A)/(∂x)`
To make it even simpler for the computer to perform the simulation for eq. (1), (2) and (3), we shall compute the flux vectors `F_1, F_2,` and `F_3` with respect to solution vectors `U_1, U_2,` and `U_3`. Therefore,
`F_1= U_2`
`F_2= (U_2^2)/U_1 +(γ-1)/γ (U_3-γ/2 (U_2^2)/U_1 )`
`F_3= (γU_2 U_3)/U_1 +(γ(γ-1))/2 ((U_3^2)/(U_1^2 ))`
`J_2= (γ-1)/γ (U_3-γ/2 (U_2^2)/U_1 ) (∂(ln(A)))/(∂x)`
To obtain primitive values from the solution vectors, we can now break down the solutions as
`ρ=U_1/A`
`V=U_2/U_1`
`T=e= (γ-1)(U_3/U_1 -γ/2 V^2 )`
`p=ρT`
Non-Conservation Form
Continuity:
`(∂ρ)/(∂t)=-V (∂ρ)/(∂x)-ρ (∂V)/(∂x)-ρV (∂(ln(A)))/(∂x)` (4)
Momentum:
`(∂V)/(∂t)=-V (∂V)/(∂x)-1/γ [(∂T)/(∂x)+T/ρ (∂ρ)/(∂x)]` (5)
Energy:
`(∂T)/(∂t)=-V (∂T)/(∂x)-(γ-1)[(∂V)/(∂x)+V (∂(ln(A)))/∂x]` (6)
Finite Difference Method
Having equipped with the appropriate governing equation, we shall now proceed to solve it using MacCormack’s explicit technique. MacCormack method follows a predictor-corrector method which can be understood better with the following flow steps
(a) Predictor Step
- Discretize the governing equation using the forward difference method.
- Update the primitive terms.
(b) Corrector Step.
- Discretize the governing equation using the backward difference method and using the updated primitive terms from the predictor step.
- Compute the average of the governing equation solved by the predictor and corrector method.
- Update the primitive terms with the average computed in the last step
Calculation of timesteps
We can see our governing equations are hyperbolic and we will be using explicit methods to solve it, and hence computing the timestep is very critical as explicit time-integration method is sensitive to timestep. Having a wrong timestep will blow up our solution in no time thereby making our solution unstable, illogical and unusable. Therefore, we shall be using the Courant number to calculate the timestep for a stable solution. We know, for a linear PDE, the Courant number needs to be less than 1 or at most equal to 1, to achieve a stable solution. The governing equation we have is a non-linear PDE, yet the above criterion, even though won’t represent the exact stability, can be used to give a picture of what could be a stable solution; in fact, upon simulation, we shall see this working out pretty good. Therefore, we use
`∆t=C[(∆x)/(alpha + V)]` (7)
where ∆t = time step
`alpha` = speed of sound
∆x = grid size
V = velocity
Area
The area for our problem was calculated using
`A = 1+2.2(x-1⋅5)^2`
Boundary Conditions
Conservation Form
Inflow Condition (at sub-sonic region i.e. i = 1)
`U_(1_((i=1) ) )=A_((i=1) )`
`U_(2_((i=1) ) )=2U_(2_((i=2 ) )``- U_(2_((i=3 ) )`
`U_(3_((i=1) ) )=U_(1_((i=1) ) ) (T_(1_((i=1) ) )/(γ-1)+γ/2 V_(1_((i=1) ))^2 )`
Outflow Condition (at super-sonic region i.e. i = n)
`U_(1_((i = n) ) )=2U_(1_((i=n - 1 ) )``- U_(1_((i=n - 2 ) )`
`U_(2_((i = n) ) )=2U_(2_((i=n - 1 ) )``- U_(2_((i=n - 2 ) )`
`U_(3_((i = n) ) )=2U_(3_((i=n - 1 ) )``- U_(3_((i=n - 2 ) )`
Non-Conservation Form
Inflow Condition (at sub-sonic region i.e. i = 1)
`ρ=1 - 0.3146x`
`T=1 - 0.2314x`
`V=(0.1+1.09x) √T`
Outflow Condition (at super-sonic region i.e. i = n)
`ρ_((i = n) )=2ρ_((i=n - 1 ))``- ρ_((i=n - 2 ) )`
`T_((i = n) )=2T_((i=n - 1 ))``- T_((i=n - 2 ) )`
`V_((i = n) )=2V_((i=n - 1 ))``- V_((i=n - 2 ) )`
Now, as we have prepared everything required to solve, we will jump directly into solving it and attempt to establish an understanding of the solution.
2. Main Program
clear all
close all
clc
% Parameters
n = 31; % Number of grids
throat = (n-1)/2; % Throat grid position
x = linspace(0,3,n); % Grid creation
dx = x(2)-x(1); % Grid size
nt = 1400; % Time
C = 0.5; % Courant Number
gamma = 1.4; % Specific heat ratio
% Area distribution
a = 1 + 2.2*(x-1.5).^2;
%% User Input
form_selector = input('Select the type of form you want to solve:\n1. Conservation\n2. Non-Conservation\n')
if form_selector == 1
[rho, vel, temp, dt, tc, mdot, rho_throat, temp_throat, p_throat, M_throat] = conservation_form(n, x, dx, nt, C, gamma, a, throat);
elseif form_selector == 2
[rho, vel, temp, dt, tc, mdot, rho_throat, temp_throat, p_throat, M_throat] = nonconservation_form(n, x, dx, nt, C, gamma, a, throat);
else
error('Please select from the options');
end
%% Comparison between Conservation, Non-Conservation and Exact solutions
% % Activate it to compare mass flow rate %
% % Deactivate user_input when this is active %
% % Deactivate all plots from each of the form solvers when this is active %
%
%
% % Exact Solution
% rho_ex = 0.634*ones(1,n);
% temp_ex = 0.833*ones(1,n);
% p_ex = 0.528*ones(1,n);
% M_ex = 1*ones(1,n);
%
% % Conservation
% [rho, vel, temp, dt, tc, mdot_c] = conservation_form(n, x, dx, nt, C, gamma, a, throat);
%
% %Nonconservation
% [rho, vel, temp, dt, tn, mdot_n] = nonconservation_form(n, x, dx, nt, C, gamma, a, throat);
%
% % figure(5)
% figure('units','normalized','outerposition',[0 0 1 1])
% plot(x,mdot_c,'r', 'LineWidth',2)
% hold on
% plot(x,mdot_n,'g','LineWidth',2)
% % axis([0 3 0.575 0.585])
% line([0 3],[0.579 0.579],'Color','k','LineWidth',2);
% title(sprintf('Mass flow variation comparison\n Number of grids = %d',n))
% lgd = legend('Conservation', 'Non-Conservation', 'Exact', 'Location', 'best')
% xlabel('x/L')
% ylabel('\rhoAV / \rhoA^*a_0')
% saveas(gcf,'rho comparison_61.png')
%% Form Functions
% Conservation Form
function [rho, vel, temp, dt, tc, mdot, rho_throat, temp_throat, p_throat, M_throat] = conservation_form(n, x, dx, nt, C, gamma, a, throat)
% Initial condition
tic;
for i = 1:n
if x(i) >= 0 & x(i) <= 0.5
rho(i) = 1;
temp(i) = 1;
vel(i) = 0.59/(rho(i)*a(i));
elseif x(i) > 0.5 & x(i) <= 1.5
rho(i) = 1 - (0.366*(x(i) - 0.5));
temp(i) = 1 - (0.167*(x(i) - 0.5));
vel(i) = 0.59/(rho(i)*a(i));
elseif x(i) > 1.5 & x(i) <= 3.5
rho(i) = 0.634 - (0.3879*(x(i) - 1.5));
temp(i) = 0.833 - (0.3507*(x(i) - 1.5));
vel(i) = 0.59/(rho(i)*a(i));
end
U2(i) = 0.59;
end
U1 = rho.*a;
U3 = rho.*((temp./(gamma - 1)) + (((gamma/2)*vel.^2))).*a;
F1 = U2;
F2 = ((U2.^2)./U1) + (((gamma - 1)/gamma) * (U3 - ((gamma/2)*((U2.^2)./U1))));
F3 = ((gamma*U2.*U3)./U1) - (((gamma*(gamma -1))/2) * (((U2).^3)./((U1).^2)));
% Time step calculation
sp_so = sqrt(temp); % Local speed of sound
dt=min(C*(dx./((sp_so + vel))));
for t = 1:nt
rho_in = rho;
vel_in = vel;
temp_in = temp;
U1_in = U1;
U2_in = U2;
U3_in = U3;
F1_in = F1;
F2_in = F2;
F3_in = F3;
% Predictor
for i = 2:n-1
dA_dx_p = (a(i+1) - a(i))/dx;
J2(i) = (1/gamma)*rho(i)*temp(i)*dA_dx_p;
% Governing equations
% Continuity
dU1_dt_p(i) = - (F1(i+1) - F1(i))/dx;
% Momentum
dU2_dt_p(i) = (- (F2(i+1) - F2(i))/dx) + J2(i);
% Energy
dU3_dt_p(i) = - (F3(i+1) - F3(i))/dx;
% Updating solutions
U1(i) = U1(i) + (dU1_dt_p(i)*dt);
U2(i) = U2(i) + (dU2_dt_p(i)*dt);
U3(i) = U3(i) + (dU3_dt_p(i)*dt);
% Updating primitive values
rho(i) = U1(i)/a(i);
temp(i) = (gamma - 1) * ((U3(i)/U1(i)) - ((gamma/2)*(U2(i)/U1(i))^2));
% Updating flux terms
F1 = U2;
F2 = ((U2.^2)./U1) + (((gamma - 1)/gamma) * (U3 - ((gamma/2)*((U2.^2)./U1))));
F3 = ((gamma*U2.*U3)./U1) - (((gamma*(gamma -1))/2) * (((U2).^3)./((U1).^2)));
end
% Corrector
for i = 2:n-1
dA_dx_c = (a(i) - a(i-1))/dx;
J2(i) = (1/gamma)*rho(i)*temp(i)*dA_dx_c;
% Governing equations
% Continuity
dU1_dt_c(i) = - (F1(i) - F1(i-1))/dx;
% Momentum
dU2_dt_c(i) = (- (F2(i) - F2(i-1))/dx) + J2(i);
% Energy
dU3_dt_c(i) = - (F3(i) - F3(i-1))/dx;
end
% Computing average time derivaive
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 Update
for i = 2:n-1
U1(i) = U1_in(i) + (dU1_dt(i)*dt);
U2(i) = U2_in(i) + (dU2_dt(i)*dt);
U3(i) = U3_in(i) + (dU3_dt(i)*dt);
end
% Updating primitive values
rho = U1./a;
vel = U2./U1;
temp = (gamma - 1) * ((U3./U1) - ((gamma/2)*(U2./U1).^ 2));
% Applying Boundary Condition
% Inlet
U1(1) = a(1);
U2(1) = (2*U2(2)) - (U2(3));
U3(1) = U1(1)*((1/(gamma - 1)) + ((gamma/2)*(vel(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);
% Defining pressure as the local static pressure
p = rho.*temp;
% Defining pressure as the local Mach number
M = vel./sqrt(temp);
% Local mass density
mdot = (rho.*a.*vel);
% Throat solution
rho_throat(n) = rho(throat);
vel_throat(n) = vel(throat);
temp_throat(n) = temp(throat);
p_throat(n) = p(throat);
M_throat(n) = M(throat);
% Mass flow rate along domain at different times
figure(1);
hold on;
if t == 50;
plot(x,mdot,'r','Linewidth',1.5);
elseif t == 100;
plot(x,mdot,'g','Linewidth',1.5);
elseif t == 200;
plot(x,mdot,'b','Linewidth',1.5);
elseif t == 400;
plot(x,mdot,'m','Linewidth',1.5);
elseif t == 800;
plot(x,mdot,'k','Linewidth',1.5);
plot(x(end),mdot(end),'rd','Linewidth',5);
end
end
tc = toc;
% Mass flow rate along domain at different times
figure(1)
title({'CONSERVATION FORM' ; 'Instantaneous distribution of mass flow rate along domain at different times'})
xlabel('x/x_0')
ylabel('\rhoAV / \rho_0 A^*a_0')
legend('50 timesteps','100 timesteps','200 timesteps','400 timesteps','800 timesteps','1400 timesteps','Location','best')
% Grid vs. Primitive terms
figure (2)
subplot(5,1,1)
plot(x,vel,'LineWidth',2)
title(sprintf('CONSERVATION FORM\n Grid size = %.3f | No. of grids = %d | No. of time steps = %d | Computational time = %.3fs', dx,n,t,tc))
xlabel('x/x_0')
ylabel('V / a_0')
grid on
subplot(5,1,2)
plot(x,rho,'LineWidth',2)
xlabel('x/x_0')
ylabel('\rho / \rho_0')
grid on
subplot(5,1,3)
plot(x,temp,'LineWidth',2)
xlabel('x/x_0')
ylabel('T / T_0')
grid on
subplot(5,1,4)
plot(x,M,'LineWidth',2)
xlabel('x/x_0')
ylabel('M / M_0')
grid on
subplot(5,1,5)
plot(x,p,'LineWidth',2)
xlabel('x/x_0')
ylabel('p / p_0')
grid on
% Iteration vs. Primitive terms
figure(3)
subplot(4,1,1)
plot(p_throat,'LineWidth',2)
title(sprintf('CONSERVATION FORM\n Grid size = %.3f | No. of grids = %d | No. of time steps = %d | Computational time = %.3fs', dx,n,t,tc))
xlabel('No. of time steps')
ylabel(sprintf('p / p_0'))
grid on
subplot(4,1,2)
plot(rho_throat,'LineWidth',2)
xlabel('No. of time steps')
ylabel('\rho / \rho_0')
grid on
subplot(4,1,3)
plot(temp_throat,'LineWidth',2)
xlabel('No. of time steps')
ylabel('T / T_0')
grid on
subplot(4,1,4)
plot(M_throat,'LineWidth',2)
xlabel('No. of time steps')
ylabel('M / M_0')
grid on
end
% Non-conservation form
function [rho, vel, temp, dt, tn, mdot, rho_throat, temp_throat, p_throat, M_throat] = nonconservation_form(n, x, dx, nt, C, gamma, a, throat)
% Initial condition
rho = 1 - 0.3146*x;
temp = 1 - 0.2314*x;
vel = (0.1 + 1.09*x).*temp.^(0.5);
% Time step calculation
dt = min(C.*(dx./(sqrt(temp) + vel)));
tic;
for t = 1:nt
rho_in = rho;
vel_in = vel;
temp_in = temp;
for i = 2:n-1
% Forward differencing the spatial terms
drho_dx_p = (rho(i+1)- rho(i))/dx;
dV_dx_p = (vel(i+1)- vel(i))/dx;
dlnA_dx_p = (log(a(i+1))- log(a(i)))/dx;
dT_dx_p = (temp(i+1)- temp(i))/dx;
% Governing Equation
% Continuity
drho_dt_p(i) = - (vel(i)*drho_dx_p) - (rho(i)*dV_dx_p) - (rho(i)*vel(i)*dlnA_dx_p);
% Momentum
dV_dt_p(i) = - (vel(i)*dV_dx_p) - ((1/gamma)*(dT_dx_p+((temp(i)/rho(i))*drho_dx_p)));
% Energy
dT_dt_p(i) = - (vel(i)*dT_dx_p) - (((gamma - 1)*temp(i)) * (dV_dx_p + (vel(i)*dlnA_dx_p)));
% Updating solutions
rho(i) = rho(i) + (drho_dt_p(i)*dt);
temp(i) = temp(i) + (dT_dt_p(i)*dt);
vel(i) = vel(i) + (dV_dt_p(i)*dt);
end
% Corrector method
for j = 2:n-1
% Backward differencing the spatial terms
drho_dx_c = (rho(j)- rho(j-1))/dx;
dV_dx_c = (vel(j)- vel(j-1))/dx;
dlnA_dx_c = (log(a(j))- log(a(j-1)))/dx;
dT_dx_c = (temp(j)- temp(j-1))/dx;
% Governing Equation
% Continuity
drho_dt_c(j) = - (rho(j)*dV_dx_c) - (rho(j)*vel(j)*dlnA_dx_c) - (vel(j)*drho_dx_c);
% Momentum
dV_dt_c(j) = - (vel(j)*dV_dx_c) - ((1/gamma)*(dT_dx_c+((temp(j)/rho(j))*drho_dx_c)));
% Energy
dT_dt_c(j) = - (vel(j)*dT_dx_c) - (((gamma - 1)*temp(j)) * (dV_dx_c + (vel(j)*dlnA_dx_c)));
end
% Computing average time derivaive
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 update
for k = 2:n-1
rho(k) = rho_in(k) + (drho_dt(k)*dt);
vel(k) = vel_in(k) + (dV_dt(k)*dt);
temp(k) = temp_in(k) + (dT_dt(k)*dt);
end
% Applying boundary conditions
% Inlet
vel(1) = (2*vel(2))-vel(3);
rho(1) = 1;
temp(1) = 1;
% Outlet
vel(n) = (2*vel(n-1)) - vel(n-2);
rho(n) = (2*rho(n-1)) - rho(n-2);
temp(n) = (2*temp(n-1)) - temp(n-2);
% Defining pressure as the local static pressure
p = rho.*temp;
% Defining pressure as the local Mach number
M = vel./sqrt(temp);
% Local mass density
mdot = (rho.*a.*vel);
% Throat solution
rho_throat(n) = rho(throat);
vel_throat(n) = vel(throat);
temp_throat(n) = temp(throat);
p_throat(n) = p(throat);
M_throat(n) = M(throat);
% Mass flow rate along domain at different times
figure(1);
hold on;
if t == 0;
plot(x,mdot,'-.','Linewidth',1.5);
elseif t == 50;
plot(x,mdot,'r','Linewidth',1.5);
elseif t == 100;
plot(x,mdot,'g','Linewidth',1.5);
elseif t == 200;
plot(x,mdot,'b','Linewidth',1.5);
elseif t == 400;
plot(x,mdot,'m','Linewidth',1.5);
elseif t == 800;
plot(x,mdot,'k','Linewidth',1.5);
plot(x(end),mdot(end),'rd','Linewidth',5);
end
end
tn = toc;
% Mass flow rate along domain at different times
figure(1)
title({'NON-CONSERVATION FORM' ; 'Instantaneous distribution of mass flow rate along domain at different times'})
xlabel('x/x_0')
ylabel('\rhoAV / \rho_0 A^*a_0')
legend('50 timesteps','100 timesteps','200 timesteps','400 timesteps','800 timesteps','1400 timesteps','Location','best')
% Grid vs. Primitive terms
figure(2)
subplot(5,1,1)
plot(x,vel,'LineWidth',2)
title(sprintf('NON-CONSERVATION FORM\n Grid size = %.3f | No. of grids = %d | No. of time steps = %d | Computational time = %.3fs', dx,n,t,tn))
xlabel('x/x_0')
ylabel('V / a_0')
grid on
subplot(5,1,2)
plot(x,rho,'LineWidth',2)
xlabel('x/x_0')
ylabel('\rho / \rho_0')
grid on
subplot(5,1,3)
plot(x,temp,'LineWidth',2)
xlabel('x/x_0')
ylabel('T / T_0')
grid on
subplot(5,1,4)
plot(x,M,'LineWidth',2)
xlabel('x/x_0')
ylabel('M / M_0')
grid on
subplot(5,1,5)
plot(x,p,'LineWidth',2)
xlabel('x/x_0')
ylabel('p / p_0')
grid on
% Iteration vs. Primitive terms
figure(3)
subplot(4,1,1)
plot(p_throat,'LineWidth',2)
title(sprintf('NON-CONSERVATION FORM\n Grid size = %.3f | No. of grids = %d | No. of time steps = %d | Computational time = %.3fs', dx,n,t,tn))
xlabel('No. of time steps')
ylabel(sprintf('p / p_0'))
grid on
subplot(4,1,2)
plot(rho_throat,'LineWidth',2)
xlabel('No. of time steps')
ylabel('\rho / \rho_0')
grid on
subplot(4,1,3)
plot(temp_throat,'LineWidth',2)
xlabel('No. of time steps')
ylabel('T / T_0')
grid on
subplot(4,1,4)
plot(M_throat,'LineWidth',2)
xlabel('No. of time steps')
ylabel('M / M_0')
grid on
end
3. Output
Conservation Form
For 31 grids
For 31 grids, even though the variation in the primitive terms through the domain length correctly represent the characteristics of a subsonic-supersonic flow, we can clearly notice it is not exactly smooth. Furthermore, from the timewise variation of primitive terms, we can see a lot of oscillation at the early time zone, and continuing towards the end of the time we see an asymptotically approached steady state. This is explained by the 'rubber-band effect'; this basically represents how during the early time zone the rubber band is stretched tightly and hence the primitive terms being driven at a rather stronger potential leading to rapid variations than during the end of the time where the rubber band is stretched lightly leading to a steady state or no drive, thereby having negligible variation.
The above plot is another way to study the mechanics of the timewise variation of the flow along the domain length. Here we have the mass flow through the domain length and see how it changes with time. After 50 timesteps, it can be seen that the oscillation is high. After 100 timesteps, the mass flow distribution has changed radically; the mass flow keeps flopping around inside the nozzle due to the transient variation of density, velocity and other primitive variables. But once it hits 200 timesteps, we can see the mass flow settling down and finally at around 800 timesteps, indicating the mass flow has converged to a constant, steady-state value throughout the nozzle.
For 61 grids
For 61 grids, we can evidently notice that it is a lot smoother flow than the one with 31 grids. Furthermore, from the timewise variation of primitive terms, we can see a similar oscillation trend at the early time zone, and continuing towards the end of the time we see an asymptotically approached steady state. But when compared to the same for 31 grids, we can notice that the oscillation with 61 grids is smoother than that of the 31 grids. Also, the damping is lower in case of 61 grids, hence taking longer to reach the steady-state; taking about 900 timesteps to reach steady state as compared to 800 for 31 grids.
We continue to study the mechanics of the timewise variation of the flow along the domain length. After 50 timesteps, it can be seen that the oscillation is high. After 100 timesteps, the mass flow distribution is still radically high and continues to stay comparatively highly varied until it hits 400 timesteps, which is when we can see the mass flow starting to settle down and finally at around 800 timesteps, the mass flow settles down indicating the mass flow has converged to a constant, steady-state value throughout the nozzle.
Non-Conservative Form
For 31 grids
For 31 grids, the variation of primitive terms through the domain length correctly represent the characteristics of a subsonic-supersonic flow, we can clearly notice the smooth flow throughout the domain. Furthermore, from the timewise variation of primitive terms, we can see a trend of oscillation at the early time zone, and continuing towards the end of the time we see an asymptotically approached steady state. This, as introduced earlier, is explained by the 'rubber-band effect'; this basically represents how during the early time zone the rubber band is stretched tightly and hence the primitive terms being driven at a rather stronger potential leading to rapid variations than during the end of the time where the rubber band is stretched lightly leading to a steady state or no drive, thereby having negligible variation; taking about 600 timesteps to reach the steady-state.
We perform the study of the mechanics of the timewise variation of the flow along the domain length by having the mass flow through the domain length and see how it changes with time. After 50 timesteps, it can be seen that the oscillation is high. After 100 timesteps, the mass flow distribution has changed radically; the mass flow keeps flopping around inside the nozzle due to the transient variation of density, velocity and other primitive variables. But once it hits 200 timesteps, we can see the mass flow settling down and finally at around 800 timesteps, the mass flow settles down indicating the mass flow has converged to a constant, steady-state value throughout the nozzle; and from the earlier plot we now know that it reaches steady state way before 800 timesteps.
For 61 grids
For 61 grids, we can notice no difference in the smoothness of the flow, even though we will later see the main significance of the finer mesh. Furthermore, from the timewise variation of primitive terms, we can see a similar oscillation trend at the early time zone, and continuing towards the end of the time we see an asymptotically approached steady state. But when compared to the same for 31 grids, we can notice that the oscillation with 61 grids is smoother than that of the 31 grids. Also, the damping is lower in case of 61 grids, hence taking longer to reach the steady-state; taking about 900 timesteps to reach steady state as compared to 600 for 31 grids.
We continue to study the mechanics of the timewise variation of the flow along the domain length with the mass flow through the domain length and see how it changes with time. After 50 timesteps, it can be seen that the oscillation is high. After 100 timesteps, the mass flow distribution is still radically high and continues to stay comparatively highly varied until it hits 400 timesteps, which is when we can see the mass flow starting to settle down and finally at around 800 timesteps, the mass flow settles down indicating the mass flow has converged to a constant, steady-state value throughout the nozzle.
Comparison of steady-mass flow rate between conservation and non-conservation form
From the flow, it can be deduced that the conservation form has closer to constant steady-mass flow distribution as compared to the non-conservation form, which has high variations specifically at the inflow and the outflow boundaries. Moreover, it is evident that the steady-mass flow is closer to the exact solution that the non-conservation form.
For 61 grids, the steady-mass flow distribution for the conservation and non-conservation form follows a similar trend except that the peaks are sharper for 61 grids with a lower peak point i.e. ranging between 0.585 and 0.586 for 61 grids as compared to the range of 0.594 to 0.596 for 31 grids. For the solution convergence, we get closer to the exact solution for both the conservation and non-conservation form as compared to that of the 31 grids.
4. Inference
From the table, the following observation can be deduced
a) Fastest convergence: Non-Convservation form with 61 grids
b) Best convergence (with least deviation from exact solution): Conservation form with 61 grids.
Moreover, the conservation form produces a better mass flow distribution; in other words, the conservation form performs the best job in conserving mass. The non-conservation form, however, produces a faster convergence due to less computational terms as compared to the conservation form. Conservation form requires marginally more work to than the non-conservation form.
To conclude, there is clearly no superiority of either form in terms of convergence and accuracy, both did a great job in representing the flow characteristics.
Note
The values will differ depending on the precision of your computer. The precision for my system was four decimal points, hence a clear difference in the developed solution can be visualized when compared to the three decimal point precision on the system used in the reference book.
Reference
Fig [a]; {Reference Book} John D. Anderson; Computational Fluid Dynamics the Basics With Applications [pg 285].