## Numerical Solution of Quasi 1D Nozzle Flows: MacCormack's Technique

• I am going to write a program in MATLAB for solving Quasi one dimensional nozzle flows using MacCormack's Technique implementing both of the conservative and non-conservative forms of governing equations and perform the grid dependence test. Note, I will use the set of initial conditions provided in Anderson's textbook (chapter 7) along with the nozzle's shape and length details.

Problem Setup:
1. Form finite-difference expressions applicable to MacCormack's technique.
2. Chosen Courant number is 0.5.
3. Grid points used for project purpose is 30.
4. For grid dependence test another grid point used is 60.

Governing flow equations

1. Non-conservative form:

1. Continuity equation: (delrho)/(delt) = -rho(delV)/(delx) - rhoV(del(lnA))/(delx) - V(delrho)/(delx)

2. Momentum equation: (delV)/(delt) = -V(delV)/(delx) - 1/gamma((delT)/(delx) + T/rho (delrho)/(delx))

3. Energy equation: (delT)/(delt) = -V(delT)/(delx) - T(gamma-1)[(delV)/(delx) + V(del(lnA))/(delx)]

2. Conservative form:

1.
Continuity equation: (delU_1)/(delt) = - (delF_1)/(delx)

2. Momentum equation: (delU_2)/(delt) = - (delF_2)/(delx) + J_2

3.
Energy equation: (delU_3)/(delt) = - (delF_3)/(delx)

In the above equations,

Solution vectors are defined as:
U_1 = rhoA, U_2 = rhoAV, U_3 = rho(T/(gamma-1)+gamma/2V^2)A

Flux vectors are defined as:
F_1 = U_2, F_2 = (U_2^2)/(U_1) + 1/gammaU_1(gamma-1)[U_3/U_1-gamma/2(U_2/U_1)^2], F_3 = gamma(U_2U_3)/(U_1)-(gamma(gamma-1))/2(U_2^3)/U_1^2

Source term is defined as:
J_2 = 1/gammarhoT(delA)/(delx)

NOTE: All of the above variables are dimensionless.

MATLAB PROGRAM:

Function files

1. 'fd_dx.m':

% fd_dx - calculates forward difference of a input function

function f1 = fd_dx(f,dx,i)

f1 = (f(i+1) - f(i))/dx;

end

2. 'rd_dx.m':

% rd_dx - calculates reward differences of the input function

function f1 = rd_dx(f,dx,i)

f1 = (f(i) - f(i-1))/dx;

end

3. 'pred_corr.m':

% pred_corr - calculates predictor/corrector step

function [f1, f2, f3] = pred_corr(v, t, a, r, y, rho, V, T, dx, i)

term1 = -rho(i) * v;
term2 = -rho(i) * V(i) * a;
term3 = -V(i) * r;
term4 = -V(i) * v;
term5 = -y^-1 * (t + T(i)/rho(i)  * r);
term6 = -V(i) * t;
term7 = -(y-1) * T(i) * (v + V(i) * a);

f1 = term1 + term2  + term3;
f2 = term4 + term5;
f3 = term6 + term7;

end

4. 'conservative.m':

%conservative - Solves conservative form of govering equations

function [m, M, V, T, p, rho, temperature, density, pressure, mach_n] = conservative(x, dx, n, C, A, nt, y, throat)

%% Initial conditions
for i = 1 : length(x)
if x(i) >= 0 && 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

% Defining velocity and pressure
V = 0.59 ./ (rho .* A);
p = rho .* T;

% Defining solution vectors
U1 = rho .* A;
U2 = rho .* A .* V;
U3 = rho .* A .* (T / (y - 1) + y / 2 * V.^2);

% Iteration loop
for k = 1 : nt

% Copies of solution vectors
U1_old = U1;
U2_old = U2;
U3_old = U3;

% dt calculations
dt = min(C * dx ./ (sqrt(T) + V));

% Flux vectors
F1 = U2;
F2=(U2.^2 ./ U1) + ((y - 1) / y) * (U3 - (y * U2.^2) ./ (2 * U1));
F3 =(y * U2 .* U3 ./ U1) - (y * (y - 1) / 2) * (U2.^3 ./ U1.^2);

% Predictor step
for i = 2 : n - 1
J2(i) = 1 / y * rho(i) .* T(i) * fd_dx(A, dx, i);
dU1dt_p(i) = -fd_dx(F1, dx, i);
dU2dt_p(i) = -fd_dx(F2, dx, i) + J2(i);
dU3dt_p(i) = -fd_dx(F3, dx, i);

% Updating solution vectors
U1(i)= U1(i) + dU1dt_p(i) * dt;
U2(i)= U2(i) + dU2dt_p(i) * dt;
U3(i)= U3(i) + dU3dt_p(i) * dt;
end

% Updation of primitive variables and flux vectors
rho = U1 ./ A;
T = (y - 1) * (U3 ./ U1 - (y / 2) * (U2 ./ U1).^2);
p = rho .* T;
F1 = U2;
F2=(U2.^2 ./ U1) + ((y - 1) / y) * (U3 - (y * U2.^2) ./ (2 * U1));
F3 =(y * U2 .* U3 ./ U1) - (y * (y - 1) / 2) * (U2.^3 ./ U1.^2);

% Corrector step
for j = 2 : n - 1
J2(j) = 1 / y * rho(j) .* T(j) * rd_dx(A, dx, j);       % Updated source term
dU1dt_c(j) = -rd_dx(F1, dx, j);
dU2dt_c(j) = -rd_dx(F2, dx, j) + J2(j);
dU3dt_c(j) = -rd_dx(F3, dx, j);
end

% Averaging corrector step and predictor step
dU1dt_a = 0.5 * (dU1dt_p + dU1dt_c);
dU2dt_a = 0.5 * (dU2dt_p + dU2dt_c);
dU3dt_a = 0.5 * (dU3dt_p + dU3dt_c);

% Updating solution vectors after averaging
for m = 2 : n - 1
U1(m) = U1_old(m) + dU1dt_a(m) * dt;
U2(m) = U2_old(m) + dU2dt_a(m) * dt;
U3(m) = U3_old(m) + dU3dt_a(m) * dt;
end

% Boundary conditions
% Inlet
U2(1) = 2*U2(2) - U2(3);

% 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);

% Updation of flow field variables
rho = U1 ./ A;
V = U2 ./ U1;
T = (y - 1) * (U3 ./ U1 - (y / 2) * (U2 ./ U1).^2);
p = rho .* T;

% Defining mass flow and mach number
m = rho .* A .* V;
M = V ./ sqrt(T);

% Capturing the values at throat for each iterations
mach_n(k) = M(throat);
pressure(k) = p(throat);
density(k) = rho(throat);
temperature(k) = T(throat);

% Ploting mass flow rate at different times
if n == 31

v = figure(1);
hold on
if k == 100
plot(x,m, 'b', 'linewidth', 2)
elseif k  == 200
plot(x, m, 'r', 'linewidth', 2)
elseif k  == 300
plot(x, m, 'm', 'linewidth', 2)
elseif k  == 550
plot(x, m, 'k', 'linewidth', 2)
title('Non - Dimensional Mass Flow Distribution at different Time Steps', 'FontAngle', 'italic', 'FontSize', 13, 'FontWeight', 'bold')
legend('100Deltat', '200Deltat', '300Deltat', '550Deltat')
ylim([0.5 0.7])
xlabel('x/L', 'FontWeight', 'bold', 'FontSize', 11)
ylabel('rhoVA/rho_{o}V_{o}A_{0}', 'FontWeight', 'bold', 'FontSize', 11)
grid minor
saveas(v, 'cons_mass_flow.png')
end

end

end

end

5. 'non_conservative.m':

%nonconservative - Solves non-conservative form of govering equations

function [m, M, V, T, p, rho, temperature, density, pressure, mach_n] = non_conservative(x, dx, n, C, A, nt, y, throat)

%% Initial conditions for the flow field variables
rho = 1 - 0.3146 * x;                                                   % Density
T = 1 - 0.2314 * x;                                                      % Temprature
V = (0.1 + 1.09 * x) .* T.^0.5;                                     % Velocity

for k = 1 : nt

% Copy of Field Variables
rho_old = rho;
V_old = V;
T_old = T;

% Time step
dt = min(C * dx ./ (T.^0.5 + V));

% Predictor step
for i = 2:n-1
vf = fd_dx(V,dx,i); tf = fd_dx(T,dx,i); af = fd_dx(log(A),dx,i); rf = fd_dx(rho,dx,i);
[drdt_p(i), dVdt_p(i), dTdt_p(i)] = pred_corr(vf, tf, af, rf, y, rho, V, T, dx, i);

%solution update
rho(i) = rho(i) + drdt_p(i)*dt;
V(i) = V(i) + dVdt_p(i)*dt;
T(i) = T(i) + dTdt_p(i)*dt;
end

% Corrector step
for j = 2:n-1
vr = rd_dx(V,dx,j); tr = rd_dx(T,dx,j); ar = rd_dx(log(A),dx,j); rr = rd_dx(rho,dx,j);
[drdt_c(j), dVdt_c(j), dTdt_c(j)] = pred_corr(vr, tr, ar, rr, y, rho, V, T, dx, j);
end

% Taking average of corrector and predictor step
drdt_a = 0.5 * (drdt_c + drdt_p);
dVdt_a = 0.5 * (dVdt_c + dVdt_p);
dTdt_a = 0.5 * (dTdt_c + dTdt_p);

% Flow field variables at next time step
for s = 2:n-1
rho(s) = rho_old(s) + drdt_a(s) * dt;
V(s) = V_old(s) + dVdt_a(s) * dt;
T(s) = T_old(s) + dTdt_a(s) * dt;
end

% Boundary Conditions
% subsonic inflow
V(1) = 2*V(2) - V(3);
% supersonic outflow
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);

% Defining pressure, mass flow and mach number
p = rho .* T;
m = rho .* A .* V;
M = V ./ sqrt(T);

% Capturing values at throat for each iterations
mach_n(k) = M(throat);
pressure(k) = p(throat);
density(k) = rho(throat);
temperature(k) = T(throat);

% Mass flow at different time steps
if n == 31

v = figure(2);
hold on
if k  == 50
plot(x, m, 'r', 'linewidth', 2)
elseif k  == 100
plot(x, m, 'm', 'linewidth', 2)
elseif k  == 150
plot(x, m, 'c', 'linewidth', 2)
elseif k  == 700
plot(x, m, 'k', 'linewidth', 2)
title('Non - Dimensional Mass Flow Distribution at different Time Steps', 'FontAngle', 'italic', 'FontSize', 13, 'FontWeight', 'bold')
legend('50Deltat', '100Deltat', '150Deltat', '700Deltat')
xlabel('x/L', 'FontWeight', 'bold', 'FontSize', 11)
ylabel('rhoVA/rho_{o}V_{o}A_{0}', 'FontWeight', 'bold', 'FontSize', 11)
grid minor
saveas(v, 'non_cons_mass_flow.png')
end

end

end

end

The driver script 'quasi_one_D_nozzle_flows.m' is as follows

% Quasi 1D Nozzle Flows

clear all
close all
clc

%% Setup
L = 3;                                               %Length of the nozzle
n = input('Choose number of grid points: ') + 1;     % Grid points - choose it manually
C =0.5;                                              % Courant number
x = linspace(0,L,n);                                 % Creating 1D grid
dx = L/(n-1);                                        % Adjecent grid spacing
y = 1.4;                                             % Gamma
nt = 1400;                                           % Iterations / Total time steps
throat = (n-1)/2 + 1;                                % Throat point
A = 1 + 2.2 * (x - 1.5).^2;                          % Area

%% Conservative form
tic
[Cm, CM, CV, CT, Cp, Crho, Ctemperature, Cdensity, Cpressure, Cmach_n] = conservative(x, dx, n, C, A, nt, y, throat);
cons = toc;
fprintf('nSolution time for conservative form is %0.3g secondsn', cons)

%% Non-Conservative form
tic
[Nm, NM, NV, NT, Np, Nrho, Ntemperature, Ndensity, Npressure, Nmach_n] = non_conservative(x, dx, n, C, A, nt, y, throat);
non_cons = toc;
fprintf('Solution time for non - conservative form is %0.3g secondsn', non_cons)

%% Ploting the above data
u = figure(3);
subplot(411)
hold on
plot(Cmach_n, 'b', 'linewidth', 2)
plot(Nmach_n, '--k', 'linewidth', 2)
s = legend('Conservative Form', 'Non-Conservative Form');
set(s, 'Location', 'northeastoutside')
title('Timewise Variation at the Nozzle Throat', 'FontAngle', 'italic', 'FontSize', 13, 'FontWeight', 'bold')
ylabel('M', 'FontSize', 11, 'FontWeight', 'bold')
grid minor

subplot(412)
hold on
plot(Cpressure, 'm', 'linewidth', 2)
plot(Npressure, '--k', 'linewidth', 2)
e = legend('Conservative Form', 'Non-Conservative Form');
set(e, 'Location', 'northeastoutside')
ylabel('p/p_{0}', 'FontSize', 11, 'FontWeight', 'bold')
grid minor

subplot(413)
hold on
plot(Cdensity, 'r', 'linewidth', 2)
plot(Ndensity, '--k', 'linewidth', 2)
r = legend('Conservative Form', 'Non-Conservative Form');
set(r, 'Location', 'northeastoutside')
ylabel('rho/rho_{0}', 'FontSize', 11, 'FontWeight', 'bold')
grid minor

subplot(414)
hold on
plot(Ctemperature, 'g', 'linewidth', 2)
plot(Ntemperature, '--k', 'linewidth', 2)
t = legend('Conservative Form', 'Non-Conservative Form');
set(t, 'Location', 'northeastoutside')
ylabel('T/T_{o}', 'FontSize', 11, 'FontWeight', 'bold')
xlabel('Number of Iterations', 'FontSize', 11, 'FontWeight', 'bold')
grid minor
saveas(u,'nozzle_throat.png')

b = figure(4);
subplot(411)
hold on
plot(x, CM, 'b', 'linewidth', 2)
plot(x, NM, '-k+', 'linewidth', 1, 'markersize', 5, 'markerfacecolor', 'k')
o = legend('Conservative Form', 'Non-Conservative Form');
set(o, 'Location', 'northeastoutside')
grid minor
hold on
title('Qualitative Aspects of Quasi - 1D Nozzle Flow', 'FontSize', 13, 'FontWeight', 'bold', 'FontAngle', 'italic')
ylabel('Mach number', 'FontWeight', 'bold', 'FontSize', 11)

subplot(412)
hold on
plot(x, CT, 'g', 'linewidth', 2)
plot(x, NT, '-k+', 'linewidth', 1, 'markersize', 5, 'markerfacecolor', 'k')
N = legend('Conservative Form', 'Non-Conservative Form');
set(N, 'Location', 'northeastoutside')
ylabel('Temperature ratio', 'FontWeight', 'bold', 'FontSize', 11)
grid minor

subplot(413)
hold on
plot(x, Crho, 'r', 'linewidth', 2)
plot(x, Nrho, '-k+', 'linewidth', 1, 'markersize', 5, 'markerfacecolor', 'k')
D = legend('Conservative Form', 'Non-Conservative Form');
set(D, 'Location', 'northeastoutside')
ylabel('Density ratio', 'FontWeight', 'bold', 'FontSize', 11)
grid minor

subplot(414)
hold on
plot(x, Cp, 'm', 'linewidth', 2)
plot(x, Np, '-k+', 'linewidth', 1, 'markersize', 5, 'markerfacecolor', 'k')
F = legend('Conservative Form', 'Non-Conservative Form');
set(F, 'Location', 'northeastoutside')
ylabel('Pressure ratio', 'FontWeight', 'bold', 'FontSize', 11)
xlabel('Non-Dimensional Length of the Nozzle', 'FontWeight', 'bold', 'FontSize', 11)
grid minor
saveas(b, 'qualitative_aspects.png')

q = figure(5);
hold on
plot(x, Cm, 'b', 'linewidth', 2)
plot(x, Nm, 'r', 'linewidth', 2)
line([0 3], [0.579 0.579], 'color', 'm', 'linewidth', 2.5)
Z = legend('Conservative Form', 'Non-Conservative Form', 'Exact Solution');
set(Z, 'Location', 'north')
grid minor
xlabel('x/L', 'FontWeight', 'bold', 'FontSize', 11)
ylabel('rhoVA/rho_{o}V_{o}A_{0}', 'FontWeight', 'bold', 'FontSize', 11)
title('Non-Dimensional Mass Flow Distribution', 'FontWeight', 'bold', 'FontSize', 13, 'FontAngle', 'italic')
saveas(q, 'mass_flow_comparission.png')

Results

Comparision of steady-state results; conservative vs. non-conservative form -

rho/(rho_o)            T/(T_o)          p/(p_o)           M

Exact analytical solution:
0.634        0.833       0.528      1.000

Non-conservative form:       0.639        0.836       0.534       0.999
(30 points)

Non-conservative form:       0.638        0.835       0.533      1.000
(60 points)

Conservative form:               0.650        0.840       0.546      0.982
(30 points)

Conservative form:               0.638        0.835       0.533      0.999

(60 points)

Command window display -

For 30 points:

Solution time for conservative form is 22.9 seconds
Solution time for non - conservative form is 34.7 seconds

For 60 points:

Solution time for conservative form is 53.7 seconds
Solution time for non - conservative form is 87.9 seconds

Plots -

1. For 30 points:

1. Qualitative aspects: 2. Timewise variation at the throat: 3. Mass flow distribution at different times: Conservative form: 4. Mass flow distribution at different times: Conservative form: 5. Mass flow rate comparison: 2. For 60 points:

6.
Mass flow rate comparison plot - Conclusions

1. For 30 points, as seen from the above plots 4 and 5 the minimum number of cycles for convergence in approximately 550 for conservative form and approximately 700 for non-conservative form.

2. It totally depends on what circumstances you are performing the simulation. If a very minor error is not going to create an issue to your task then 30 grid points are enough else you can increase the number of grid points for better yield of the solution.

3. In terms of accuracy, there is no clear supremacy in either of the forms. Though the results for each of the forms are certainly satisfactory.

4. As from the above plot 5 and 6 its seen that conservative form results overall are pretty much closer to the exact analytical solution (0.579). Whereas in non-conservative form there are spurious oscillations observed at inflow and outflow boundaries.

5. The conservative form fetches a better mass flow distribution. It simply does a better of conserving mass which is obviously guessed from the name of the form itself.

6. Although it takes more efforts in computing the conservative form, it eventually produces results which are faster than the non-conservative form.

### Planck's Radiation Law | Quantum Physics Rohan Shah · 2019-08-24 13:25:40

I was going through Max Planck's radiation law for a couple of days and the way that Planck had derived his beautiful equations just by assuming that the energy is not continuous rather it is an integral multiple of h (Planck constant) just amazed me. I couldn't re Read more

### Determining minimum cushion pressure to break the thickness of ice field | Newton - Raphson Method Rohan Shah · 2018-10-16 01:16:41

I am going to write a code in python for determining the minimum cushion pressure required to break the given thickness of ice by implementing the Newton-Raphson technique. I will also be running a setup to find the optimum relaxation factor for the thickness of 0.6 fee Read more

### Constraint Minimization Problem Rohan Shah · 2018-08-06 07:08:12

I am going to minimize a non-linear function by using Lagrange Multipliers. Problem statement:Minimize the function 5-(x-2)^2 -2(y-1)^2 subject to x + 4y - 3 = 0. Solution:    L = f - lambdag,    where  f = 5-(x-2)^2 -2(y-1 Read more

### Data Visualizer Tool Rohan Shah · 2018-06-17 21:44:39

I am going to write a program in python to parse information from the CONVERGE output file (IC engine simulation results) and use that data for basic performance calculations for engine and generating appropriate plots. Aim:1. Plot the data for any columns Read more

### A 2R Robotic Arm simulator Rohan Shah · 2018-06-16 12:40:49

I am going to write a program in python to simulate a 2R Robotic Arm. I've used online GIF maker to stitch all of my figures and generate the GIF. Note that, I have just used it to stitch the figures, everything else is done by the program. All the figures (200) have be Read more

### Demonstrating the Curve Fitting Rohan Shah · 2018-06-15 02:03:07

I am going to write a program in python to demonstrate the curve fit on the data set provided by increasing the degree of polynomials namely first degree, third degree and fourth degree. Program: # Program to demonstrate Curve Fit by increasing the degree of the pol Read more

### Otto cycle - Calculating the efficiency of the cycle and generating P - V Diagram Rohan Shah · 2018-06-14 14:31:48

I am going to write a program in python to simulate the Otto cycle. I will be calculating the thermal efficiency of the engine and be generating the P - V diagram.   Problem Setup and geometric parameters:1. Gamma is set to 1.4.2. The pressure at state p Read more

### Calculating Drag Force over a Bicycle Rohan Shah · 2018-06-14 09:43:42

I am going to write a program in python which will calculate the drag force over a bicycle and plot it against velocity and drag coefficients.   Problem Setup:1. Velocities to be used are 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12 and 132. Drag Coefficients to be u Read more

### Spectral radius Rohan Shah · 2018-05-12 21:47:10

I am going to code in MATLAB for solving linear system Ax = b` using iterative methods namely, Jacobi, Gauss-Seidel, and Successive over-relaxation. I will be calculating the eigenvalues and spectral radius of iteration matrix. The aim of this project is to n Read more