I have run a series of 3D Poiseuille flow simulations with the following parameters:
tau = 0.92 (fixed), Re = 1 (fixed) and additionally Ma proportional to lattice spacing
I use full-way BB at the side walls and velocity boundary conditions at inlet and outlet (both fully developed flows).
My aim is to check the second order accuracy of the LBM and its convergence to the NS equations, when the Mach number and the lattice spacing are reduced in the same way.

The overall L2-error I calculate is defined like this:
E = sqrt[(sum ||vel_an - vel_sim||^2) / (sum ||vel_an||^2)]
The sum takes into account the entire lattice including inlet / outlet (except the side walls).

In total I have 13 data points ranging between lattice sizes 6x12x48 and 42x84x336.
The result is that the velocity error has a slope of -2.33 instead of 2. What does it mean? What could be the reason?

I have a rectangular channel with cross section ratio 2:1 and length to shortest width ratio 8:1, e. g. 336x42x84. The flow enters at x = 1 (inlet) and exits at x = 336 (outlet). The laminar steady state solution is known analytically, but it contains an infinite series, see for example here.

However, I think I have made a mistake with the velocity boundary conditions. I am about to recheck the results and compare it to other boundary conditions.

I have tested the 3D Poiseuille flow for a body force driven flow and with velocity boundary conditions (both inlet and outlet) using the Skordos approach (BC4 in your paper). As stated above, I have computed the overall L2 error of the numerical result, compared to the analytic solution.
For a body force driven flow with periodic boundary conditions I get a -2.01-slope for the velocity error, which is a very nice value. For the VBC profile, however, there is still a -2.31-slope, although I am now pretty sure that the implementation is correct. As I can see from your publication (Fig. 8), not all error curves for the 3D Poiseuille flow have a -2-slope. Can you tell me which slopes the curves have? I am not sure, if the deviation (-2.31 instead of -2) is due to the model or if I have done anything wrong.
Or do you have other suggestions? What might be the reason for the deviation?

I am not sure, but -2.3 sounds like a pretty steep curve to me. Is it possible that you are running at too high Mach number, and that the slope is asymptotically -2 as you lower the Mach number?

Anyway, I am sending you the data points of the paper in a PM. You can also reproduce the points by running the program referenced in the article.

Thank you Jonas for the data. I will have a look at it in short.
But first: The simulation parameters are the same for the periodic and the velocity, starting at Ma = 0.02 and ending at Ma = 0.0028 obeying \Delta x \propto Ma. So I do not think that the Mach number is too large.
The strange thing is that for the periodic BC I really have 2.01, so the effect should be due to the boundary conditions themselves or my lack of implementing them correctly. And even if I only take into account the data points with the smallest Mach numbers (say four or five data points), there is nearly no change in the power law index. Thus, the power law index is constant for my data set.

@Jonas
Now I have checked the slope in your paper for the finite difference method: it is very close to -2.
Is it possible that the interplay of BB on the side walls and VBC on the inlet / outlet could cause the deviation?
You did not use BB on the side walls, but VBC with zero velocity, is it true?

That’s true we use wet-node boundaries everywhere.

There’s an optimistic and a pessimistic way of interpreting your results. The optimistic one states that your simulation is good, because the error decreases at a better-than-second-order rate as you increase the resolution. The pessimistic one states that your simulation is bad, because the error increases at a worse-than-second-order rate as you decrease the resolution. There are two ways of knowing which one is right. The first is to run the code in our article and compare the accuracy with yours, in absolute numbers (not the asymptotic rate). If yours is better, this indicates that the optimistic interpretation is (probably) right, and if its worse, you’ll tend towards the pessimistic one. The other approach is to keep increasing the lattice resolution. In the optimistic interpretation, your slope should keep on being around -2.3, whereas in the pessimistic interpretation, it seems to me that it should approach the value -2 at some point.

An explanation for the optimistic case would be that LB is better than second-order accurate on a Poiseuille flow with BB boundaries. That might be possible, who knows: even 3D Poiseuille has many symmetries which are not present in generic flows, and it wouldn’t be surprising that LB does a better job than expected in this case. An explanation for the pessimistic case would be that some additional errors appear during usage of a too coarse grid, which you don’t see in the asymptotic analysis. Or that your code has bug. Such as, misinterpretation of the position of the wall with BB, or whatever other possibility one may choose among the many pitfalls of LB.

Maybe I have another explanation (not really an explanation in fact). As you noticed and posted on this forum the accuracy of the LB simulations are quite sensitive to the value of relaxation time you use. I guess that the choice of tau=0.92 is not innocent, is it? Maybe you could try with another value (farther away from 1) and see what happens…

Thanks Jonas and Orestis.
I had already the same idea with the “good” and the “bad” interpretation of the results.
So, Orestis, you think that the convergence may depend on tau. May I ask which values of tau you have used for your benchmark tests in the straight velocity paper?
There is another idea: I have used an aspect ratio of 2:1 for the channel cross-section. You have used 1:1. Right now I test the channel with cross-section ratio 1:1. This may also be a reason. What is the length of the 3D channel in your benchmarks?

All values, including the resolution and the aspect ratio of the channel are described in the paper. As usual, the viscosity is computed by using the definition of the Reynolds number, Re = U*N/nu .

I see, so there are some minor differences.
My channel is much longer than yours, making the simulations running longer. I should change that.
My error definition is slightly different, I should also check that.
So, for now, I have plenty of work to do. Maybe I find an explanation for the behavior I have observed.
Thanks a lot again!