How to start with porous media simulation ?


I am going to perform a simulation of 1-phase flow in porous media (permeability assessment) . The geometry will be imported form uCT images. I plan to use d3q19 model
Can I have some basic questions?

  1. Which code will be best suitable (OpenLB or Palabos ?)
  2. Which example I should start from ?
  3. How to introduce pressure gradient in order to calculate permeability according to Darcy’s law.

I know that in case of porous media the bounce back bc will be best suitable. Has anyone any experience in this kind of simulation ? I will be very grateful for any suggestions.

Best regards,

The simplest thing to do is to convert your geometry file into a simple ASCII format: 1 for fluid nodes and 0 for solid nodes. The file examples/codesByTopic/io/loadGeometry.cpp shows how to read the geometry and instantiate the bounce-back nodes. To create a pressure gradient and test Darcy’s law, use standard pressure boundary conditions on the inlet and the outlet, as explained in the Palabos user’s guide. Note: if the solid nodes overlap with the inlet and/or the inlet, you should first create the inlet/outlet pressure boundaries, and thereafter instantiate the bounce-back nodes.

All of this is really straightforward to do and works well for various applications in geophysics. We’ve tested it multiple times. Just ask further questions if there’s an issue.

1 Like

Thank you jlatt for your help! I have already written simple software to threshold and convert ct images to binary data. I will try use loadGeometry as a template, but I am interested in 3D case (D3Q19) - I hope that modifications will be quite simple. I am also confused about lattice size. Where in the “loadGeometry” lattice size is defined? Sorry for my basic question, but I am not familiar with Palabos yet.

Best wishes

Yes, translating from 2D to 3D should be straightforward (the syntax is identical). Attention: the geometry must be in an ASCII file, not a binary file. This is a simple text file which contains the characters “1” and “0”.

The dimensions of the lattice are defined at the point where the lattice is created:

MultiBlockLattice2D<T, DESCRIPTOR> lattice (
              parameters.getNx(), parameters.getNy(),
              new BGKdynamics<T,DESCRIPTOR>(parameters.getOmega()) );

Here, the object “parameters” is used to compute the dimensions of the lattice and the value of the relaxation parameter, to automatically convert between physical and lattice units. You don’t need to use this object if you don’t want. For example, to get a 100x100 lattice and a relaxation parameter omega=1 for the BGK model, use the following command:

MultiBlockLattice2D<T, DESCRIPTOR> lattice (
              100, 100, new BGKdynamics<T,DESCRIPTOR>(1.) );


Thanks jlatt for help !

Can I have a few further question ?

  1. At first I am a little bit confused about size of geometry.dat file. It consists of 9881 characters (without spaces). The size of lattice is 40640*1 = 9600. Am I right that something is wrong with this file ?

  2. What kind of fluid model would be suitable for this kind of simulation? Would Single-relaxation-time BGK be fine ?

  3. What about periodicity ? Should I make my lattice periodic in any direction? (I reckon that in this case it should be periodic in every direction. Shouldn’t it be ?)

  4. Finally, what about initial condition? Initializing all cells to an equilibrium distribution with the constant velocity value will be fine?

I would be grateful for any help.


  1. The size of the lattice is (406+1) * (401+1) = 9881. The +1 comes from the fact that the simulation assumes on-lattice boundary conditions. Basically, if you want a lattice which spans over 2 lattice cells, you need three lattice nodes (node 0, node 1, and node 2), because the walls (or the inlet/outlet) are exactly on node 0 and on node 2, and there are 2 cell-widths between them.

  2. I’ve seen plenty of people doing this kind of simulation with single-relaxation-time BGK, and it certainly works well. I’d suggest that you start with BGK and see what you get. In case you find a particular reason to switch to MRT, you can still do this later on (it’s just about changing one line in the program).

  3. It doesn’t seem to me that you have periodicity between inlet and outlet. After all, they don’t have the same pressure, right? Technically speaking, it doesn’t hurt if you turn on periodicity in the lattice, though, because the periodicity is explicitly interrupted when you define the pressure inlet and outlet. Thus, you get the same result, whether periodicity is on or off in the inlet-outlet direction. In the other two directions, it depends on your physics whether you want periodicity or not. If you just want to verify Darcy’s law, I’d say that periodic boundaries are an excellent choice.

  4. Yes, this sounds fine (a zero velocity is probably just as good and spares you the trouble of thinking about the value of the constant velocity). You shouldn’t run into problems, even if the initial condition is bad. A bad initial condition in this case just means that the simulation takes more time to converge to a stationary state. If you want to speed up the convergence, you could probably care about varying the pressure in the initial condition in such a way that it decreases linearly between the inlet and the outlet value.


Thanks again!

Now I am a few step closer to the solution. I have written “jpg -> 0 and 1 converter” I can load 3d geometry, and produce VTK files. Everything seems to be fine, but the problem is in inlet and outlet boundary conditions. I have started with simple tunnel on 10x10x15 lattice. The part of my code where I define bc is fallowing


boundaryCondition->setVelocityConditionOnBlockBoundaries (lattice, Box3D(14, 14, 0,9,0,9 ) ); 
setBoundaryDensity(lattice, lattice.getBoundingBox(), 1.2 );

boundaryCondition->setVelocityConditionOnBlockBoundaries (lattice, Box3D(0, 0, 0,9,0,9 ) );
setBoundaryDensity(lattice, lattice.getBoundingBox(), 0.8 );

Do I properly define boundary conditions in order to achieve pressure gradient ? I understand that in BGK model pressure and density are connected according to ideal gas law and there is no additional method to define pressure?

I obtained some results ( simulation runs) but physical interpretation is rather strange or difficult.
I will be very grateful for any further help.

Best regards


To get a pressure condition, you must use “setPressureCondition…” instead of “setVelocityCondition”.

Right, the ideal gas law states that p = 1/3 rho. However, when you compute p, you need to convert everything to lattice units. I am pretty sure that 0.8 vs 1.2 is a too large pressure difference. You can read the following document about unit conversion:

Or search the forum for “unit conversion”.


Sorry for my another silly question. I have read the article and searched out the forum…

does pressure difference convention go like this:

DeltaP_{LB}=1/3+DeltaP_{Phy}*dt_{phy}^2/dx_{phy}^2 ?

I just want to be sure that everything is correct.

Can anyone give me some hints how to use the “converter” in Palabos code, please?

Hi, if it’s about pressure difference, then the 1/3 goes away. For the rest, the formula looks OK. As an example, you can use a density of rho=1 at the inlet and rho = 1 - 3*(DeltaP_{Phy}*dt_{phy}^2/dx_{phy}^2) at the outlet, if DeltaP is positive.


Thanks a lot!

I am still confused.
In my simulation:

  1. I have to define Re but I don’t know u0 which I here understand as fluid speed (correctly?) - this is exactly what I need to calculate for Darcy law. Maybe there is different way to calculate Re before simulation ?

  2. A problem with the converter. I try to use converter i.e

LBunits<T> converter(uMax, Re, 50., 1., 1., 1.); 

but I have got:
loadGeometry.cpp:61: error: ‘LBunits’ was not declared in this scope
during compliation. What is wrong ?
Anyway I still can’t use the converter as long as I don’t know Re number.

Best regards

Hi, the LBunits class is now called IncomprFlowParam.

Now I see :slight_smile: thanks a lot ! Any hint to Re number I should use ?


You can circumvent the problem by working entirely in lattice units; you don’t need to know the Reynolds number to verify Darcy’s law. Start by doing the unit conversion from physics to lattice units so you have a rough estimate of an appropriate value of deltaP in lattice units (and use Darcy’s law to know the velocity, which you put into the definition of the Reynolds number. Then simply vary the pressure gradient, and verify Darcy’s law right away, working in lattice units.


Am not sure if I understand. My goal is to find permeability of some sandstones, not only general verification of Darcy law, so physical units are necessary, I reckon.

Dear All,

  1. How can I compute dt_phy ? I have read that dt should be proportional to dx^2 in dimensionless units. What about dt_phy ?

  2. I can’t also work out how to compute filtration velocity in Palabos (discharge per unit area, with units of length per time). I need it to Darcy law. Any hints ?



I am not sure if it is a proper approach but in order to compute the filtration velocity I am going to compute norm of average velocity component perpendicular to inlet and outlet. Then filtration velocity should be equal to this norm times porosity ? Am I right ? Has anyone any experience in porous media simulation in Palabos and could help me ? I would appreciate and kind of help. Thanks in advance.


Maybe I try to ask more specific question.
My physical porous medium:

lx = 800 [um]
ly = 800 [um]
lz = 800 [um]

average single pore size = about 50 [um]
porosity = about 25%
fluid = water => nu_phy = 10^-6 [m^2/s]

I have ASCII file 200x200x200 characters which is obtained from uCT measurement.
I have imposed pressure boundary condition on inlet and outlet. I added 25 free fluid surface on inlet and outlet.

So, I suppouse that deltaX_phy=deltaY_phy=deltaZ_phy = 800[um]/200 = 4[um]

and deltaX_lu= … 1/200 ? right?

  1. How should I define variables in IncomprFlowParam? So far I use
    lx=ly=lz =200
    Re =1
    Umax =0.02
    but shouldn’t it be rather like:

lx=ly=lz = 1 … ?

I am still confused about it.

  1. How can I estimate Re ? I want to simulate flow of water so I know nu_phy, but how about characteristic velocity.

  2. I want to change value of pressure difference in order to obtain permeability from Darcy law, but how can I compute discharge of fluid and convert it into physical unit.

I would appreciate any kind of help.
Sorry for maybe quite silly quesitons.

Best regards

Kuba Wrote:


Sorry for my another silly question. I have read
the article and searched out the forum…

does pressure difference convention go like this:


I just want to be sure that everything is

Can anyone give me some hints how to use the
“converter” in Palabos code, please?

On which paper did you see this formula:DeltaP_{LB}=1/3+DeltaP_{Phy}*dt_{phy}^2/dx_{phy}^2?
i simulate simple component multphase flow in porous media use Shan-Chen model.i have a problem about the pressure.when i use the BGK model the pressure inlet is many times greater than the Shan-Chen model when other parameter are same.did anyone have the same question?

I´m working on something quite similar, although I did not move to 3D so far.

Right now I´m also busy with implementing physical units, so i can´t help you there (yet).

For the boundary conditions and the applied pressure gradient, i simply apply an external force on every lattice cell, the code and my answers to finally get it working in 2D can be found in,2111
I also use periodic boundaries, and the applied force in x-direction represents a higher pressure on the right and a lower on the left border.

Maybe this helps a bit,