I’m new on this forum and also in LB methods and I have a question about accuracy of BGK and MRT methods in a simple case : the lid driven cavity. I simulated the lid-driven cavity for Reynods = 100 with the two codes and my problem is the following~: when I plot the evolution of the error ((my_solution - benchmark_solution)/benchmark_solution) over spatial discretisation, I have a 1st order accuracy (ie. the error is divided by 2 when the the spatial step is divided by 2). Maybe I’m wrong but I read that LB schemes are 2nd order in space (and time) so I don’t understand why is this case, there is a fisrt order accuracy (the BC are quite the same than Zou & He and they have a 2nd ordrer accuracy in a Poiseuille case).
Although the LBE is supposed to be second order accurate in space and time, the restriction on the Mach number usual means it is more like first order in time (see the work of people like Junk et al and Holdych et al). Also (and this may only be semi-relevant), the “usual” LBE is a first order in time approximation to the discrete Boltzmann equation (that is, the Boltzmann equation with a finite set of 9 particle velocities). This is because the LBE can be obtained from an Euler integration of the discrete Boltzmann equation. For a second order implementation of the DBE you can discretise in space and time using the trapezium rule (or by Crank-Nicholson).
Then there is the issue of Zou and He boundary conditions. I think it gives an incorrect mass flux, so the scheme will not be second order accurate for pressure and could infect your velocity field. This will not show up in Poiseuille flow because there is no vertical component of velocity.
If you are using the streamfunction for comparison (which is the usual benchamark test), you have to use a suitably accurate scheme numerical approximation to calculate this (and be careful near the walls).
By benchmark solution you mean the solution by Ghia e.a., which provides the velocity in a few chosen points, right?
The problem here is that the LB method (like most CFD methods) converges uniformly, but not point-wise, with second-order accuracy. To measure the rate of convergence of your method, you should monitor an integral quantity, such as, the L2-norm of the difference between numerical and analytical solution over the full domain.
If no analytical solution is available (as in the lid-driven cavity), your best bet is to compute a numerical reference solution at a very high resolution: your Matlab codes will most likely not be able to produce such a result. Plus, it is rarely a good idea to measure the rate of convergence towards your own solution, as your code might converge rapidly, but to a wrong solution. It is actually not an easy task to measure the convergence rate of a code in a problem without analytic solution…
In the case of the lid-driven cavity, and if you are serious about wanting to measure the convergence rate of your Matlab codes, I would suggest that you generate reference data (in the whole domain) by means of a spectral CFD method. These methods are simple to apply to the square geometry of the cavity, and very precise, as they converge with spectral accuracy. You can for example take the OpenSpeculoos code, which implements a spectral element method, and which provides the lid-driven cavity as one of the example applications. You need however to dig a bit into the code, because in the end you will need to re-interpolate the results onto the equally-spaced LB grid.
Last but not least, remember that lattice Boltzmann converges in a diffusive limit. The time-step must evolve as the square of the space step. To put it differently, when you double the space resolution, you must also divide the velocity, as measured in lattice units, by two, to get a second-order rate of convergence.
Yes, that’s a very good point and I’d agree that an L2 norm comparison with a high-order spectral method (over the entire flow) is the best way to judge. It may quickly be becoming a controversial paper but I think this is what Luo et al (2011) have done. Ignoring the ELBE stuff, the paper shows order of accuracy and convergence of the LBE compared to spectral methods (if I remember correctly).
Regarding the order (or accuracy, for the moment) of the LBE in time; if we take the diffusive scaling epsilon^2=dx^2=dt, then we do indeed see the method is a second order approx. to the incompressible Navier-Stokes equations, with respect to epsilon and dx. But this is then only first order in dt, hence why we must half the velocity if we double the resolution. Yes, doing so gives us better convergence properties but requiring order (dx^2) timesteps is quite inefficient, especially for high Re numbers (and the error terms have strange dependencies on dt/tau, which can alter the convergence properties). A second order discretisation of the discrete velocity model allows you to use larger timesteps.
I certainly don’t know this for sure but I have a feeling you’ll have problems simulating the cavity flow with Re more than, say, 2000 on a reasonably sized mesh using Zou and He boundaries (with BGK, anyway). If I am wrong (and I often am!) then please correct me but I think (half-way) bounce back, for all its flaws, is more stable than Zou and He for this flow.
I fact, I compared my results with a spectrale method solution (Botella and Peret 1998, Benchmark spectral results on the lid-driven cavity flow). I modified the orginal matlab code and introduced half-way BB and on the moving wall, instead of Zou e.a. bc i used a bc found in the Latt et e.a. courses (http://moodle.epfl.ch/course/view.php?id=4501, lesson 7). At first, my first objective was to understand full-way a half-way bb schemes.
For Re = 100, it shows good results and a second order convergence for quantities such as max u at mid-width. (I get a 1st order with full BB).
Calculations for re = 1000 are in progress.
For both cases, I observe oscillations on the pressure field near the 4 corners (with large grid resolution, with a finner grid, it disapears). Maybe i have to smooth the imposed velocity profile in order to reduce the generated noise.
But, I agree with you, a more rigourous way is to compare the full field obtained with a spectral method with the LB one, maybe i will try this later.