function grains (grain_arrangement_file, Nx, Ny, nu, accel, t_max)
% grains.m
%
% Created by Jeff Cunningham while on sabbatical from USF and visiting
% CSIRO in Floreat, Western Australia, 2013-2014
%
% Uses lattice Boltzmann modeling to simulate 2D fluid flow through a
% porous medium. The user specifies the geometry of the porous medium by
% providing the file grain_arrangement_file.
%
% Flow is left to right.
% Uses periodic flow, i.e., flow exiting the right boundary re-enters on the
% left boundary.
%
% Some parts of this lattice Boltzmann code are modeled after codes posted
% on-line by Jonas Latt and co-workers. I also have received assistance
% from Shadab Anwar and Mike Sukop as I have developed this code.
%
% grain_arrangement_file is the name of a file containing data that specify
% the geometry of the porous domain. A 1 indicates that the node is
% occupied by a grain; a 0 indicates that the node is fluid.
% Nx and Ny are the size of the domain in lattice units. Nx and Ny should
% agree with the data in grain_arrangement_file.
% nu is the fluid viscosity in LBM units (lu^2/ts).
% Recommend nu = 1/6.
% accel is the acceleration acting on the fluid, units lu/(ts^2).
% accel = dP/(rho*L) if dP/L is an applied pressure gradient.
% Recommend a low value of accel, perhaps 1.E-7, to ensure low enough
% fluid velocities. This code is intended to simulate laminar flow.
% t_max is the number of time steps for the simulation to run.
% t_max must be high enough that the system reaches steady state.
% As a rough guideline, use t_max = Nx*Ny/10 .
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% Geometry of the problem %%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
fid = fopen(grain_arrangement_file, 'r');
pore_wall = fscanf(fid, '%i', [Nx, Ny]);
fclose(fid);
is_fluid = 1 - pore_wall; % if it's not a pore wall, then it is fluid
[y, x] = meshgrid(1:Ny, 1:Nx);
figure(1);
%imagesc(pore_wall');
surf(pore_wall'); shading flat;
axis('equal');
axis([1 Nx 1 Ny]);
xlabel('x coordinate (lattice units)', 'FontSize', 16);
ylabel('y coordinate (lattice units)', 'FontSize', 16);
title('Map of the porous medium', 'FontSize', 16);
drawnow;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% Set some LBM variables %%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
density = 1; % density in LBM space
omega = 1. / (3*nu + 1./2.); % relaxation parameter
tau = 1/omega;
u_avg = (accel *(Ny/2)^2)/(3*nu);
Re = (Ny/2)*u_avg/nu;
disp([nu omega tau accel u_avg Re]);
green_light = true; % parameter to tell if we should proceed
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% D2Q9 lattice constants %%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Use this numbering scheme for the directions:
% 6 2 5
% 3 9 1
% 7 4 8
%
w = [1/9, 1/9, 1/9, 1/9, 1/36, 1/36, 1/36, 1/36, 4/9];
opp = [ 3, 4, 1, 2, 7, 8, 5, 6, 9];
%e_x = [ 1, 0, -1, 0, 1, -1, -1, 1, 0];
%e_y = [ 0, 1, 0, -1, 1, 1, -1, -1, 0];
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% Pre-allocate the f matrices %%%%
%%%% This makes Matlab happy %%%%
%%%% Some of the matrices don't need pre-allocation %%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
f = zeros(9, Nx, Ny);
f_temp = zeros(9, Nx, Ny);
%f_eq = zeros(9, Nx, Ny);
%u_x = zeros(Nx, Ny);
%u_y = zeros(Nx, Ny);
%rho = zeros(Nx, Ny);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% Initialize the distribution function f %%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for k = 1:9
f(k, 1:Nx, 1:Ny) = w(k) * density;
end
initial_fluid_mass = sum(sum( is_fluid .* reshape(sum(f),Nx,Ny) )); % Track the fluid mass over time.
disp(initial_fluid_mass); % Verifies code is running properly.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% Begin main loop %%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
t = 0;
while (green_light) && (t < t_max)
t = t+1;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% Calculate macroscopic variables %%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[rho, u_x, u_y] = calc_macroscopic_variables_5 ( Nx, Ny, f );
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% Every 50 time steps, check the mass balance %%%%
%%%% Also display the velocity field %%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if (mod(t, 50) == 0)
green_light = check_mass_balance(t, initial_fluid_mass, rho, is_fluid);
u = reshape(sqrt(u_x.^2 + u_y.^2), Nx, Ny);
figure(2);
%imagesc(u');
%surf(u'); colorbar;
pcolor(u');
axis('equal'); axis([1 Nx 1 Ny]); shading interp;
drawnow;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% Calculate equilibrium distribution function %%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
f_eq = calc_f_eq ( Nx, Ny, is_fluid, u_x, u_y, rho, accel, tau);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% Collision or relaxation %%%%
%%%% Use standard BGK method %%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for i = 1:Nx
for j = 1:Ny
if is_fluid(i,j) % fluid nodes relax to equilibrium
for k = 1:9
f_temp(k,i,j) = f(k,i,j) - (f(k,i,j) - f_eq(k,i,j))/tau;
end
else % solid nodes use bounce-back
for k = 1:9
f_temp(k,i,j) = f(opp(k),i,j);
end
end
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% Streaming step %%%%
%%%% This code assumes periodic boundaries. %%%%
%%%% Fluid exiting on the right re-enters on the left. %%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
f = stream_2(Nx, Ny, f, f_temp);
f = periodic_left_right(Nx, Ny, f, f_temp);
f = periodic_top_bottom_2(Nx, Ny, f, f_temp);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% Compute streamlines %%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
u_x = u_x'; % Need to transpose these to get the
u_y = u_y'; % figure to turn out properly
% I have hard-wired the starting locations for the streamlines. That might
% not be a great way to do this, but I have found that it seems to work
% pretty well for most possible grain configurations.
y_start = 3:6:Ny;
vec_length = length(y_start);
x_start = ones(vec_length, 1);
x_start = [x_start', 6:6:Nx, 6:6:Nx];
new_length = length(6:6:Nx);
y_start = [y_start, ones(1,new_length), (Ny-1)*ones(1,new_length)];
XY = stream2(x', y', u_x, u_y, x_start, y_start);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% Graph streamlines %%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure(3);
orient(3, 'portrait');
set(3, 'PaperPosition', [1 2.5 6.5 6]);
h = streamline(XY);
set(h, 'Color', 'red');
set(h, 'LineWidth', 1.);
axis('equal');
axis([1 Nx 1 Ny]);
xlabel('x coordinate (lattice units)', 'FontSize', 16);
ylabel('y coordinate (lattice units)', 'FontSize', 16);
title('Streamlines of flow through a porous medium', 'FontSize', 16);
hold on;
contour(x, y, pore_wall);
hold off;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% Graph velocity vectors %%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure(4);
orient(4, 'portrait');
contour(x, y, pore_wall);
hold on;
h = quiver(x(4:4:Nx,4:4:Ny), y(4:4:Nx,4:4:Ny), (u_x(4:4:Nx,4:4:Ny))', (u_y(4:4:Nx,4:4:Ny))' );
set(h, 'Color', 'red');
xlabel('x coordinate (lattice units)', 'FontSize', 16);
ylabel('y coordinate (lattice units)', 'FontSize', 16);
title('Visualization of velocity vectors', 'FontSize', 16);
axis('equal');
axis([1 Nx 1 Ny]);
hold off;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% End the program %%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end