Lattice Boltzmann simulation of 2D flow in a random porous medium

Background for the project

The calculation of Fluid flow inside a porous medium is one of the important problems in the field of environmental fluid dynamics, mostly in the context of aquifer contamination and movement of solutes in the sub-soil layers. The important quantity to be calculated is the velocity vector field of the fluid. Most methods in computational fluid dynamics employ the well known Navier-Stokes equations of momentum conservation in some manner of space and time discretization of the continuous PDE's to find the velocity field, but in a porous medium the boundary conditions are usually too complex to solve by such an approach.

The Lattice Boltzmann method on the other hand, is based not on a continuum approach but a discretized particle approach, and therefore is very suitable for use in the realm of porous medium, since every point in space is either occupied by a "fluid particle" or a "solid particle", without any need to describe a continuous boundary. It also lends itself naturally to a parallel computing structure, which is another important feature in regards to the high level of detail that is necessary to describe the complex geometry of the medium.

The Lattice Boltzmann method

Press Here for a very clear explanation about 2D Lattice Boltzmann, including governing equations, boundary conditions and algorithm implementation.

Notice that in this project, the boundary condition schemes that are used are the ones described in the article - the bounce-back B.C. on solid boundaries and the Zou-He velocity and pressure B.C. in the inlet and outlet boundaries.

If you are interested in the development of the Zou-He boundary conditions, you can read the original article published by the authors, on which I based the development of the outlet velocity calculation.

The program

The program is basically divided into 3 parts -

1. Pre processor - Generates a random porous medium.

2. Processor - Solves the velocity field in the generated porous medium.

3. Post processor - Creates a GIF animation of the time evolution of the velocity field.

The random porous medium generator is written in matlab, since it has a better random number generator than fortran. The porous medium is represented by a binary matrix in which '1' represents a solid particle and '0' represents a fluid particle. To create the porous medium, the matrix starts as solid (all ones) and the program tunnels through it by creating zeros at several fixed starting points and then creating more zeros around them in a random order, until a wanted porosity of the matrix is achieved.

The Lattice Boltzmann solver is coded in fortran for higher performance speed. The algorithm is structured after the following scheme:

1. Read porous domain from generator.

2. Read user input for run time and pressure gradient.

3. Initialize distribution function, density and velocity field.

4. Calculate new distribution function by streaming step, bounce-back and Zou-He pressure B.C.

5. Calculate new density and velocity field.

6. Calculate new distribution function by collision step.

7. Write velocity field to output file.

8. Repeat 4-7 until time specified by user.

9. Write domain to output file.

10. Write printing commands to output file.

Graphic results are displayed using matlab. The processor creates a matlab file with the results, including the commands needed to turn them into a GIF animation. The velocity field vectors are represented as blue arrows, created by the matlab function 'quiver'. The size of the arrows represents their relative magnitude, in relation to the maximum velocity vector.

Files and commands

There are only 2 files needed to run the program - LBM.f90 and generator.m.

The user needs to have access to matlab (with graphics) and to a fortran 90 compiler. If you are working from home and don't have matlab, you can connect to aluf server using 'ssh -X -l username' and run matlab on the terminal.

The procedure for running the program is as follows:

1. Download LBM.f90 and generator.m to the same folder. You can copy them to aluf using:

scp /local folder/filename

2. Run generator.m in matlab - you will be asked to choose length and width of the porous domain (in mm) and porosity. Recommended values appear in parentheses. Remember that the larger the domain, the longer is the computation time. When choosing porosity (which is always between 0-1) it is preferable to choose above 0.5, or you might get a domain with no connection between the inlet (left side) and the outlet (right side). The generator prints the domain on the screen when it is done, and gives you the option of creating a new domain if you are not pleased with the one generated. If you choose to create a new domain, you can also change the porosity. If not, it will be saved for the processor to use.

3. Now you have a file named domain.txt in your folder - it contains the data of your generated porous medium.

4. Compile LBM.f90 in the same folder - if you are in phelafel or aluf, you can use the line:

gfortran LBM.f90 -o LBM.ex

Which will create an executable file named LBM.ex

5. Run the executable file and choose time and pressure gradient, again the recommended values are in parentheses. Choose the file name for the GIF animation (up to 8 characters). The program will calculate results and tell you when it is done. The output to the folder is a matlab file named print.m

6. Now all you have to do is run print.m in matlab to get a GIF animation of the flow.

IMPORTANT - when the file is running in matlab, do not go to other windows, since in order to make the GIF animation, matlab is taking snapshots of each graphic figure - if you move to another window it will take a snapshot of that window instead.

Some results

Flow at a pressure gradient of 0.0001 atm, domain size 20 mm X 20 mm, porosity 0.6, time 100 seconds

flow animation

Flow at a pressure gradient of 10 atm, same domain, time 100 seconds

flow animation

Flow at a pressure gradient of 10 atm, domain size 50 mm X 50 mm, porosity 0.65, time 200 seconds

flow animation

It is also possible to create other domains using matlab, for example to simulate flow around an infinite cylinder

flow animation

Parallel code structure

Creating a parallel code for the Lattice Boltzmann method is not very complicated, as suggested earlier. The solver is practically the same; the only thing that has to be done is to decompose the domain into subdomains before sending each subdomain to a different processor, together with its matching distribution function. It is important that every subdomain has a set of boundary conditions - therefore there must be overlapping between the subdomains (each one is calculated based on points surrounding it that were calculated in the previous iteration). This system of overlapping decomposition is best understood graphically:


This picture represents an 18 X 18 matrix, divided into 9 subdomains, sized 6 X 6, each one surrounded by a red line. The blue squares are the original domain's boundaries (they are outside of the domain), the yellow squares are the overlapping boundaries of the subdomains. When calculating, each processor receives an 8 X 8 subdomain, and returns the answer for the internal 6 X 6 subdomain. The main processor then connects the subdomains back to reconstruct the original domain, and distributes them again to the other processors.

The algorithm for the parallel code would be as follows:

1. Create random porous medium with the generator.

2. Initialize the distribution function.

3. Decompose the domain and the distribution function.

4. Send each subdomain with matching distribution function to a sub processor.

5. Calculate new distribution function and velocity field at each subdomain.

6. Send back to main processor and reconstruct distribution function and velocity field for full domain.

7. Print velocity field to output file.

8. Repeat 3-7 until time specified by user.

You can try the decomposition for yourself - Use the generator to create a 20 X 20 domain (with any porosity), then download the file overlap.f90 to the folder containing domain.txt, compile it with gfortran and run it. You will see how the program breaks down the full domain into 9 overlapping subdomains of size 8 X 8.

Contact me at for any questions about the project.