Multiphase flow in porous media for digital rocks

Hi all,

I recently created a repository with my code for simulating multiphase flow through porous material (steady and unsteady approaches) Using the Shan Chen model. Feel free to visit it and let me know what you think.


Comments and suggestions are appreciated



I have recently been studying multiphase flow in porous media,I read your code on guthub,You wrote it very well, but I have a question, has palabos incorporated the Gads parameters, the fluid and solid forces.I didn’t find the code for the Gads parameter call in palabos. Do I need to write the fluid and solid forces by myself?
Thank you for your reply

Hi @liuqi,

Thank you for your kind words. I have been working very hard to make a clean interface for people to build on top of it. For porous media applications, the G_ads parameter is very important, since it controls the wettability of the system, by setting a specific contact angle at the solid surface. Different contact angles can be set by printing different numbers at the solid surface. I currently support up to 5 different contact angles, by I can add more if necessary.

To do this, you will need to assign the numbers 1,3,5,6 through your geometry (as done here, and then add the contact angles that you want in those voxels here

An deeper explanation of the process is provided in this readme:

Thanks for your interest

@liuqi if you want a simple example to start experimenting with different contact angles, change this value:

The ranges of it should be from -0.4 to 0.4 (non-wetting to wetting, drainage->imbibition) according to Huang et al., 2007.

I’ve figured out how to simulate forces by setting a virtual density on the boundary.Thank you very much, I learned a lot from your code, and I have a question about the two-phase boundary condition, which boundary is better for export in palabos when the import is constant speed?Or can you provide me with some code for boundary conditions that I can learn about?Thank you for your guidance and reply

Hello, I’ve found your code base, and I think it’s good. It fits my needs.But I didn’t run it successfully. Can you help me?Thanks

Sorry @Taihengliu I did not get a notification for this.
Have you tried running the matlab code that creates the geometry?

Hi, what am I going to do at the end of the unsteady_relperm_spherepack operation? I just looked at the generated data using ParaView.
Thank you for your reply.

Hallo guys, i am new to palabos and i am still trying to figure out some codes. i have encountered a hurdle as i was trying to run an old code on which runs well on the older version of palabos but failing to build using cmake on the recent palabos version. i am getting the following response on the terminal. please help.

Hello lkadzungura,

welcome tho the Palabos Forum!

Try to substitute the deprecated auto_ptr with the unique_ptr (or the shared_ptr), I think it could solve the issue.

I hope this helps,


Hie Mars

you have helped a lot, thanx. it worked but now i have to look for the executable file that it makes. that is not a problem though. thanx a lot.

Hi,why do I use the new version of the package and compute slower than the old version?

I have run your code and found that when the density ratio and viscosity ratio are 1. omega1 and omega2 are both 1, the flow process conforms well to the requirements of Laplace’s Law, but when the viscosity ratio increases, omega1=1.6, omega2 =1, the pressure gradient is 0, and the flow also occurs. Why?Can you help me to explain the problem?Looking forward to your reply!

That’s a common set-back of the Shan-Chen method. Let’s wait until the color-gradient method it’s released.

How is the surface tension of the two phase interface determined in the program? thanks

Hi Javier
How is the surface tension value 0.15 obtained?
Thank you for your reply.

Howdy @Taihengliu,

Feel free to check Chapter 5 of my thesis:


1 Like

I use your program, if I don’t change Gc, the interfacial tension doesn’t change either, it’s still 0.15.
Thank you for your reply.

Can your program be modified to make fluid 2 drive fluid 1?If so, how should I modify it?
Thank you for your reply!

Howdy @Taihengliu,

Fluid 1 and 2 are interchangeable labels. I believe that you would like to change the contact angle and that can be achieved by changing the G_ads of the appropriate surface.