Numerical Solution of Quasi 1D Nozzle Flows: Conservative Form
Problem Setup:
- Form finite-difference expressions applicable to MacCormack's technique [Conservative Form].
- Chosen Courant number is 0.5.
- Grid points used for project purpose is 30.
- For grid dependence test another grid point used is 90.
Matlab Program:
%%=========Conservative form of governing equations=========%%
clear all
close all
clc
%% PROBLEM INITIALIZATION
tic
Lx = 3; % Length along x axis
nx = 31; % Number of Nodes along x axis
dx = Lx/(nx-1); % Grid Spacing
gamma = 1.4; % Gamma for Air
nt = 1400; % Number of time step
C = 0.5; % Courant Number
throat = (nx-1)/2 + 1; % Throat point
%% Grid Creation
x = linspace(0,Lx,nx);
%% Initial conditions
a = 1 + 2.2*(x - 1.5).^2;
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 / (gamma - 1) + gamma / 2 * v.^2);
dt = min(C * dx ./ (T.^0.5 + v));
%% Copy for Initial Profile
rho_initial = rho;
v_initial = v;
T_initial = T;
%% Solver
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) + ((gamma - 1) / gamma) * (U3 - (gamma * U2.^2) ./ (2 * U1));
F3 =(gamma * U2 .* U3 ./ U1) - (gamma * (gamma - 1) / 2) * (U2.^3 ./ U1.^2);
% Predictor step
for i = 2 : nx - 1
J2(i) = 1 / gamma * rho(i) .* T(i) * (a(i+1)-a(i))/dx;
dU1dt_p(i) = -(F1(i+1)-F1(i))/dx;
dU2dt_p(i) = -(F2(i+1)-F2(i))/dx + J2(i);
dU3dt_p(i) = -(F3(i+1)-F3(i))/dx;
% 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 = (gamma - 1) * (U3 ./ U1 - (gamma / 2) * (U2 ./ U1).^2);
p = rho .* T;
F1 = U2;
F2=(U2.^2 ./ U1) + ((gamma - 1) / gamma) * (U3 - (gamma * U2.^2) ./ (2 * U1));
F3 =(gamma * U2 .* U3 ./ U1) - (gamma * (gamma - 1) / 2) * (U2.^3 ./ U1.^2);
% Corrector step
for j = 2 : nx - 1
J2(j) = 1 / gamma * rho(j) .* T(j) * (a(j)-a(j-1))/dx; % Updated source term
dU1dt_c(j) = -(F1(j)-F1(j-1))/dx;
dU2dt_c(j) = -(F2(j)-F2(j-1))/dx + J2(j);
dU3dt_c(j) = -(F3(j)-F3(j-1))/dx;
end
% Averaging corrector step and predictor step
dU1dt_avg = 0.5 * (dU1dt_p + dU1dt_c);
dU2dt_avg = 0.5 * (dU2dt_p + dU2dt_c);
dU3dt_avg = 0.5 * (dU3dt_p + dU3dt_c);
% Updating solution vectors after averaging
for m = 2 : nx - 1
U1(m) = U1_old(m) + dU1dt_avg(m) * dt;
U2(m) = U2_old(m) + dU2dt_avg(m) * dt;
U3(m) = U3_old(m) + dU3dt_avg(m) * dt;
end
% Boundary conditions
% Inlet
U2(1) = 2*U2(2) - U2(3);
% Outlet
U1(nx) = 2*U1(nx-1) - U1(nx-2);
U2(nx) = 2*U2(nx-1) - U2(nx-2);
U3(nx) = 2*U3(nx-1) - U3(nx-2);
% Updation of flow field variables
rho = U1 ./ a;
v = U2 ./ U1;
T = (gamma - 1) * (U3 ./ U1 - (gamma / 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);
% Mass flow at different time steps
figure(1);
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 == 300
plot(x, m, '+r', 'linewidth', 2)
elseif k == 500
plot(x, m, 'g', 'linewidth', 2)
elseif k == 700
plot(x, m, '--k', 'linewidth', 2)
elseif k == 1400
plot(x, m, 'ob', 'linewidth', 2)
title(['Non - Dimensional Mass Flow Distribution at different Time Steps [Conservative]-' num2str(nx-1) 'points'], 'FontAngle', 'italic', 'FontSize', 13, 'FontWeight', 'bold')
legend('50\Deltat', '100\Deltat', '150\Deltat', '300\Deltat', '500\Deltat', '700\Deltat', '1400\Deltat')
xlabel('x/L', 'FontWeight', 'bold', 'FontSize', 11)
ylabel('\rhoVA/\rho_{o}V_{o}A_{0}', 'FontWeight', 'bold', 'FontSize', 11)
grid minor
end
end
%% Ploting the above data
figure(2)
subplot(411)
hold on
plot(mach_n, '--r', 'linewidth', 2)
% s = legend('Conservative Form');
% set(s, 'Location', 'northeastoutside')
title(['Timewise Variation at the Nozzle Throat [Conservative]-' num2str(nx-1) 'points'], 'FontAngle', 'italic', 'FontSize', 13, 'FontWeight', 'bold')
ylabel('M', 'FontSize', 11, 'FontWeight', 'bold')
grid minor
subplot(412)
hold on
plot(pressure, '--r', 'linewidth', 2)
% e = legend('Conservative Form');
% set(e, 'Location', 'northeastoutside')
ylabel('p/p_{0}', 'FontSize', 11, 'FontWeight', 'bold')
grid minor
subplot(413)
hold on
plot(density, '--r', 'linewidth', 2)
% r = legend('Conservative Form');
% set(r, 'Location', 'northeastoutside')
ylabel('\rho/\rho_{0}', 'FontSize', 11, 'FontWeight', 'bold')
grid minor
subplot(414)
hold on
plot(temperature, '--r', 'linewidth', 2)
% t = legend('Conservative Form');
% set(t, 'Location', 'northeastoutside')
ylabel('T/T_{o}', 'FontSize', 11, 'FontWeight', 'bold')
xlabel('Number of Iterations', 'FontSize', 11, 'FontWeight', 'bold')
grid minor
figure(3)
subplot(411)
hold on
plot(x, M, '-b+', 'linewidth', 1, 'markersize', 5, 'markerfacecolor', 'k')
% o = legend('Conservative Form');
% set(o, 'Location', 'northeastoutside')
grid minor
title(['Qualitative Aspects of Quasi - 1D Nozzle Flow [Conservative]-' num2str(nx-1) 'points'], 'FontSize', 13, 'FontWeight', 'bold', 'FontAngle', 'italic')
ylabel('Mach number', 'FontWeight', 'bold', 'FontSize', 11)
subplot(412)
hold on
plot(x, T, '-b+', 'linewidth', 1, 'markersize', 5, 'markerfacecolor', 'k')
% N = legend('Conservative Form');
% set(N, 'Location', 'northeastoutside')
ylabel('Temperature ratio', 'FontWeight', 'bold', 'FontSize', 11)
grid minor
subplot(413)
hold on
plot(x, rho, '-b+', 'linewidth', 1, 'markersize', 5, 'markerfacecolor', 'k')
% D = legend('Conservative Form');
% set(D, 'Location', 'northeastoutside')
ylabel('Density ratio', 'FontWeight', 'bold', 'FontSize', 11)
grid minor
subplot(414)
hold on
plot(x, p, '-b+', 'linewidth', 1, 'markersize', 5, 'markerfacecolor', 'k')
% F = legend('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
figure(4);
hold on
plot(x, m, 'r', 'linewidth', 2)
line([0 3], [0.579 0.579], 'color', 'k', 'linewidth', 2.5)
Z = legend('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-' num2str(nx-1) 'points'], 'FontWeight', 'bold', 'FontSize', 13, 'FontAngle', 'italic')
non_cons = toc;
fprintf('Solution time for conservative form is %0.3g seconds\n', non_cons)Computation Time:
For 30 points:
Solution time for conservative form is 2.6 seconds.
For 60 points:
Solution time for conservative form is 3.1 seconds.
For 90 points:
Solution time for conservative form is 3.8 seconds.
Plots:
For 30 Points:
Mass flow distribution at different times: Conservative form:
Timewise variation at the throat:
Qualitative aspects:
Mass flow rate comparison:
For 90 Points:
Mass flow distribution at different times: Conservative form:
Timewise variation at the throat:
Qualitative aspects:
Mass flow rate comparison:
Results:
Conclusion:
- Time of computation for 30 points < Time of computation for 60 points < Time of computation for 90 points.
- No much variation in results for 30, 60 and 90 points respectively.