# Lattice Boltzmann methods vs Navier stokes/ other eulerian methods for *water* simulation

Note, there is already a question here, however the answers don't answer the original question, let alone specific considerations when dealing with nearly in-compressible fluids (water). Another similar question here, but does not speak to water, and really only talks about performance considerations. I'm not too concerned about performance right now.

I've been trying to figure out what method to use for CFD with water. I've tried simulating with LBM, and no matter what I do, I get massive instability with infinities and NaN's everywhere after some time, it seems inevitable (using float32's, as non scientific Nvidia Gpus can have 1/32 double precision throughput). The thing is LBM is extremely easy to implement, and Navier Stokes methods are not.

So what I'm trying to figure out is LBM even the right tool for the job? I see people talk about LBM in terms of compressible fluid, but for water and the scale I'm dealing with this would appear moot, it isn't necessary to deal with. So if I'm dealing with water, and I don't want to deal with compressibility, should I just stick with Navier Stokes? Or is there some way to get LBM to behave for water?

1. Note, I'm primarily concerned about the force of water on macroscopic objects, and what effect the water has on that single object's movement.
2. I'll be using a GPU, and I will be implementing everything (and have already done so for LBM on the GPU), my simulations are specific enough that other peoples frameworks will not work (though I also don't need the precision of simulation that they need)
• Hi there. When you say you want to simulate water, what kind of properties do you want to include in your model. You already named compressibility. If you want to study the interaction of water flow, then most likely the effects of compressibility will not be very important (except if you are in some high velocity setting). I would start out with dismissing compressibility! – MPIchael Oct 4 '19 at 9:30
• @MPIchael Yes compressibility is not important ( I mention this in my question ), in fact the only thing I'm actually concerned about is the effect of water on the movement of a moving rigid body object within water. The object will move which will also displace water, which will also have an effect of the movement of said object. These kinds of effects are the ones I'm hoping to model, ie like fish movement through water. – whn Oct 4 '19 at 12:09
• I use LBM for almost three years in my PhD research and before that for many years I used FEM and FVM for CFD. In my opinion: LBM is only suitable for creepy flows in structured grid geometries. Otherwise, even in incompressible regime, if you have intermediate Reynolds number ~ 100-200 you will be in trouble... Unless you have a really powerful cluster with many core hours queue reserved specially for you to do simulation. I strongly recommend to use FEM or FVM for any kind of CFD simulations (incompressible, compressible, etc.). – Alone Programmer Oct 4 '19 at 13:00
• After three years days and nights using LBM for more than 1000 really large-scale CFD simulations (my number of lattice points are 100 Millions to 200 Millions and number of timesteps ~ 3000000, so imagine the nightmare that I have right now just for simulating a flow in intermediate laminar regime of Re ~ 500 to 1000): Don't waste your time with LBM... Yes it's really easier to implement and parallaelize it (even it's not true really, I can refer you to my code which has ~ 2 millions C++ lines parallelized with MPI, which only god knows what's going on there...) but it doesn't worth it... – Alone Programmer Oct 4 '19 at 13:03
• And about GPU: LBM and GPU is just waste of time... Unless your model is small enough to be able to fit it in a single GPU without any need to have hybrid MPI+GPU architecture. MPI+GPU that we recently implemented in our code with hope of getting better performance: It really sucks... why? Because of really expensive handshaking of communications between GPUs. In fact, recently we realized that communication between GPUs through MPI is about 1000 times slower... So, don't waste your time and go ahead with many powerful FEM or FVM based CFD solvers like OpenFOAM and Nek5000. Believe me bro... – Alone Programmer Oct 4 '19 at 13:05

First, I'm not sure why you emphasize on water? I mean I understand that you are looking for a CFD scheme that works for your special case, but you need to know that water fluid is not a special fluid at all. Water is a incompressible fluid and you can simulate its movement by using incompressible Navier-Stokes equation as long as your Mach number is not high, typically smaller than 0.3. So, by knowing your Mach number and the Reynolds number for your specific case, you can pick incompressible Navier-Stokes equation to simulate your water flow.

You have two options here:

1. Solve Incompressible Navier-Stokes equation: $$\nabla \cdot \mathbf{u} = 0$$

$$\rho\frac{\partial \mathbf{u}}{\partial t} + \rho \mathbf{u} \cdot \nabla \mathbf{u} = \nabla \cdot \tau + \mathbf{f}$$

For Newtonian fluid: $$\tau = -P\mathrm{I} + \mu(\nabla \otimes \mathbf{u} + (\nabla \otimes \mathbf{u})^{T})$$

Finally: $$\rho \frac{\partial \mathbf{u}}{\partial t} + \rho \mathbf{u} \cdot \nabla \mathbf{u} = -\nabla P + \mu \nabla^{2} \mathbf{u} + \mathbf{f}$$

You can solve it by using FEM or FVM methods, which is successfully implemented and verified in OpenFOAM and Nek5000 packages. I recommend Nek5000 which even has industrial scale accuracy.

1. Use lattice Boltzmann method (LBM): $$f_{i}(\mathbf{r}+\Delta t \mathbf{c}_{i},t+\Delta t) = f_{i}(\mathbf{r},t) - \frac{f_{i}^{eq}-f_{i}}{\tau}$$

For ensuring conditional stability you need to have relaxation time $$\tau > 0.5$$ and $$\mathrm{Mach} < 0.02$$. Note that LBM $$\mathrm{Mach}$$ is not equal to real Mach number. These are the conditions which may work for your case or for countless cases will not work at all.

If I was in your shoe and I had the knowledge about LBM that I have right now three years ago, I would definitely choose the first option. Any sane person will go for the first option I believe. The most famous argument that we LBM people tell other people including traditional CFD people is: LBM is easy to implement cause it just contains two operations of collision and streaming and in collision step everything could be calculated locally and in streaming step you need to just copy everything in certain directions in comparison to traditional CFD techniques that you have nonlinear momentum acceleration term ($$\mathbf{u} \cdot \nabla \mathbf{u}$$) which is a pain to discretize and also for the pressure you need to solve time-consuming Poisson pressure equation: $$\nabla^{2} P = - (\nabla \otimes \mathbf{u})^{T} : (\nabla \otimes \mathbf{u})$$. Also nowadays, we argue that meshing step is more time-consuming in comparison to LBM, cause we call LBM a mesh free method. All of them are true for creepy flows ($$Re << 1$$) in structured grids like cube, which of course it just limits us to toy model and is not suitable for real world engineering applications. In fact, one iteration of LBM is much much faster than one iteration of FEM or FVM but the problem is that in order to satisfy those stability conditions for LBM you need so many iterations which is at least three order of magnitude higher than conventional FEM or FVM flow solvers. If your adviser is a LBM expert and you are forced to use LBM and of course if you have lots of computational resources, go for it, tough luck! But, if you don't, I recommend to pick a famous open source or commerical CFD solver and more focus on engineering aspects instead of thinking why your simulation become unstable. Already, you see that it is really difficult to get a stable LBM simulation, so just write down the cons pros of these two approaches and take the one that is more inclined to your knowledge, resources, etc.

Update:

Let me put some calculation here to show that indeed LBM needs a huge number of meshes and timesteps to work properly at mid range Re flow regime. Let's say we have a water pipe ($$\nu = 10^{-6}$$ $$\frac{\mathrm{m}^{2}}{\mathrm{s}}$$) with 1.0 mm diameter and 1 $$\frac{\mathrm{m}}{\mathrm{s}}$$ velocity in it, which is not really a big deal to solve with FEM or FVM. The Reynolds number is: $$Re = \frac{D u}{\nu} = \frac{10^{-3} \times 1}{10^{6}} = 1000$$. For LBM stability, a rule of thumb is: $$Mach < 0.02$$ and $$\tau = 0.54$$ for Mach number and relaxation time respectively. So:

$$\tau - \frac{1}{2} = \frac{3 \nu \Delta t}{\Delta x^{2}}$$

$$Mach = \frac{\sqrt{3} u \Delta t}{\Delta x}$$

If $$Mach < 0.02$$:

$$\frac{\Delta t}{\Delta x} < \frac{0.02}{\sqrt{3} u} = 0.012$$

But:

$$0.04 = \frac{3 \times 10^{-6} \Delta t}{\Delta x^{2}}$$

So:

$$\frac{\Delta t}{\Delta x^{2}} = 13333.33$$

Finally:

$$\frac{\Delta t}{\Delta x} = 13333.33 \Delta x$$

Then:

$$13333.33 \Delta x < 0.012$$

$$\Delta x < 10^{-6}$$

I pick the upper bound as $$\Delta x = 10^{-6}$$ m. So:

$$\Delta t = 1.3 \times 10^{-8}$$

First, look at the insane timestep that we calculated. It means if you want to simulate water flow in a pipe for 1 second in reality it needs almost 100 millions timesteps! Wow!

Let's say your 1 mm diameter pipe has 10 mm or 1 cm length, which is not that uncommon. So, the total volume of your pipe is:

$$V = \frac{\pi}{4}D^{2} L = \frac{\pi}{4} (10^{-3})^{2} \times 10 \times 10^{-3} = 8 \times 10^{-9}$$

Remember LBM works with voxels and it means your mesh is just a tiny cube and has same size everywhere in your domain (ridiculous right?! forget about adaptive mesh refinement in LBM...). So, each voxel volume is:

$$V_{voxel} = \Delta x^{3} = 10^{-18}$$

Finally, you would find the number of voxels that you need to fill this 1 mm diameter and 10 mm length pipe as:

$$N = \frac{V}{V_{voxel}} = \frac{8 \times 10^{-9}}{10^{-18}} = 8 \times 10^{9}$$

Another insane finding... You need 8 billions voxels to do this simulation and still you can't be 100% sure that your simulation will remain stable or not... Come on... FEM or FVM people do this simulation as a piece of cake less than a couple of seconds on their crappy old laptops by using ANSYS Fluent or any other open source code if you want to argue that ANSYS is commercial. Want to say something about MRT? That's garbage... It's just good for publishing nonsense articles in simple 2D rectangular geometries to show that they were able to simulate lid driven cavity problem with higher Reynolds number. Just look at the CFD market and count besides XFlow and Exa who produce LBM software for serious industrial usage? For your information Exa is acquired by Dassault Systemes and I'm pretty sure they just bought Exa to close it.

• This in combination with everything else I've read, tells me that LBM is a no go for the kind of simulation I'm doing. I think I'm going to use PIC/FLIP methods, as I need to be fast (not necessarily 100% accurate), though all my simulations are submerged. Also those resources you've linked to will be a good resource for me to verify the accuracy (or lack of) of my solutions! – whn Oct 4 '19 at 20:03

I disagree with the answer given by Alone Programmer. His reasoning is not objective and seems to be based on a very primitive lattice-Boltzmann model with BGK collision operator. LBM has significantly matured over the years and is in my opinion a very attractive numerical solver. Nonetheless it depends on the application if it really fits your needs.

I have been programming LBM multi-threaded OpenMP algorithms now for 2 years and I am still convinced that the method or a modification of it will be a very attractive tool for transient incompressible flows and in the future also for compressible simulations. The following thoughts consider positive and negative aspects of LBM:

• LBM is based on the kinetic theory of gases. In the kinetic theory of gases one creates fictive dilute gas models of discrete particles or particle clusters interacting with each other and observes their dynamics. Boltzmann introduced a probability distribution over phase space, an extended concept of the traditional density that is only strictly valid for dilute gases. Surprisingly the continuous Boltzmann equation even with simplified collision models such as the BGK approximation preserves the complete set of compressible Navier-Stokes equations (Chapman-Enskog expansion). This means the theoretical framework can be used to make predictions about compressible and incompressible flows. One can make use of the law of similarity, meaning it does not even matter if you simulate a liquid or a fluid: As long as the assumptions about the material parameters such as density and the dimensionless parameters of interest in the differential equations such as the Reynolds and Mach-number are identical, the behaviour is similar. Traditional LBM methods apply a discretisation of the space with distinct single-speed (they only connect neighbouring cells) velocity subsets- so called lattices. These have a certain number of speeds and thus a number of "symmetries". These symmetries are connected to tensors and these again to conserved quantities and conservation equations. The common discretisations have too little speeds in order to preserve the full Navier-Stokes equations and you are left with the conservation equations for incompressible flow with an error term that rises with the model Mach number. The model is therefore used to preserve the incompressible conservation equations asymptotically. You can simulate substances in the limit of a low-Mach number incompressible approximation wherever that is an incompressible fluid like water or the flow around a car. That being said there is room for improvement: There has been continuous effort in particular from the research group of Prof. Ilya Karlin in Zürich to make modifications to the basic idea and apply it to (highly) compressible flows. The more basic nature of LBM - namely the distribution functions- as supposed to the direct application of the continuum hypothesis and the Navier-Stokes equations - might allow for a simulation of highly compressible fluids with shock waves where the continuum hypothesis clearly breaks down. If such an algorithm will really be faster than finite-volume-based methods has still to be seen but I have the confidence it can be.

• Incompressible LBM struggle with medium Reynolds numbers already when used with the standard formulation with BGK operator and basic boundaries such as the 2D Zou/He boundaries. But there are significantly more stable formulations such as MRT, TRT, cascaded and entropic collision operators and furthermore several modifications of the Smagorinsky turbulence model and wall functions that can be implemented efficiently. These formulations also require for more sophisticated boundary conditions and strategies like block-structures mesh refinement are available. This means while the basic algorithm is simple incorporating a stable model is significantly more difficult, doing so effectively even more so. Nonetheless the limitations discussed by Alone Programmer are not there. You can simulate highly turbulent flows even on desktop setups. The limitations arise from the nature of the models (DNS, LES, DES vs RANS) and not from LBM as a CFD solver!

• LBM traditionally works on a regular structured grid. This means meshing is easy but boundaries are stair-cased if one uses the traditional bounce-back algorithm. Implementing more precise boundaries like interpolated bounce-back violates though exact mass-conservation.

• In LBM you have limited options for turbulent inlet conditions. You have to use experimentally determined velocity profiles, trip wires to trigger turbulence, pressure-periodic boundary conditions or a periodic box.

• The LBM is comparably simple and predictable. Thus the algorithm is very fast and its parallel scalability is excellent, unrivaled by finite-volume based methods if you apply it to transient problems. You can find implementations in the literature on CUDA GPUs and Xeon Phi coprocessors that update up to 1 billion lattice nodes a second. On a 12-core processor I have archieved 150 million lattice node updates with a D3Q19 lattice. Do not compare apples to oranges: lattice-Boltzmann is inherently transient and requires a lot of times steps you can control the accuracy and temporal resolution of the algorithm by setting the macroscopic velocity and therefore a particular model Mach number. The model Mach number should be lower than 0.05 across the majority of the flow field to obtain usable results. You can obtain faster results with finite-volume method if you calculate steady-state processes and you can use RANS models but let's face it: Eddy-viscosity RANS models such as $$k$$-$$\epsilon$$ or $$k$$-$$\omega$$, have severe drawbacks. They are inherently isotropic and thus struggle with anisotropy such as in swirling flows. $$k$$-$$\epsilon$$ further overestimates the turbulent intensity rendering it pretty much useless for most tasks. On the other hand you have very expensive RSM models that have 6 distinct transport equations for the turbulent stresses. They handle anisotropy way better but they have stability issues and so only few modern publications on this topic can be found. URANS formulations altogether offer very limited insight into transient turbulent processes. DNS is almost impossible for technical applications and LES is still restricted to large-scale clusters. DES + LES is computationally very expensive but delivers good results. Now where does LBM go? It is a solver like FVM and FDM. The traditional formulation is basically DNS but there exist RANS and as mentioned LES formulations. Without a doubt LBM is way faster than FVM when performing DNS or LES but the question is would a RANS model suffice? There is though a reason truck aerodynamics where separation prevails are dominated by LBM: If one needs a transient resolution of turbulent structures like vorticity LBM is clearly to be preferred! Just look at what XFlow can do and look at how much of that traditional commercial software such as Ansys Fluent simply can't because of the different modelling paradigm.

• Nonetheless LBM is currently actively researched, it lacks a wide acceptance and readily available code. There are only few optimised options available such as the open-source projects Palabos, OpenLB and the commercial softwares XFlow and ProLB that allow somebody to use LBM without programming it by yourself from scratch. LBM codes are often used in research, custom made, tailored to a particular platform by skilled personnel.

So in my opinion the reason why LBM is becoming increasingly more popular is because of its capabilities to resolve transient flows accurately and effectively. Writing a well performing LBM algorithm from scratch requires though time and I would not advise anybody to do so for their PhD thesis. Most supervisors will expect you to model physics and likely won't give you the time you would have to put into a solver itself. If you already have a well-performing framework then it is an excellent choice for transient flows in particular in complex geometries. But you should ask yourself: Does a cheaper model suffice? If you do not need the transient information of LBM simulations then every time step is a burden! If the smeared time-averaged results of cheap RANS models are sufficient turn to OpenFOAM.

• Unfortunately, I should disagree. You pretty much described LBM and its various applications/implementations, etc., which are not necessarily interesting to OP or even answering this question. The main point in your answer is that you provided inherent transient nature of LBM as some sort of a unique feature, which I think you just exaggerating it... There are a couple of well defined, validated, verified, transient CFD frameworks based on FEM or FVM, so LBM does not add anything here or does not have any unique feature. – Alone Programmer Oct 30 '19 at 1:50
• You talked about more complicated collision operators such as MRT, TRT, etc. which mathematicians can understand what's going on with them, but in my experience as someone use all of these collision operators not just simple LBGK, in my opinion these are just nonsense for engineering applications. If you are a mathematicians or a physicist that study a fundamental problem of fluid dynamics, probably LBM with its all LBGK, MRT, TRT, filtering, BFL boundary conditions, etc might be a good option for you, but for engineering application, never LBM... – Alone Programmer Oct 30 '19 at 1:53
• I agree about kinetic theory, Chapman-Enskog expansion, etc that might give a better understanding of fluid dynamics fundamental aspects or probably from a non-equilibrium thermodynamics point of view, but I can't understand why an engineer should use LBM when there are accurate, fast, scalable, multiphysics, parallel implementation of CFD based on FEM or FVM such as OpenFOAM or Nek5000. This debate will never get stopped cause some people just love LBM... But my opinion is based on my three years of experience working days and nights with LBM. Nothing is in LBM right now that I did not try... – Alone Programmer Oct 30 '19 at 1:58
• So again: if your flow is creepy ($Re << 1$) and steady state and you don't want to spend money to buy license of Fluent or CFX, the best option is to pick an open source LBM implementation such as OpenLB or Palabos. Other than that using LBM is just torturing yourself... I have a friend that her research is about creepy flows in porous media and she's pretty OK with using LBM cause she could fit her simulations even in her desktop computer without the need for any cluster, but for me to deal with 100 to 200 millions lattice points with 3000000 timesteps, it's just an endless nightmare... – Alone Programmer Oct 30 '19 at 2:13
• @Alone Programmer I have listed advantages and disadvantages of the method and left it to the questioner to decide. The transient nature and the small required time steps are unique features of LBM which help it for transient applications. You are comparing a framework to a method. LBM as a method is excellent for transient applications, just the easily accessible frameworks are limited. And for certain applications there are no alternatives. E.g. find me a solver that is able to mesh complex geometries such as realistic porous media. Nobody will use FVM for that for sure. – 2b-t Oct 30 '19 at 8:20