porous media geometries in openlb

nice work on the open source code, web home and forum!

i’m interested in running LB on 3D CT data to determine permeability/ porosity relations for sea ice (brine inclusions within pure ice matrix). input data will be segmented voxels: brine or ice.

looking through the lessons in the intro pdf, it seems that it will be quite a task to specify each corner, edge and face of the very complex geometry - and this would ned to be done for multiple data sets of ~ megavoxels. i can imagine a script to do this and generate the required code. does one exist already, and what are your thoughts on porous media applications of openlb?


Porous media – that’s a difficult topic for which you’ll not find a lot of support in the current OpenLB version. The velocity boundary conditions in OpenLB are definitely not appropriate for modeling porous media, for reasons of efficiency, accuracy and code engineering. Your best bet is to use bounce-back for your solid nodes. That works actually pretty well if you are happy with no-slip walls (porous media people often are) and if it’s sufficient for you to represent the boundaries through a stair-case approximation (this seems to be your case, given that you were ready to use the velocity boundaries to represent your geometry).

Full-way bounce-back is delightfully easy to use, as it is independent of orientation: whatever the orientation of your boundary, you can always use the same type of dynamics. See the OpenLB example on the 2D flow past a cylinder to get an illustration (the same approach works in 3D).

Actually, Orestis and Wim have already modeled porous media in OpenLB using this approach. They may have some hints or ideas they could share with us here.

as Jonas said you should use bounceback for wall nodes in your porous media (and noDynamics on “in-wall” nodes). The thing we did for our simulations was to assign a bounceBackDynamics on the wall nodes inside the domain and a pressure BC on the inlet and outlet domain.


Jonas, Orestis - thanks a lot for your suggestions. I have also been in touch with Wim. I am encouraged that my planned work seems to follow closely his poster presentation at AGU07 - some of it already done, and for sea ice not volcanic material.

the obvious thing to need to be able to do is import the CT files (tifs in my case) before assigning dynamics to wet/ wall / in-wall nodes. is this capability already within OpenLB, or did you just write this directly into your code?

thanks again.

there is no feature to read tiff files in OpenLB for the moment.
In fact what we do with wim is to first convert the slides of figures (our are bmp files) into a format that is directly readable by c++ (pure text) and assign a number to each node (0=bulk, 1=wall, 2=no dynamics). Then when we launch the simulation we have a simple loop that define the dynamics dependeing on this number.

ok, your approach sounds like what i’d figured it would be, i’ll go for something similar then.

I forgot to tell you that we used matlab to convert our images to text.

I am also interested in modeling porous media, but the problem that I am facing is that my boundary conditions can change during the simulation. I want to couple OpenLB to another program that calculates the rate of growth/decay of the pores. Can boundary conditions be changed dynamically within OpenLB? Could not make out from the documentation whether this is possible.


Oh, good point. This should be made more clear in the documentation then. Yes, you can change the value of the velocity resp. the density on boundaries dynamically (you can also change the type of boundary condition dynamically, but I think that’s not what you are talking about). Once a lattice site has been defined to be a velocity boundary node by use of the method defineVelocityBoundary(), the value of the velocity can be modified at any time by a command of the type lattice.get(iX,iY).defineVelocity(vel) . If it’s a pressure boundary, you may modify the value of the density in the same way.

Thanks. I will look into those functions. I really like the way you all have set this up for open-source access. This is truly a tremendous resource for the community. BTW, in an earlier post you had mentioned some other work already done on porous media using OpenLB. Would you know if this has been published?

Hi gersappe,
in fact the work on porous media using OpenLB has not been published yet. We plan to do so when we finish all this. There is not even a preprint for the moment.


Hi Orestis and Wim

I am using the LBM for porous media (252252252) simulation. I have already started my simulation and my gif images are looking good, even though the simulation is taking a lot of time. I am wondering if you know of any method/formula that can be used to directly calculate the permeability of a medium from the particle distribution functions.

I found a paper that used the body force instead of pressure difference but the permeability formula that was given there did not make sense. I could send to you if you wish.



Hello everybody,

I have been recently trying to use OpenLB for flow simulation in porous media (geology, in my case). I have started with cylinder2d example as the template. I have modified it as follows:

// Definition of the obstacle (bounce-back nodes)
for (int iX=0; iX<nx; ++iX) {
for (int iY=0; iY<ny; ++iY) {
T u[2] = {poiseuilleVelocity(iY),0.};
T rho = (T)1 + poiseuillePressure(iX) *

     //my code
     if( (rand()%10)+1<2) // set rand number from 1 to 10 to the lattice$
     // instead of this
      //  if ( (iX-obst_x)*(iX-obst_x) +   (iY-obst_y)*(iY-obst_y) <= obst_r*obst_r )


in order to do first step and create random structure.

but I got:

av energy=nan; av rho=nan
convert: Invalid pixel ./tmp/p000100.ppm' @ coders/pnm.c/ConstrainPixel/141. convert: Invalid pixel./tmp/u000100.ppm’ @ coders/pnm.c/ConstrainPixel/141.

and so on … and gifs look strange. Only in the first column i can see fluid.

Could anyone give me some sugestions ? Thanks in advance.
Best regards,



See the code below…When I use the portion that has been commented out, I get the same error as you. But When I use the code you suggest, everything worked out fine.


        if ( (iX-obst_x)*(iX-obst_x) +
             (iZ-obst_z)*(iZ-obst_z) <= obst_r*obst_r )
            lattice.defineDynamics( iX,iX,iZ,iZ,
                    &instances::getBounceBack<T,DESCRIPTOR>() );
        else if (iZ < nz/4 - 1)
              if (rand()%10+1<8)
                lattice.defineDynamics( iX,iX,iZ,iZ,
                    &instances::getBounceBack<T,DESCRIPTOR>() );
	/*if (iX%2==0 && iZ%2==1)
            lattice.defineDynamics( iX,iX,iZ,iZ,
                    &instances::getBounceBack<T,DESCRIPTOR>() );
            */ //End the very ordered bounce back rule