I have written an LBM code for the case of 2D Taylor-Green-vortex decaying problem using the following parameters:
delta t = delta x, Re = 10, U=0.01
Domain: [-1 1] x [-1 1]
Mesh size: 21 x 21, 41 x 41, 81 x 81, 161 x 161.
The error was measured using the L2-norm of the difference between the numerical solution and the analytical solution, at time = 1.
When I chose U = 0.01 or below, the numerical solution exhibits 2nd order accuracy. However, this is not the case for U>0.01 (although the solutions for those cases remained close to the analytical solution).
I am wondering if this is normal, or could it be due to some mistakes in my program?
The LBE, as you usually see it, is formally second order accurate in space and first order in time. We get (spatial) second order accurate predictions for the incompressible Navier-Stokes equations despite using a first order discretisation of the (discrete) Boltzmann equation (i.e. the Boltzmann equation with a discrete set of particle velocities) if we use a diffusive scaling , i.e. dx^2~dt. See the work of Junk for a detailed analysis.
Increasing the velocity (U) you use in your simulations probably means you are increasing the time step if you use the same number of points, so it is perhaps to be expected. Second order accuracy under the acoustic scaling (dt~dx) is more likely to be obtainable if you use an LBE that is obtained by, for example, a Crank-Nicolson discretisation (i.e. a second order in time integration of the discrete Boltzmann equation), but note you are solving the compressible Navier-Stokes equations, although you still need to keep the Mach number small.