Appropriate wall boundary condition for shear stress

Hello everyone,

first, I am very happy to have found such a nice forum dealing with LB.

My aim is to calculate the deviatoric stress tensor near the walls of my system. How shall I implement the no-slip boundary conditions so that the accuracy for the deviatoric shear stress is preferably large? A local approach would be nice.
My first benchmark using a simple bounceback scheme results in a relative error of up to 40 % compared to the analytic solution for a rectangular Poiseuille flow at the last fluid node layer near the wall. The velocity distribution however is reproduced very nicely.
Until now I have found nothing really useful in the literature.
Does anybody have a good idea?

Thanks a lot,

Hi Timm,

On a bounce-back node, you can’t compute the deviatoric stress from local information. This is because bounce-back nodes aren’t normal lattice Boltzmann nodes. The result from the Chapman-Enskog analysis which states that the off-equilibrium stress (pi = sum_i c_i c_i f_i - rho u u) is proportional to the deviatoric stress S is therefore no valid for these nodes.

I suggest that you have a look at the page on boundary conditions. This page distinguishes between “wet-node” boundaries and “bounce-back” nodes. Pi is proportional to S on wet-node boundaries, but not on bounce-back nodes.

This being said, there exist many boundary conditions in the wet-node category, among which you may pick out your favorite. As you prefer a purely local approach, I would suggest for example the Inamuro, Zou/He or Regularized boundary conditions. A comparison of these boundary conditions, as well as details of their implementation, can be found in this article.

Hi Jonas,

thanks for your very fast reply! I have read your paper comparing the five approaches to implement the velocity boundary conditions, and I like your own approach very much, since it has all the advantages I am interested in (3D, local).
Regarding this, I have one question. I have understood how to implement this approach for a planar wall, but how about edges, e. g.


where ‘X’ is obstacle and ‘O’ is fluid. How can I deal with the distribution functions f_i at the edge (bottom left)? And more complicated: What do I do if I have inclined walls with an arbitrary angle or even irregular shapes?

Thanks a lot again,

Unfortunately, there is no easy answer to either of your questions. To implement inclined boundaries, you need to use an “off-lattice boundary condition”, for which the boundary evolves across lattice nodes. This is required to resolve the shape of the boundary with sufficient accuracy, consistent with the accuracy of the LB scheme. The article by Guo e.a. suggests for example such a boundary condition within the framework of wet wall boundaries. You should however be aware that, from a software engineering point of view, the implementation of such boundary conditions requires substantially more efforts.

Corner nodes aren’t easy to deal with either. Although it is in principle possible to determine the stress tensor from incoming particle populations, and thus implement the regularized boundary condition, the usual trick to determine the density rho does not work on corner nodes. In OpenLB for example, corner nodes are for this reason always computed according to the Skordos algorithm, based on an interpolation scheme. Furthermore, the density is extrapolated from neighboring nodes.

Thanks Jonas,

I get your point. For the beginning I should just take the usual BB scheme, because my geometry is not trivial and contains not only inclined walls but also cylindrical elements. However, since I am at first interested in static problems (i. e. a fully developed flow in the geometry), are there any simplifications possible? In other words: Is there for a stationary flow a more accurate boundary scheme than BB which is simply to implement even at corner nodes?


the LB scheme is intrisically time dependent. Therefore there are no “simplification” which would allow better accuracy for stationary flows.

Actually, both full-way and half-way bounce-back are quite accurate. I suspect that there is a misunderstanding which makes you think they have lower accuracy than other boundary conditions. A common misconception is to think the wall is located on top of a lattice node, as in wet wall boundaries, while it really is located half-way between nodes. This is explained on the boundary condition page.

Of course, you can’t measure the strain-rate directly on a BB node, but then again, this node is not part of the fluid anyway. It is therefore not surprising that fluid quantities are ill defined. By computing the strain-rate, did you intend to compute the drag/lift force acting on an obstacle? In that case, with BB you should consider using the momentum-exchange method instead, as explained for example in this paper.

Thanks for your answers!
Actually I think I did not make my intention very clear. I am interested in the shear stress field itself. But what interests me most is the shear stress near the walls. I know that using half-way BB the obstacle surface is between two lattice nodes. And in this case I want to know the shear stress 0.5 lattice nodes away from the wall.
My preliminary results now are that the error of the shear stress (compared to the analytic solution) near the walls is much larger than say at the center of the flow. This in unfortunate. And that is the reason why I was looking for a modified BB scheme possibly decreasing the relative error of the shear stress near the wall.
However, now I will first perform a thorough error analysis and see if the my results are sufficient for my needs. Most probably I will use the usual BB scheme, because the implementation of another approach is too time-consuming for my geometry.
But this discussion has been very illustrative.


Dear Jonas,

in a former post in this thread you told me that the Skordos algorithm is used for interpolating the distribution functions on edges and that the density has to be extrapolated at those points.
Can you suggest one paper dealing with the latter issue? I would be very delighted.

Thank you,

I don’t know of any paper dealing with this issue. It’s just that the ususal approach (I think it was first published in the Zou/He article), which deduces the density value from the incoming particle populations and the imposed velocity, does not work on corner nodes (and on 3D edge nodes), where too little information is available. Most people circumvent this problem by restricting their attention to cases where they know the pressure on corner/3d-edge nodes from analytical considerations. Another approach, which is currently used in OpenLB, is to extrapolate the pressure value on these nodes from nearest neighbors. Yet another possibility would be not to compute the on-wall density directly, but use another closure criterion instead. This is for example done in the paper on mass-conserving boundaries, which Orestis suggested to you in another thread.

Dear Jonas and Orestis,

when testing the boundary conditions, have you ever taken a look at the shear stress? I know that you implement the boundary conditions in such a way that the overall behavior is second order accurate. But did you explicitly check if the shear profiles on the boundaries are correctly reproduced by the boundary conditions, or have you just investigated the velocity behavior?


I would say that the only problem in which the accuracy for velocity was explicitly investigated is the dipole-wall collision. In this case, the value compared with the literature is the enstrophy. As shown by the results, the shear is well reproduced by the wall, as the observed accuracy is superior to the one found e.g. with a finite difference method.

What do you actually mean by “shear profiles on the boundaries”? Are you referring to the stress tensor Pi? Or to velocity gradients computed by finite differences? Do you have unexpected results in your simulations?

Hi Jonas,

with “shear profiles” I mean the shear stress on the inlet and the outlet of the flow, using velocity boundary conditions.
I have compared three types of simulations (3D Poiseuille flow): your regularized approach, Skordos’ method with finite differences and a body force driven flow with periodic boundary conditions. The results for the velocity are all (more or less) as one would expect, but the accuracy of the deviatoric shear stress, i. e. \sum_i c_i c_i f_i^neq, strongly depends on the type of boundary condition.
My finding is that for larger Mach numbers (> 0.04) your regularized method produces slightly inaccurate results as compared to Skordos’ method - but only on the boundary itself. This is of course not problematic, since the velocity is behaving nicely and only few people are interested in the correct shear behavior. So my question is, if you have observed something like this in your simulations.
I have to admit that I use fullway BB on the side walls - in contrast to your benchmarks. I wonder if the interplay BB at side walls <=> velocity boundary conditions at inlet / outlet influences the flow behavior significantly. But I do not have enough time for testing every possible setup.


The thing is, I don’t think that one can expect the linear relationship between shear stress (Pi^neq = \sum_i c_i c_i f_i^neq) and strain rate (S = 1/2 rho(d_alpha u_beta + d_beta u_alpha)) to hold with second-order accuracy. To recover the Navier-Stokes equations in a Chapman-Enskog analysis, the relation between Pi^neq and S is required to hold at first order only. And it seems to me that the asymptotic analysis of Michael Junk e.a., which is independent of Chapman-Enskog, draws similar conclusions.

In practice, it may be that better-than-first-order accuracy is observed, but I don’t think that you can expect this to hold in general. My interpretation is that second-order terms with respect to the grid spacing are “spoiled” by non-hydrodynamic terms in the multi-scale expansion. To measure velocity gradients in a simulation, I would therefore guess that you are better off using a finite difference scheme rather than computing Pi^neq. This doesn’t mean that a boundary condition using Pi^neq is less accurate than a boundary condition using finite differences, because the boundary condition only needs to be self-consistent and does not necessarily depend on an explicit relationship between Pi^neq and S. All I am saying is that if you want to explicitly measure velocity gradients, finite differences may perform a better job than local lattice variables.

This is a very interesting post, and I tend to agree.
Getting velocity gradients from FD is a sound method with high accuracy and returns comparable if not better results in all my simulations. Yet, I am interested in the basic understanding of the shear stress and its role in the LBM. And I have found out that this is anything but a straightforward task.