%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Developed and updated by Assoc. Prof. Dr. Eng. Kiril Shterev.
% All rights reserved by Assoc. Prof. Dr. Eng. Kiril Shterev.
% Match, 21st, 2022.
%
%
% Please cite my papers if you find this information useful:
%
% K. Shterev and S. Stefanov, Pressure based finite volume method
% for calculation of compressible viscous gas flows, Journal of
% Computational Physics 229 (2010) pp. 461-480, doi:10.1016/j.jcp.2009.09.042
%
% K. S. Shterev and S. K. Stefanov, A Parallelization of Finite Volume Method
% for Calculation of Gas Microflows by Domain Decomposition Methods, 7th
% Internnational Conference - Large-ScaleScientific Computations, Sozopol,
% Bulgaria, June 04-08, 2009. Lecture Notes in Computer Science Volume 5910,
% 2010, DOI: 10.1007/978-3-642-12535-5, SJR 0.295.
%
% Kiril S. Shterev, GPU implementation of algorithm SIMPLE-TS for calculation
% of unsteady, viscous, compressible and heat-conductive gas flows,
% URL: https://arxiv.org/abs/1802.04243, 2018.
%
%
% This is program to calculate One-dimentional unsteady, viscous,
% compressibnle, isotermal fluid flow in a duct.
%
% The program calculates the flow using two tipes of coefficients:
%
% 1. The coefficients used in SIMPLE-like methods: SIMPLE, SIMPLER, SIMPLEC
% and SIMPLEST-ANL.
% These coefficients do NOT substitude the density into unsteady term of
% equation for pressure and term (rho - rho_pr) exists.
% 2. The coefficients used in SIMPLE-TS.
% Here are substituted density into unsteady term of pressure equation
% using equation of state.
%
% The target is to compare pressure profiles obtained by two approaches
% at the end of first iteration.
%
% The results are used in paper:
% K. Shterev and S. Stefanov, Pressure based finite volume method
% for calculation of compressible viscous gas flows, Journal of
% Computational Physics 229 (2010) pp. 461-480, doi:10.1016/j.jcp.2009.09.042
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
FigureTitle='One_dimensional_flow_in_a_duct';
Figure_FontSize=17;
% Define the problem
Nx = 101; % Number of cells in OX direction
x_b = 0; % the begin coordinate of the mesh
x_e = 10; % the begin coordinate of the mesh
hx = x_b + (x_e - x_b) / (Nx - 1) % Step on OX direction
ht = 0.1; % Step on time
Nt = 1; % Number of time step to calculate
time_begin = 0; % start time of calculation
N_Iter = 20000; % Maximum nimbre of iterations to calculate one time step
p_in = 3; % Pressure on the inlet of the duct
p_out = 1; % Pressure on the outlet of the duct
% Set initial boundary conditions
% Initialize arrays
for i = 1:1:Nx
x_v(i) = x_b + hx * (i - 1);
u(i) = 0;
p(i) = p_out;
au0(i) = 0;
du(i) = 0;
apx(i) = 0;
bpx(i) = 0;
ap0_slash(i) = 0;
bp0_slash(i) = 0;
end
p(1) = p_in;
rho = p;
u_pseudo = u;
% Here are accepted for simplicity:
hxht = 1; % hx = ht
Dux = 1;
A = 1;
%
% Therefore:
au0(2) = (8.0 / 3.0) * Dux + 2 * hxht;
du(2) = A / au0(2);
apx(2) = 2.0 * (A / au0(2));
ap0_slash(2) = (2.0 * (A / au0(2)) + (A / (au0(2) - hxht))) * ht;
for i = 3:1:(Nx)
au0(i) = au0(2) - hxht;
du(i) = A / (au0(2) - hxht);
apx(i) = A / (au0(2) - hxht);
end
for i = 3:1:(Nx - 1)
ap0_slash(i) = 2.0 * (A / (au0(2) - hxht)) * ht;
end
% Calculate using coefficients without substitution of density woth pressure
% (standard SIMPLE-like methods).
for Iter = 1:1:N_Iter
residual = 0.0;
for i = 2:1:(Nx - 1)
p_old = p(i);
p(i) = ((apx(i) * p(i-1) + apx(i+1) * p(i+1)) * ht) / ap0_slash(i);
if(abs(p_old - p(i)) > residual)
residual = abs(p_old - p(i));
end
end
end
residual
% Plot data for SIMPLE-like methods
fig = figure('visible','on', 'Position', [1 1 1280 360]);
plot(x_v, p, 'k--', 'LineWidth', 3);
%% Calculate using coefficients with substitution of density with pressure
% (SIMPLE-TS method).
% Set initial boundary conditions
% Initiate arrays
for i = 1:1:Nx
p(i) = p_out;
ap0(i) = 0;
bp(i) = 0;
end
p(1) = p_in;
p_pr = p;
for i = 2:1:(Nx - 1)
ap0(i) = hx + (apx(i) + apx(i+1)) * ht;
bp(i) = p_pr(i) * hx;
end
% Start calculation on time step
for Iter = 1:1:N_Iter
residual = 0.0;
for i = 2:1:(Nx - 1)
p_old = p(i);
p(i) = ((apx(i) * p(i-1) + apx(i+1) * p(i+1)) * ht + bp(i)) / ap0(i);
if(abs(p_old - p(i)) > residual)
residual = abs(p_old - p(i));
end
end
end
hold on;
grid on;
residual
plot(x_v, p, 'k-', 'LineWidth', 3);
grid on;
% Format figure
set(gca,'FontSize', Figure_FontSize);
% Se labels
xlabel('x');
ylabel('pressure');
% Set labal
legend('SIMPLE-like methods (without substitution of density with pressure)', 'SIMPLE-TS (with substitution of density with pressure)');
% Save figures
saveas(fig, strcat(FigureTitle, '.fig'));
saveas(fig, strcat(FigureTitle, '.png'));
% Remove numbers fom axes, but do not remove grid lines
%set(gca, 'XTick', [1; 2; 3; 4; 5; 6; 7; 8; 9;], 'XTickLabel', []);
%set(gca, 'YTick', [1.5; 2.0; 2.5], 'YTickLabel', []);