**6. Numerical Implementation and Test Problems**

The model is implemented similarly to the conventional LB *D*2*Q*9 model [3]. Firstly, we perform the collision step, then post-collision distribution functions are streamed at appropriate directions. It is well-known from the LB theory that the discretization of space–time affects the viscosity. The DV Boltzmann model discretized in a similar form as LB model reads as

$$f\_i(t + \delta t, \mathbf{r} + \mathbf{c}\_i \delta t) - f\_i(t, \mathbf{r}) = I\_i[f\_1, \dots, f\_N](t, \mathbf{r})\delta t,\tag{42}$$

by applying the Taylor expansion this equation can be rewritten as

$$\left(\frac{\partial}{\partial t} + \mathbf{c}\_{\dot{t}} \frac{\partial}{\partial \mathbf{r}}\right) f\_{\mathbf{i}}(t, \mathbf{r}) \delta t + \frac{1}{2} \left(\frac{\partial}{\partial t} + \mathbf{c}\_{\dot{t}} \frac{\partial}{\partial \mathbf{r}}\right)^2 f\_{\mathbf{i}}(t, \mathbf{r}) \delta t^2 + O(\delta t^3) = l\_{\mathbf{i}}[f\_1, \dots, f\_N](t, \mathbf{r}) \delta t, \mathbf{r}$$

then

$$I\_i\left(\frac{\partial}{\partial t} + \mathbf{c}\_i \frac{\partial}{\partial r}\right) f\_i(t, r) = I\_i[f\_1, \dots, f\_N](t, r) - \frac{1}{2} \left(\frac{\partial}{\partial t} + \mathbf{c}\_i \frac{\partial}{\partial r}\right)^2 f\_i(t, r) \delta t + O(\delta t^2), \quad i = 1, 2$$

therefore, we can conclude that the scheme (42) led to the hydrodynamic equations in which the contributions from <sup>−</sup><sup>1</sup> 2 *∂ <sup>∂</sup><sup>t</sup>* <sup>+</sup> *ci <sup>∂</sup> ∂r* 2 *fi*(*t*,*r*)*δt* + *O*(*δt* <sup>2</sup>) = <sup>−</sup>*δ<sup>t</sup>* 2 *d*2 *dt*<sup>2</sup> *fi* + *<sup>O</sup>*(*δ<sup>t</sup>* <sup>2</sup>) are present. The additional terms for the Navier–Stokes equations can be obtained with the application of the Chapman–Enskog expansion [3]. Note that the terms *O*(*δt* <sup>2</sup>) do not affect the Navier–Stokes equations, since they contain third order derivatives, which, in the Chapman–Enskog multiple-scale expansion, enter the equations for the moments at the Burnett level. For the Navier–Stokes equation, the additional viscosity terms result from <sup>−</sup>*δ<sup>t</sup>* 2 *d*2 *dt*<sup>2</sup> *f eq <sup>i</sup>* , its contribution to *<sup>∂</sup> <sup>∂</sup><sup>r</sup>* · *<sup>P</sup>*(1) is *<sup>δ</sup><sup>t</sup>* <sup>2</sup> <sup>∑</sup>*<sup>i</sup> cici* · *<sup>d</sup> dr d dt f eq <sup>i</sup>* , remembering that *<sup>d</sup> dt f eq <sup>i</sup>* can be expressed by (40), we eventually obtain *<sup>δ</sup><sup>t</sup>* <sup>6</sup> ( *<sup>∂</sup>*<sup>2</sup> *<sup>∂</sup>x*<sup>2</sup> <sup>+</sup> *<sup>∂</sup>*<sup>2</sup> *<sup>∂</sup>y*<sup>2</sup> )*u*.

Then, for the DV *D*2*Q*9 Boltzmann model in the form (42), the viscosity is given by

$$\nu = \frac{3}{4a} - \frac{\delta t}{6}.$$

In the simulations, the parameters are taken as follows

$$\alpha = \frac{3}{4\left(\nu + \frac{\delta t}{6}\right)'}, \quad \beta = 0.25a, \quad \gamma = 4a - 4\beta = 3a, \quad \lambda = 0,\tag{43}$$

i.e., we have six different collisions.

To validate the second-order convergence of the presented scheme, we estimate the simulation error defined as

$$error = \frac{\sqrt{\sum\_{z} (u\_m(z) - u\_{bmch}(z))^2}}{\sqrt{\sum\_{z} u\_{bmch}(z)^2}},\tag{44}$$

where *z* denotes the spatial variable, *um*, *ubench* are the modeled variable (velocity) and the benchmark solution, respectively. The convergence rate is evaluated by fitting the values of log(*error*) for the various log(*h*) = log( <sup>1</sup> *<sup>N</sup>* ) (*N* is the number of the lattice nodes, *h* is proportional to the lattice spacing) using the linear regression, the second-order convergence is achieved if the regression slope coefficient is close to 2.

Compared to LB *D*2*Q*9 model, the scheme (29)–(37) differs only in the collision term and the expression for the viscosity. This means that the computation time for (29)–(37) implemented in the form (42) is approximately the same as for LB *D*2*Q*9 model.
