Simulation Verification
The discussion on this page is based off of our work in our Dimensionless Verification IPython Notebook.
Verification via Pipe Flow
After our Lattice Boltzmann code initially ran and looked qualitatively correct (i.e. no gigantic shock waves at slow flow speeds), we wanted to make sure that it was quantitatively correct. We therefore needed a test case that could be used to simulate and compare with analytical predictions. We decided to use 2-dimensional flow in a pipe because the analytical steady-state profile is well known. Simplifying the Navier Stokes equations (as discussed in What is the Lattice Boltzmann Technique?) for a 1-dimensional flow in a pipe, the equation governing horizontal velocity in the pipe is
Solving it using the no-slip (u=0) boundary conditions on the walls, we find
where y is the distance away from the bottom of the pipe and D is the diameter of the pipe. The shape is quadratic, and, unsurprisingly, maximum in the center of the pipe.
We implemented a Lattice Boltzmann simulator for a pipe flow and imposed a no-slip boundary condition along the walls of the pipe. To implement zero velocity on the walls of the pipe, we used a first-order accurate bounce-back technique; see the boundary condition section of this link for details. To implement a pressure difference across the pipe, we used a more subtle technique discussed in our Boundary Conditions page. We essentially specified an appropriate density difference across the pipe. This can be obtained from the pressure in the Lattice Boltzmann model through the equation of state P=ρcs2. Because of the algebraic manipulations to solve for the boundary conditions, this took a lot of time to get right. The resulting Poiseuille flow from our numerical can be seen below .
Dimensionless Verification & Resolution
Note that the horizontal velocity is in terms of "dimensionless units." What does this mean? From the What is the Lattice Boltzmann Technique? page, we know that the Navier-Stokes equation is
Since we solved for u above, we can calculate umax; the characteristic time scale is thus
Following the procedure outlined by the Palabos guide on non-dimensionalizing Lattice Boltzmann simulations and non-dimensionalizing every parameter with units, we find that the non-dimensional Navier-Stokes equation is given as
where u-tilde is a dimensionless velocity and P-tilde is a dimensionless pressure. For example u is related to u-tilde.
The only parameter that changes when there are different inputs is non-dimensional parameter Re, the Reynold's Number
Since the non-dimensional Navier-Stokes PDE is equivalent to the dimensional one after rescaling variables, two simulations are ultimately the same if they have the same Reynold's number and geometry. This suggests a convenient way to run simulations: take the physical parameters of your simulation, determine T and L and the corresponding Re, then specify a resolution parameter N which will run the simulation with more grid-points but will ensure that the Reynold's number of the simulation remains the same. The accuracy of the simulation should increase with N for two models with the same Reynold's number.
Note that this "dimensionless resolution" idea is discussed in greater detail in the Palabos guide on non-dimensionalizing Lattice Boltzmann simulations.
Dimensionless Simulation Convergence
We implemented our simulations in a dimensionless manner as discussed above and slowly increased the resolution of our simulation. As expected, as N increased, the steady-state profile of the pipe approached a constant shape. However, based on our input parameters, the velocity profile in the pipe did not exactly match theory (see the "Theory" curve below). If we multiply the theoretical curve by a constant value, however, it always has the correct parabolic shape (see the "Scaled Theory" curve below). The factor that the theoretical curve must be multiplied by is never greater than a factor of 2, however, and seems to depend on the input viscosity.
We guess that the difference between theory and simulation is caused by slightly incorrect dissipation at the walls because we used a boundary condition that was only first-order accurate, not second-order accurate in terms of enforcing zero velocity. We have plans to fix this problem in the future. Since the factor multiplying the theoretical curve was not too large (i.e. less than two) and the shape of the simulation curve was correct, we guess that this problem did not drastically alter the behavior of our simulations.