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

Question

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)

Answer 1

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.

2. 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.

Answer 2

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. Irrespectively of the numerical method chosen you have to think about whether you really want to maintain your own code base or if it is not better to use and extend existing software tools.

---

I have been programming LBM multi-threaded OpenMP algorithms now for 2 years (update 2022: now over 4 years off and on) 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 (update 2022: this has proven true, a lot of LBM start-ups and companies have emerged in the meanwhile and XFlow has proven to be very competitive, e.g. see their simulation of the Bugatti Chiron). 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. (Update 2022: With increasing processing power of CPUs and GPUs modern implementations are able to update between 300 million cells on a modern CPU and 3.5 billion cells on a modern accelerator card, see e.g. the benchmark of the University of Geneva) 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.3$ across the majority of the flow field to obtain usable results (see e.g. Aslan et al.). 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. *(update 2022: This has gotten better, in particular Palabos and OpenLB have worked on their GPU support! XFlow and SimScale have gained traction.).

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.

Answer 3

Just happened to have a look at this old post and happened to notice that the given calculation in the accepted answer is scientifically completely (!) wrong such that I had to outline it in this post:

• The given computational burden is wrong by four orders of magnitude: A calculation which actually takes less than 2 hours and a half on a modern GPU would take close to three years and close to 40 years on a modern CPU with his assumptions! Cringe!
• The answer does not differentiate between a model, a numerical method and an implementation giving a skewed view and rendering the given answer virtually worthless.

Even though the author claims to have extended experience with the method, I have the impression that his/her frustration (which is visible throughout the answer resulting in a highly subjective answer) is caused by a lack of profound knowledge of the theoretical foundations of the method as well as its efficient implementation. I will outline this in the post below.

---

• The post lacks a clear differentiation between numerical methods and implementations: Finite element method (FEM), finite volume method (FVM) and the lattice-Boltzmann method (LBM) are methods and should be clearly separated from their implementations which are software tools/libraries such as OpenFOAM, Nek5000 or Ansys Fluent. There are also open-source (e.g. Palabos, OpenLB) and commercial solvers (e.g. Simulia XFlow by Dassault Systemes or Pacefish used by SimScale) based on the lattice-Boltzmann method which offer also advanced features like grid-refinement, large-eddy and RANS turbulence models. It is sufficient to have a look what modern LBM codes with large-eddy turbulence models can do, such as this aerodynamics simulation of the flow around a Bugatti by XFlow, to see that the comparison is completely lacking context: You can't compare your own buggy, in-house code to a commercially developed and fully tested professional framework in order to decide which method is more suitable than the other. This is in particular true for LBM where a main source of instability are low-order boundary conditions that fail to model higher-order contributions that still have an effect on $f^{(neq)}$ and therefore on the evolution of the Boltzmann equation.

• Finite element methods are not popular in computational fluid dynamics due to conservativity and are generally significantly slower. Neither Nek5000 (finite volume, finite differences, spectral elements) nor OpenFOAM (finite volume) are based on finite element.

• The answer mixes methods (see point above) and levels of modelling and turbulence models (such as Direct Numerical Simulation (DNS), large-eddy turbulence models and (unsteady) Reynolds-averaging (RANS)). Directly solving the Navier-Stokes equations (DNS) for a transient flow is quite challenging in particular with traditional numerical methods such as finite element methods. A commercial solution based on FVM will only be able to solve a pipe flow at Reynolds $1000$ reasonably fast when assuming steady flow (so not transient, which is reasonably for the given Reynolds number but would not be reasonable for higher Reynolds numbers as the flow will transition to turbulent flow for $Re_D \geq 2300$) and deploying an unstructured grid. Even then lattice-Boltzmann will be about as fast when calculating transiently, in particular if taking into consideration grid generation. Said scenario is actually the worst case for a comparison and not representative: LBM can't accelerate towards steady-state as it is inherently transient - even though the solution is steady you have to perform a transient LBM simulation.

• The local Mach number required for stability given by $0.02$ is a magnitude smaller than you would actually need (see e.g. Aslan et al. - "Investigation of the Lattice Boltzmann SRT and MRT Stability for Lid Driven Cavity Flow") making the required number of time steps a magnitude (and more) larger.

• The Mach number in lattice-Boltzmann is defined as

  $$Ma := \frac{u_{lb}^{max}}{c_s}$$

  where $c_s$ is the lattice speed of sound which for the most common single-speed lattices is given by $c_s = \frac{1}{\sqrt{3}} \frac{\Delta x}{\Delta t}$. But $\Delta x$ and $\Delta t$ are lattice units, the dimensionless system in LBM which chooses $\Delta x = 1$ and $\Delta t = 1$ (You are never ever calculating any $\Delta x$ or $\Delta t$ with the standard incompressible LBM! Just have a look at Latt or Krueger et al.! You are working on the basis of a dimensionless system with convenient dimensionless units that is related to the dimensioned system through the law of similarity. What is then done in the derivation above just hurts if you have seen any introductory video to the lattice-Boltzmann equation.). This means your highest possible over-velocity in the fluid domain is

  $$|u_{lb}^{max}| = c_s \, Ma \approx \frac{1}{3} \frac{\Delta x}{\Delta t} \, 0.3 = 0.1$$

  If we assume that there are no significant contractions inside the domain we can assume that the for a pipe flow the characteristic velocity can be chosen to the maximal admittable velocity $U_{lb} = |u_{lb}^{max}|$.

  Even keeping the very conservative stability limiy $\tau \approx 0.54$ as a stability limit (which is only required for simple collision operators like BGK without any turbulence models, while other collision operators are closer to $0.51$ and turbulences models allow arbitrary close relaxation times to $0.5$ that are only limited by the floating-point accuracy so something like $10^{-15}$ for the double data-type) this gives with the formula for the model viscosity

  $$ \nu = \left( \frac{\tau}{\Delta t} - \frac{1}{2} \right) c_s^2 \Delta t $$

  a limiting model viscosity $\nu_{lb} = 0.01\overline{3}$ which means that with you have to choose your characteristic length according to the Reynolds number $Re_D = \frac{U_{lb} \, D_{lb}}{\nu_{lb}}$ resulting in $D_{lb} = \frac{Re \, \nu_{lb}}{U_{lb}}$ which in our case gives $$D_{lb} = \frac{Re \, \nu_{lb}}{U_{lb}} = \frac{1000 \cdot 0.01\overline{3}}{0.1} = 133.\overline{3}$$ for $D_{lb}$ which is the diameter of the pipe $D$ in lattice units. Assuming a rectangular computational domain that is $C$ times longer than its diameter this gives a number of cells $$N = NX \cdot NY \cdot NZ \approx C L_{lb}^3$$ which for a mentioned $C = 10$ gives $2.37 \cdot 10^7$ nodes, so 23 million nodes which are ehhhhm, a bit less than the mentioned 8 billion. And as mentioned these stability limits are not relevant for more advanced collision operators, turbulence models and in any case resolution requirements might be reduced with grid-refinement approaches.

• Same goes for the time step: For unit conversion we can introduce the conversion factor between the dimensionless $lb$ system and the dimensioned system $d$ in SI-units: $c_L = \frac{L_d}{L_{lb}}$ for length, $c_U = \frac{U_d}{U_{lb}}$ for velocity resulting in a conversion factor for time $c_T = \frac{c_L}{c_U}$ and therefore in $T_{lb} = \frac{T_d}{c_T} = \frac{c_U \, T_d}{c_L}$. For the given numbers $c_L = \frac{10^{-3}}{133} \approx 7.5 \cdot 10^{-6}$, $c_U = \frac{1}{0.1} = 10$ and therefore $T_{lb} \approx \frac{10 \cdot 1}{7.5 \cdot 10^{-6}} = 1.\overline{3} \cdot 10^6$. So roughly 1.3 million time steps are required for a stable simulation of 1 second for the given system, sounds like a lot but it is not that much: A well-programmed three-dimensional lattice-Boltzmann implementations can update around 200-300 million nodes within one second on a single CPU (see e.g. Wittmann et al. or Bauer et al.), implementations on the old Xeon Phis can reach around 1 billion node updates per second (see Robertsén et al.) and on a graphics/accelerator card well over 3.5 billion nodes per second (see Latt et al.). This means even for what seems to be an unreasonable large number of time steps the method will be very fast!

• Therefore the simulation of one second real-time with a primitive BGK collision operator for the given Reynolds number of $1000$ resulting in a grid with $23$ million computational nodes and $10$ million time steps on a single GPU will remain stable and complete in approximately $\frac{23 \cdot 10^6 \cdot 1.3 \cdot 10^6}{3500 \cdot 10^6} = 8542.\overline{857142}$ seconds, so around 2 hours and 20 minutes which is reasonable for a transient simulation.

• If you are though implement a more stable collision operator such as multiple-relaxation times (MRT, see d'Humieres) or cumulant (see Geier), or if you add a turbulence model such as the Smagorinsky large-eddy model (see e.g. Sagaut) your stability requires significantly less computational nodes (something like $67$ or $33$ cells or even lower - likely anything less than $67$ is too little for resolving the physics with a uniform grid, you would have to calculate your dimensionless wall distance $y+$ and see, probably one would reduce the simulation domain to $5 \, D$ in favour of having a larger $D_{lb}$ anyways). Assuming a simulation domain of $67 \cdot 67 \cdot 667$ we would have $8$ times less cells, with just under $3$ million nodes, $c_L$ would drop to the half $1.49 \cdot 10^{-5}$ and so would the required time steps $666666.\overline{6}$. Our computational burden decreases by a factor of $16$ resulting in $571$ seconds or less than 10 minutes simulation for 1 second of real-time (on a single A100 GPU or less than 2 hours on a powerful desktop). This means the LBM simulation on a GPU will voxelise the geometry (sub-seconds) and complete the simulation before you have even created your computational mesh for a finite volume simulation.