Recalibration of LBM Populations for Construction of Grid Refinement with no Interpolation

1. Introduction

Lattice Boltzmann Method (LBM) is one of the state of the art methods for computational fluid dynamics. It is a curious case of how the cellular automata application [1] had evolved to model fluids, and subsequently was proven to be a rigorous numerical scheme for integrating the kinetic Boltzmann equation with the use of Gauss cubatures [2]. From the lattice gas stage of the evolution of the method, it inherited the use of rectangular meshes, which are more common in LBM simulations rather than in other hydrodynamic schemes, such as finite element and discontinuous Galerkin methods [3].

Compared to other methods in computational fluid dynamics (CFD), LBM stands out as a lightweight, efficient method with high locality and many options for parallelism [4,5]. The downside is the limited fluid parameter range where LBM, in its classical formulation, works. It is known to lack Galilean invariance, be working in the range of low Mach numbers, and for isothermal flows. Even in that range, many applications are found, such as mesoscale modeling in additive manufacturing [6]. Moreover, with high popularity of the method, many advanced variations were developed for compressible flows, thermal flows and additional physics, and LBM is used in almost every area of CFD.

In this work, we study the adaptive meshes for LBM. In a simulation of a model with several characteristic scales, where both large-scale background flow and small-scale flow variations have to be reproduced, one can use finer mesh for the areas where it is required, and save the computational cost by coarsening the mesh in the area with lower gradients of fluid properties. Adaptive meshes are also used to cover the areas with fine geometry features, such as thin channels. In the flow around the geometry with high curvature, or flat boundary that is not parallel to the Cartesian axes, the rectangular mesh has to be either sufficiently fine to resolve it, or replaced by a mesh that is conformal with the fluid-solid boundaries in the model. Since the classical LBM is constructed for rectangular meshes with a mesh step and time step equal to unity in non-dimensional units, in all these cases, there is an issue of information transfer between different mesh resolutions, different mesh geometries, and an issue of writing LBM for non-rectangular mesh.

Over the years, many researchers have solved the issues by supplementing LBM with some other numerical method. When mesh is refined in some regions, the information is transferred by interpolation in space and/or time. There are two ways of refining the lattice: node-based and cell-based [7,8]. In a node-based approach, a node is added in between each pair of cells. In a cell-based approach, a lattice node is interpreted as a cubic cell, and each cell is subdivided into smaller cubes.

In the node-based methods [7,9,10,11], the order of interpolation can be increased by using more nodes in the reconstruction of the populations [11,12], or by the use of the information on fluid moments to increase the accuracy of the relevant quantities only [13].

Cell-based methods [8,14] have the mass-conserving property. However, in the conventional method [8], the information transfer from coarse to fine mesh, which is often referred to as ’explode’, is, essentially, an interpolation with a polynomial of the order zero: the values are just copied to the nearest cells. An extension with improved accuracy was proposed [7,14]. The mass-conserving method is used in software like PARAMESH [15], waLBerla [16].

There are other methods to modify LBM so as to work with non-uniform meshes, such as finite difference LBM [17,18,19], finite volume and finite element methods [20,21,22,23,24] which can adapt LBM to curvilinear, tetrahedral and unstructured mesh geometries, Taylor series expansion and least-squares methods [25].

Since in the areas with high space resolution finer time step is desired as well, interpolation in time is also used, however, it can be avoided [26,27].

The motivation of the current study is based on the following observation. As a numerical scheme for the kinetic Boltzmann equation, LBM uses a set of discrete velocities to compute moments of the particle distribution function with Gauss quadrature rules. Various sets can be chosen as long as there are enough points in the quadrature to integrate polynomials of a given order [28]. The discretization of the Boltzmann equation in space-time is performed by integration in the direction of the characteristic [29]. The integration along the characteristic is first order, but, with a change of variables, the second order scheme is obtained. The velocity sets are chosen to form a uniform lattice in space, so that the streaming in the direction of the characteristic moves the particle from one lattice node to another.

What if we choose an arbitrary lattice, that fits the spatial and temporal scales and the geometry of our model, and choose a velocity set for each node in a way that every node is connected to its neighbors? If it is possible, than no other scheme except for the basic LBM is introduced in the model, interpolation is avoided, and the order of accuracy depends only on the proper LBM construction.

Methods, where LBM adapts to a stretched [30] or non-rectangular [31] mesh exist. The fact that a composition of different stencil was not used for adaptive grids before is due to the following issue. A lattice with varied geometry requires different discrete velocity sets in some nodes, and a discrete particle distribution function corresponds to each velocity. These discrete particle distribution functions are the populations that are updated in the LBM. We come to a situation where, in general, there are nodes in which the numerical scheme is constructed differently and operates on a different number of populations, and it is not known how these population sets can be converted into one another.

A well-known method of rescaling was used in the original grid refinement method [9], and improved in further works [32]. This was motivated by the change in the value of the time step, which required change in the collision parameter that controls the viscosity, and, in turn, the non-equilibrium part of the distribution functions. The change in the quadrature rule requires further transformation.

Our solution for the issue is motivated by a success in another area of advanced LBM variations. In [33], a revolutionary method of variable rescaling was introduced for modeling high-Mach flows. Even with classical LBM collision operators and no other modifications, flows with Mach numbers over ${10}^3$ were reproduced. The work had put a start for many new development for complex flow simulations. The scheme was extended and applied for various compressible flows [34,35,36,37,38,39] Recalibrating the populations to the flow velocity temperature for the collision step is known to improve stability [40]

Here we propose to use this method of rescaling for construction of non-uniform grids. Based on the moment matching condition from [33], we develop a population recalibration method that allows to perform streaming between LBM nodes in which the numerical parameters used in construction of LBM method (time and spatial scales, and the velocity set) are different. In general, we identify the numerical parameters that are used to discretize the Boltzmann equations, and find a method to convert a set of LBM populations from a LBM node to an LBM node where LBM was constructed with different parameters.

In the following text, in Section 2, we remind the construction of the LBM method, and what is required to recalibrate LBM variables from one set to another; in Section 3, we propose an algorithm for LBM simulation on a non-uniform grid with the use of the proposed recalibration; in Section 4, we evaluate the resulting method on the trial problems, in Section 5, we discuss the impact of the results.

2. Theoretical Background

2.1. Lattice Boltzmann Method and its Parameters

The kinetic Boltzmann equation is

$$\frac{\partial f}{\partial t}+\mathbf{\xi }\frac{\partial f}{\partial r}=\Omega ,$$

(1)

where $f\left(t,r,\mathbf{\xi }\right)$ — is a single particle distribution function, and its arguments are time, spatial coordinate and particle velocity, and $\Omega$ is the collision term.

Macroscopic flow parameters (density $\rho (t,r)$, flow speed $\mathbf{u}(t,r)$, and temperature $T(t,r)$) are the moments of $f\left(t,r,\mathbf{\xi }\right)$:

$$\rho =\underset{{\mathbb{R}}^D}{\int }f{d}^D\xi ,\rho \mathbf{u}=\underset{{\mathbb{R}}^D}{\int }\mathbf{\xi }f{d}^D\xi ,\rho (u^2+DT)=\underset{{\mathbb{R}}^D}{\int }{\xi }^{2}f{d}^D\xi ,$$

(2)

where D — is the number of spatial dimensions.

In the equilibrium, the collision term is zero, and the solution of (1), which satisfies (2), is

$$f^{\mathrm{eq}}\left(\mathbf{\xi }\right)=\frac{\rho }{{\left(\sqrt{2\pi }{\xi }_0\right)}^D}{e}^{-{\left(\mathbf{\xi }-\mathbf{u}\right)}^2/2{\xi }_0^2},T={\xi }_0^2.$$

(3)

In LBM, the integration in (2) is performed numerically. Taking into account the solution (3), the integration kernel is

$$\omega \left(\xi \right)=\frac{1}{{\left(\sqrt{2\pi }{\xi }_0\right)}^D}{e}^{-{\xi }^2/2{\xi }_0^2}.$$

(4)

Let us take a quadrature rule with the order of accuracy equal to n

$$\underset{{\mathbb{R}}^D}{\int }g\left(\mathbf{\xi }\right)d^D\xi =\underset{{\mathbb{R}}^D}{\int }\frac{g\left(\mathbf{\xi }\right)}{\omega \left(\xi \right)}\omega \left(\xi \right)d^D\xi =/\begin{array}{c}\mathbf{\xi }={\xi }_0\mathbf{v}\\ d^D\xi ={\xi }_0^D{d}^{D}v\end{array}/=\frac{1}{{\left(2\pi \right)}^{D/2}}\underset{{\mathbb{R}}^D}{\int }\frac{g\left({\xi }_0\mathbf{v}\right)}{\omega \left({\xi }_{0}v\right)}{e}^{-v^2/2}{d}^{D}v\simeq \sum _i\frac{w_{i}g\left({\mathbf{c}}_i\right)}{\omega \left(c_i\right)},$$

(5)

where ${\mathbf{c}}_i={\xi }_0{\mathbf{v}}_i$, ${\mathbf{v}}_i$ are the quadrature points, and $w_i$ are the corresponding weights. The expression is exact when $g\left(\mathbf{\xi }\right)$ is a polynomial of the order no higher than n.

With (2) and (5), we find

$$\rho =\sum _i\frac{w_{i}f(t,r,{\mathbf{c}}_i)}{\omega \left(c_i\right)}=\sum _i{f}_i,$$

(6)

where the discrete populations $f_i$ are defined as

$$f_i(t,r)\equiv \frac{w_{i}f(t,r,{\mathbf{c}}_i)}{\omega \left(c_i\right)}={\left(\sqrt{2\pi }{\xi }_0\right)}^D{w}_{i}f(t,r,{\mathbf{c}}_i)e^{c_i^2/2{\xi }_0^2}.$$

(7)

Other moments are obtained in a similar fashion:

$$\rho \mathbf{u}=\sum _i{f}_i{\mathbf{c}}_i,\rho (u^2+DT)=\sum _i{f}_i{c}_i^2.$$

(8)

The equilibrium discrete populations are

$$f_i^{\mathrm{eq}}=\rho w_i{e}^{(2{\mathbf{c}}_i·\mathbf{u}-u^2)/2{\xi }_0^2}.$$

(9)

The standard polynomial representation of $f_i^{\mathrm{eq}}$ [29] can be obtained by expanding this expression in the Taylor series:

$$f_i^{\mathrm{eq}}=\rho w_i\left(1-\frac{u^2}{2{\xi }_0^2}+\frac{{\mathbf{c}}_i·\mathbf{u}}{{\xi }_0^2}+\frac{{({\mathbf{c}}_i·\mathbf{u})}^2}{2{\xi }_0^4}\right).$$

(10)

By inserting (10) and $T={\xi }_0^2$ in (6), (8) the moments are obtained exactly if the chosen quadrature rule is at least of the fourth order of accuracy.

In this derivation, ${\xi }_0$ has been introduced as a scheme parameter in (4). It allows to scale ${\mathbf{c}}_i$, therefore, it provides a connection between time and space scales and controls the lattice geometry. Taking the expression for $f^{\mathrm{eq}}$ into account, it is also associated with the model temperature and the speed of sound. In the current work, athermal LBM is considered, and ${\xi }_0$ remains a numerical parameter in the scheme construction. For athermal LBM, the 5-th order quadrature is enough [28].

The LBM equations are obtained by letting $\mathbf{\xi }={\mathbf{c}}_i$ in (1) and multiplying by $w_i/\omega \left(c_i\right)$ for each i:

$$\frac{\partial f_i}{\partial t}+{\mathbf{c}}_i\frac{\partial f_i}{\partial r}={\Omega }_i,{\Omega }_i\equiv \frac{w_i}{\omega \left(c_i\right)}{\left.\Omega \right|}_{\mathbf{\xi }={\mathbf{c}}_i}.$$

(11)

In this work, the collision operator is taken in the Bhatnagar-Gross-Krook (BGK) form [41,42]:

$$\Omega =\frac{f^{\mathrm{eq}}-f}{\tau },{\Omega }_i=\frac{f_i^{\mathrm{eq}}-f_i}{\tau },$$

(12)

where $\tau$ is the collision parameter that controls the relaxation rate of populations and the fluid viscosity.

Let us take a uniform lattice with mesh step $\Delta x$. Then, (11) is discretized on a time step $\Delta t$ in two successive steps:

$$\begin{array}{cc}•\mathrm{streaming}\text{:}\frac{\partial f_i}{\partial t}+{\mathbf{c}}_i\frac{\partial f_i}{\partial r}=0\Rightarrow \hfill & \\ & \hfill f_i(t+\Delta t,r+\Delta r)-f_i(t,r)=0,\Delta r\equiv {\mathbf{c}}_i\Delta t;\end{array}$$

(13)

$$\begin{array}{cc}•\mathrm{collision}\text{:}\frac{\partial f_i}{\partial t}=\frac{f_i^{\mathrm{eq}}-f_i}{\tau }\Rightarrow \hfill & \\ & \hfill f_i(t+\Delta t,r)-f_i(t,r)=\frac{\Delta t}{\tau }\left(f_i^{\mathrm{eq}}(t,r)-f_i(t,r)\right).\end{array}$$

(14)

In (13), $\Delta r$ should point from one node to another, therefore, every velocity in the set $\left\{{\mathbf{c}}_i\right\}$ should be scaled accordingly. The ${\xi }_0$ parameter is used to scale $\left\{{\mathbf{c}}_i\right\}$. For example, the widely used D2Q9 set, which is

$$\left\{{\mathbf{v}}_i^{\mathrm{D}2\mathrm{Q}9}\right\}=\left\{(0,0),(0,\pm \sqrt{3}),(\pm \sqrt{3},0),(\pm \sqrt{3},\pm \sqrt{3})\right\},$$

(15)

corresponds to ${\xi }_0=\Delta x/(\sqrt{3}\Delta t)$.

Thus, according to (13), (14), LBM method consists of the two steps: streaming from one node to another in the direction of ${\mathbf{c}}_i$, and a local collision operation.

Let us define the LBM stencil as a complete set of scheme parameters that uniquely define LBM streaming: quadrature points $\left\{{\mathbf{v}}_i\right\}$ and weights, scaling coefficient ${\xi }_0$, time step $\Delta t$. We denote the stencil by inserting these parameters in the brackets, i.e. $\mathrm{D}N\mathrm{Q}M(\Delta t,{\xi }_0^2)$. The common D2Q9 LBM stencil (15) with $\Delta t=\Delta x=1$ is denoted as $\mathrm{D}2\mathrm{Q}9(1,1/3)$ hereafter.

The fluid viscosity $\nu$ is derived with the Chapman-Enskog analysis as [29]:

$$\nu ={\xi }_0^2\left(\tau -\frac{\Delta t}{2}\right).$$

(16)

Finally, given a lattice and a fluid with viscosity equal to $\nu$, to use classical athermal LBM with BGK collision, we have to ensure that each vector ${\xi }_0{\mathbf{v}}_i\Delta t$ points exactly from one lattice node to another, the quadrature set is enough to integrate with the 5th order of accuracy, and (16) is satisfied.

2.2. Recalibration of Populations

In LBM with adaptive grids, lattice nodes may have varying stencils. For the information exchange between LBM nodes that operate with different LBM stencils, we need to know how the population sets are converted into one another. To construct the method of recalibration of populations from one stencil to another, let us find out what changes are required with the change of $\Delta t$, ${\xi }_0$, and the set of quadrature points.

Let us remark that we allow the stencil, and, thus, ${\mathbf{c}}_i$, vary in time and space, and the derivation of (11) requires inserting $\omega \left({\mathbf{c}}_i\right)$ under the derivative [43]. There is no contradiction here. In each node, we construct LBM as if it is on a uniform grid. In the following, we find out how the information is converted between LBM constructed for different stencils.

2.2.1. Recalibration with $\Delta t$

According to (13), (14), the equations of LBM are:

$$f_i(t+\Delta t,r+{\mathbf{c}}_i\Delta t)=f_i(t,r)+\frac{\Delta t}{\tau }\left(f_i^{\mathrm{eq}}(t,r)-f_i(t,r)\right).$$

(17)

In other words, the populations are collided in $(t,r)$ and travel to $(t+\Delta t,r+{\mathbf{c}}_i\Delta t)$.

To find the relationship between population which travel from the node $(t,r)$ according to different stencils, let us use Chapman-Enskog analysis. Taking into account only the term that is linear in the Knudsen number $ϵ$,

$$f_i\simeq f_i^{\mathrm{eq}}+ϵf_i^{\left(1\right)},f_i^{\left(1\right)}=-\tau \left({\partial }_t^{\left(1\right)}{f}_i^{\mathrm{eq}}+{\mathbf{c}}_i{\partial }_r^{\left(1\right)}{f}_i^{\mathrm{eq}}\right),{\partial }^{\left(1\right)}=\frac{\partial }{ϵ},$$

(18)

and

$$f_i\simeq f_i^{\mathrm{eq}}-\tau \left({\partial }_t{f}_i^{\mathrm{eq}}+{\mathbf{c}}_i{\partial }_r{f}_i^{\mathrm{eq}}\right)=f_i^{\mathrm{eq}}-\tau D{f}_i^{\mathrm{eq}},D\equiv {\partial }_t+{\mathbf{c}}_i{\partial }_r.$$

(19)

Let us consider conversion between stencils with differ only in $\Delta t$: D2QM($\Delta t_c,{\xi }_0$) and D2QM($\Delta t_f,{\xi }_0$). Let us take a coarse stencil with time step $\Delta t_c$ and relaxation parameter ${\tau }_c$, and a fine stencil with $\Delta t_f$ and ${\tau }_f$. By inserting (19) into (17) we get

$$f_{i,c}(t+\Delta t_c,r+{\mathbf{c}}_i\Delta t_c)-f_i^{\mathrm{eq}}(t,r)=(\Delta t_c-{\tau }_c)D{f}_i^{\mathrm{eq}}(t,r).$$

(20)

and a similar relation for the fine stencil. By equating $D{f}_i^{\mathrm{eq}}(t,r)$ in the expression for the fine and the coarse stencils, we get a well-known relationship [9]

$$\frac{f_{i,c}(t+\Delta t_c,r+{\mathbf{c}}_i\Delta t_c)-f_i^{\mathrm{eq}}(t,r)}{f_{i,f}(t+\Delta t_f,r+{\mathbf{c}}_i\Delta t_f)-f_i^{\mathrm{eq}}(t,r)}=\frac{\Delta t_c-{\tau }_c}{\Delta t_f-{\tau }_f},$$

(21)

and the relationship between ${\tau }_c$ and ${\tau }_f$ is obtained be requiring the invariance of fluid viscosity (16) .

$$\nu ={\xi }_0^2\left({\tau }_c-\frac{\Delta t_c}{2}\right)={\xi }_0^2\left({\tau }_f-\frac{\Delta t_f}{2}\right)\Rightarrow \frac{2{\tau }_c-\Delta t_c}{2{\tau }_f-\Delta t_f}=1.$$

(22)

2.2.2. Recalibration with both $\Delta t$ and ${\xi }_0$

Let us take two stencils DNQM$(\Delta t_f,{\xi }_{0,f}^2)$ and DNQM$(\Delta t_c,{\xi }_{0,c}^2)$ which correspond to the same quadrature $\left\{{\mathbf{v}}_i\right\}$, but different ${\xi }_0$, $\Delta t$, $\tau$. We have

$$\left\{{\mathbf{c}}_{i,f}\right\}\ne \left\{{\mathbf{c}}_{i,c}\right\},D_f\ne D_c,f_{i,f}^{\mathrm{eq}}\ne f_{i,c}^{\mathrm{eq}},$$

(23)

and, instead of (22),

$$\nu ={\xi }_{0,c}^2\left({\tau }_c-\frac{\Delta t_c}{2}\right)={\xi }_{0,f}^2\left({\tau }_f-\frac{\Delta t_f}{2}\right)\Rightarrow \frac{2{\tau }_c-\Delta t_c}{2{\tau }_f-\Delta t_f}=\frac{{\xi }_{0,f}^2}{{\xi }_{0,c}^2}.$$

(24)

Taking (23) into account, (20) leads to the following expression:

$$\frac{f_{i,c}(t+\Delta t,r+{\mathbf{c}}_{i,c}\Delta t)-f_{i,c}^{\mathrm{eq}}(t,r)}{f_{i,f}(t+\Delta t,r+c_{i,f}\Delta t)-f_{i,f}^{\mathrm{eq}}(t,r)}=\frac{(\Delta t_c-{\tau }_c)D_c{f}_{i,c}^{\mathrm{eq}}}{(\Delta t_f-{\tau }_f)D_f{f}_{i,f}^{\mathrm{eq}}}.$$

(25)

Let us find how $D{f}_i^{\mathrm{eq}}$ depends on the stencil. Let us express it through the derivatives of the physical quantities $\rho$, $\mathbf{u}$. From (10) and (19), we get

$$\begin{array}{c}D{f}_i^{\mathrm{eq}}=w_iϵ\left({\partial }_t^{\left(1\right)}\rho -\frac{1}{2{\xi }_0^2}{\partial }_t^{\left(1\right)}\rho u^2+\frac{c_{i,\beta }}{{\xi }_0^2}{\partial }_t^{\left(1\right)}\rho u_{\beta }+\frac{c_{i,\beta }{c}_{i,\gamma }}{2{\xi }_0^4}{\partial }_t^{\left(1\right)}\rho u_{\beta }{u}_{\gamma }\right.+\hfill \\ \hfill \left.c_{i,\alpha }{\partial }_{\alpha }^{\left(1\right)}\rho -\frac{c_{i,\alpha }}{2{\xi }_0^2}{\partial }_{\alpha }^{\left(1\right)}\rho u^2+\frac{c_{i,\alpha }{c}_{i,\beta }{c}_{i,\gamma }}{2{\xi }_0^4}{\partial }_{\alpha }^{\left(1\right)}\rho u_{\beta }{u}_{\gamma }+\frac{c_{i,\alpha }{c}_{i,\beta }}{{\xi }_0^2}{\partial }_{\alpha }^{\left(1\right)}\rho u_{\beta }\right).\end{array}$$

(26)

Time derivatives in (26) are replaced by the space derivatives with the use of the conservation laws

$${\partial }_t^{\left(1\right)}\rho +{\partial }_{\alpha }^{\left(1\right)}\left(\rho u_{\alpha }\right)=0,{\partial }_t^{\left(1\right)}\left(\rho u_{\alpha }\right)+{\partial }_{\gamma }^{\left(1\right)}{\Pi }_{\alpha \gamma }^{\mathrm{eq}}=0,$$

(27)

where

$${\Pi }_{\alpha \gamma }^{\mathrm{eq}}=\rho u_{\alpha }{u}_{\gamma }+\rho {\xi }_0^2{\delta }_{\alpha \gamma }.$$

(28)

The time derivative is expanded,

$${\partial }_t^{\left(1\right)}\rho u_{\beta }{u}_{\gamma }=u_{\beta }{\partial }_t^{\left(1\right)}\rho u_{\gamma }+u_{\gamma }{\partial }_t^{\left(1\right)}\rho u_{\beta }-u_{\beta }{u}_{\gamma }{\partial }_t^{\left(1\right)}\rho .$$

(29)

Finally, with the use of (27) and (28), (26) becomes

$$\begin{array}{c}D{f}_i^{\mathrm{eq}}=w_iϵ\left(-\rho {\partial }_{\alpha }^{\left(1\right)}{u}_{\alpha }+\frac{c_{i,\alpha }{c}_{i,\beta }}{{\xi }_0^2}\rho {\partial }_{\alpha }^{\left(1\right)}{u}_{\beta }\right.-\left.\frac{c_{i,\beta }}{{\xi }_0^2}{\partial }_{\gamma }^{\left(1\right)}\rho u_{\beta }{u}_{\gamma }-\frac{c_{i,\alpha }}{2{\xi }_0^2}{\partial }_{\alpha }^{\left(1\right)}\rho u^2+\frac{c_{i,\alpha }{c}_{i,\beta }{c}_{i,\gamma }}{2{\xi }_0^4}{\partial }_{\alpha }^{\left(1\right)}\rho u_{\beta }{u}_{\gamma }+\mathcal{O}\left(u^3\right)\right).\hfill \end{array}$$

(30)

According to (18), here, we can use ∂ instead of ${\partial }^{\left(1\right)}$, and insert ${\mathbf{c}}_i={\xi }_0{\mathbf{v}}_i$ to find how $D{f}_i^{\mathrm{eq}}$ depends on ${\xi }_0$.

$$\begin{array}{c}D{f}_i^{\mathrm{eq}}=w_i\left(-\rho {\partial }_{\alpha }{u}_{\alpha }+v_{i,\alpha }{v}_{i,\beta }\rho {\partial }_{\alpha }{u}_{\beta }\right.-\left.\frac{v_{i,\beta }}{{\xi }_0}{\partial }_{\gamma }\rho u_{\beta }{u}_{\gamma }-\frac{v_{i,\alpha }}{2{\xi }_0}{\partial }_{\alpha }\rho u^2+\frac{v_{i,\alpha }{v}_{i,\beta }{v}_{i,\gamma }}{2{\xi }_0}{\partial }_{\alpha }\rho u_{\beta }{u}_{\gamma }+\mathcal{O}\left(u^3\right)\right).\hfill \end{array}$$

(31)

Here the underlined terms do not depend on ${\xi }_0$. These terms are unchanged as long as the same quadrature rule is used. The following terms can be neglected in the low compressibility limit. Therefore, the recalibration expression is similar to (21):

$$\frac{f_{i,c}(t+\Delta t,r+{\mathbf{c}}_{i,c}\Delta t)-f_{i,c}^{\mathrm{eq}}(t,r)}{f_{i,f}(t+\Delta t,r+{\mathbf{c}}_{i,f}\Delta t)-f_{i,f}^{\mathrm{eq}}(t,r)}=\frac{\Delta t_c-{\tau }_c}{\Delta t_f-{\tau }_f},$$

(32)

and the relationship between ${\tau }_f$ and ${\tau }_c$ is (24).

2.2.3. Recalibration with a change in quadrature

Let us construct the conversion between populations in D2Q$M_1$$(\Delta t,{\xi }_0^2)$ and in D2Q$M_2$$(\Delta t,{\xi }_0^2)$.

The conversions in (21), (32) are between the stencils of the same quadrature $\left\{{\mathbf{v}}_i\right\}$, but different $\tau$. Therefore, these are the conversions between LBMs with different collision operators, which are essentially different LBM equations, and the conversions are giving the correct non-equilibrium part of the populations.

When the $\tau$ parameter remains the same, the moment matching condition [33,34,35,36] is used for conversion. It is required that the change in the LBM populations does not lead to the change in the physical moments:

$${\left.\left(\sum _i{c}_{i,x}^p{c}_{i,y}^q{f}_i\right)\right|}_{\mathrm{D}2\mathrm{Q}{M}_1(\Delta t,{\xi }_0^2)}={\left.\left(\sum _i{c}_{i,x}^p{c}_{i,y}^q{f}_i\right)\right|}_{\mathrm{D}2\mathrm{Q}{M}_2(\Delta t,{\xi }_0^2)}.$$

(33)

The number of equations in the linear system (33) depends on the order of approximation of the quadrature and its symmetries [28,29,44]. For athermal fluid physics, we require at least

$$p+q\le 5p,q\in {\mathbb{N}}_0.$$

(34)

Additionally, an explicit relationship can be found for $f_0$ (which corresponds to the zero discrete velocity vector) with (7):

$${\left.(f_0/w_0)\right|}_{\mathrm{D}2\mathrm{Q}{M}_1(\Delta t,{\xi }_0^2)}={\left.(f_0/w_0)\right|}_{\mathrm{D}2\mathrm{Q}{M}_2(\Delta t,{\xi }_0^2)}.$$

(35)

Similar relationship can be written whenever stencils have a share discrete velocity.

Relations (33), (35) are a system of linear equations, that has to be solved to find the populations of one stencil, while the populations of another stencil are known. In general, the system can be underdetermined or inconsistent. For example, in the conversion from D2Q9 to D2Q15 we have to find 15 unknown populations from computable moments. In this case, we append the systems by equating relationships from (33) to the equilibrium moments [45].

2.2.4. Recalibration with the change of stencil

In total, in the recalibration of populations from D2Q$M_1$$(\Delta t_1,{\xi }_{0,1}^2)$ to D2Q$M_2$$(\Delta t_2,{\xi }_{0,2}^2)$ two conversions are applied.

First, the quadrature rule is fixed, and $\Delta t$, $\tau$, ${\xi }_0$ are modified with (24), (32). This converts D2Q$M_1$$(\Delta t_1,{\xi }_{0,1}^2)$ to D2Q$M_1$$(\Delta t_2,{\xi }_{0,2}^2)$. Second, the quadrature is changed with fixed $\tau$, ${\xi }_0$, and the stencil is converted to D2Q$M_2$$(\Delta t_2,{\xi }_{0,2}^2)$. The two steps can be performed in any order. In any case, in the middle of the conversion, there is a virtual LBM stencil (here, D2Q$M_1$$(\Delta t_2,{\xi }_{0,2}^2)$) for which neither streaming nor collision takes place.

3. Grid Refinement Interface without Interpolation

3.1. Grid Geometry

The recalibration method allows to construct a transition between grids with different time and space steps, if only we construct the stencils that point to/from the lattice points exactly in the uneven grid transition region.

Let us take a grid refinement boundary with coarse uniform grid ($\Delta t_c=\Delta x_c=1$) on the left side and fine uniform grid ($\Delta t_f=\Delta x_f=1/2$) on the right side (Figure 1).

LBM stencils can be used in the ’Pull’ or ’Push’ paradigm, which is relevant for the conversion purposes. In the ’Push’ method, the populations are streamed from a node according to the stencil of that node. After the streaming, the populations in any one node may come from different stencils, and there is no full set of populations in some nodes. In the ’Pull’ method, populations are converted into the stencil of the node before they are streamed into that node. The ’Pull’ method is in the original work [33], and is used in the current study. In case of the uniform grids, no stencil conversions are used and the ’Pull’ and ’Push’ methods are equivalent. The difference of rescaling incoming and outgoing populations in refined grids supplemented with interpolation is reported in [32].

3.2. Stencils and Recalibration

Let us construct an LBM algorithm, in which coarse and fine grids use the classical LBM stencil, and the transition between grids do not require interpolated values, that is, all ${\mathbf{c}}_i$ point from one node to another. This is possible by introducing several transitional LBM stencils.

Let us use D2Q9(1, 1/3) (dark-green in Figure 1) and D2Q9(1/2, 1/3) (light-green) for the coarse and the fine grids respectively. This leads to the classical unaltered D2Q9 LBM in the fine and coarse grids away from the boundary.

For the grid transition, the well-known quadrature derivation rules can be used [46]. D2Q15(1, 25/38) (see Appendix A for the derivation) is a stencil of the 5th order which uses the existing nodes, thus, its can reproduce the moments of the equilibrium distribution at least as good as D2Q9 which is used away from the grid transition. We have also tried the D2Q7(1, 1/4) (details in Appendix A). This stencil performs correct integration for less moments, but is more local. The tests with both stencils are reported in Section 4. In Figure 1, D2Q7(1, 1/4) is depicted, but D2Q15(1, 25/38) can be used instead.

To perform grid transition, we use nodes with the D2Q7(1, 1/4) (or D2Q15(1, 25/38)) stencil (blue in Figure 1) on the integer time steps, and the D2Q9(1/2, 4/3) on half-integer time steps (orange nodes in Figure 1). This is just one of the many possibilities that can be constructed on the boundary transition.

Other configurations can be used depending on the requirements of the model.

When the introduced stencils are used in the ’Pull’ paradigm, the streaming operation at a node requests populations from the neighboring nodes. At the node from where the population is requested from, the recalibration in the target stencil takes place. If it is a recalibration between different variations of D2Q9, (24) is used. If any other stencil is involved, the recalibration is performed in two steps according to Section 2.2.4.

$$\mathrm{D}2\mathrm{Q}9(\cdots )\rightleftarrows \mathrm{D}2\mathrm{Q}9(1,1/4)\rightleftarrows \mathrm{D}2\mathrm{Q}7(1,1/4),$$

(36)

$$\mathrm{D}2\mathrm{Q}9(\cdots )\rightleftarrows \mathrm{D}2\mathrm{Q}9(1,25/38)\rightleftarrows \mathrm{D}2\mathrm{Q}15(1,25/38),$$

(37)

Here D2Q9(1, 1/4) and D2Q9(1, 25/38) are the virtual stencils that are used only as an intermediate state in the conversion of populations.

3.3. Full Grid Transition Algorithm

The initial state $t=t_0$ corresponds to Figure 1(a).

• Perform streaming on the coarse grid $t_0\to t_0+\Delta t_c$ with the use of the stencils

  (a)

  D2Q7(1, 1/4) (or D2Q15(1, 25/38)) for the blue nodes;

  (b)

  D2Q9(1, 1/3) for dark green nodes. Here the incoming populations at the nodes which are exactly on the boundary are saved in a separate temporary buffer to be used in step 5, since the pre-streaming populations are still needed in the next step.

• Perform streaming on the fine grid $t_0\to t_0+\Delta t_f$ (Fig Figure 1(b)) with the use of the stencils

  (a)

  D2Q9(1/2, 1/3) at the light-green nodes;

  (b)

  D2Q9(1/2, 4/3) at the orange nodes.

• Perform collision on the fine grid at the orange and light-green nodes (Fig Figure 1(b)).
• Perform second streaming $t_0+\Delta t_f\to t_0+\Delta t_c$ to the light-green nodes of the fine grid (Fig Figure 1(a)).
• Restore the values of the boundary nodes from the buffer.
• Perform collision on all nodes according to the stencils in Fig Figure 1(a).

4. Benchmarks

The proposed algorithm was implemented in C++ with the use of the Zipped Data Structure for Adaptive Mesh Refinement (ZAMR) [47] library in the Aiwlib [48] package. The library provides convenient tools to work with binary refined grids based on the Z-order curve traversal. It allows to perform operations on all nodes of uneven grids, and to set flags on the nodes where operations differ. With it, we mark different color of nodes in Figure 1 with different flags and implement the described algorithm.

4.1. Poiseuille Flow

The stationary solution of the Navier-Stokes equation

$$\frac{\partial \mathbf{u}}{\partial t}+(\mathbf{u}\nabla )\mathbf{u}=-\frac{1}{\rho }\nabla p+\nu \Delta \mathbf{u}$$

(38)

in the $-H_x/2\le x\le H_x/2$, $0\le y\le H_y$ domain with no-slip boundary in x and periodic boundary in y:

$$\mathbf{u}\left(-\frac{H_x}{2},y\right)=\mathbf{u}\left(\frac{H_x}{2},y\right)=\mathbf{0},\mathbf{u}(x,y+H_y)=\mathbf{u}(x,y),\nabla p=(0,-g),$$

(39)

subject to the pressure gradient

$$\nabla p=(0,-g),$$

(40)

is

$$u_x=0,u_y\left(x\right)=\frac{g{H}_x^2}{8\rho \nu }\left(1-\frac{4x^2}{H_x^2}\right).$$

(41)

The flow is modeled with the use of the bounce-back boundary conditions in x [29]. The pressure gradient is simulated with an additional term $\Delta u_y=g\tau /\rho$ in the flow velocity in the computation of equilibrium distribution [49] for the collision step. The initial state is equilibrium where

$$\rho =1,u_x=u_y=0.$$

(42)

The flow evolves to its equilibrium state (41).

In this benchmark, classical LBM with BGK collision and bounce-back boundary has a well-known problem. The effective channel width $H_x$ can depend on the collision parameter $\tau$ [50]. In this benchmark we fix ${\tau }_f=(\sqrt{3}+2)\Delta t_f/4$ on the fine grid. Then the viscosity is $\nu =\sqrt{3}\Delta t_f/12$, and the boundary is considered to be at the distance $\Delta x_f/2$ from the boundary nodes [29].

We refine the grid near the no-slip boundaries. The sizes of the grid are

$$H_x^r=2^{3+r}\Delta x_c+\Delta x_f,H_b^r=2^r\Delta x_c+\frac{\Delta x_f}{2}\approx \frac{H_x^r}{8},H_y=4\Delta x_c.$$

(43)

Where r is the parameter that controls the grid resolution. With $r=1$, half of the grid geometry is in Figure 2.

The fine grid takes approximately a quarter of the domain. In our benchmarks, such grid geometry performed slightly better than if the grid is refined in the half of the simulation region.

To evaluate the order of approximation of the proposed scheme, the numerical error dependency on r is studied. With the change of r with fixed $\Delta x$, the effective channel width and velocity maximum increase. Thus, the error is normalized to ${(H_x^0/H_x^r)}^2$. Note that in evaluating the $L_p$ norm in case of non-uniform grids, the contribution of nodes is scaled to the effective area of the node $dS\left(\mathrm{node}\right)$ to ensure correct numerical representation of the finite integral.

$$L_{\infty }\left(H_x^r\right)={\left(\frac{H_x^0}{H_x^r}\right)}^2\underset{\mathrm{nodes}}{max}|u_y^{r,\mathrm{num}}-u_y^r|,$$

(44)

$$L_n\left(H_x^r\right)={\left(\frac{H_x^0}{H_x^r}\right)}^2\sqrt[n]{\frac{\sum _{\mathrm{nodes}}{|u_y^{r,\mathrm{num}}-u_y^r|}^{n}dS\left(\mathrm{node}\right)}{\sum _{\mathrm{nodes}}dS\left(\mathrm{node}\right)}},n\in \{1,2\},$$

(45)

$$dS\left(\mathrm{coarse}\right)=\Delta x_c^2,dS\left(\mathrm{fine}\right)=\Delta x_f^2,dS\left(\mathrm{border}\right)=(\Delta x_c^2+\Delta x_f^2)/2,$$

(46)

The blue nodes in Fig Figure 2 serve auxiliary purpose in the scheme and are not taken into account in the total error.

The benchmark results are reported in Figure 3, Figure 4. The results are compared to LBM simulation on a uniform fine grid (light green lines).

We observe that the proposed method of grid transition works correctly.

The D2Q7(1, 1/4) stencil performs better in the space time discretization, and gives the result almost as good as a uniform fine grid. However, the stencil does not reproduce high order velocity moments, and suffers from Galilean invariance errors. That is why its performance with grid refinement is more consistent with lower g (Figure 4).

The discretization errors for the D2Q15(1, 25/38) are higher and show the first order of accuracy. The source of the error is yet to be found. The dependencies of this stencil in the space-time grid are the longest of all other stencils in the constructed scheme. The low angle of the space-time characteristics may cause such behavior.

Let us compare the order of accuracy of our algorithm with the already known schemes for constructing non-uniform grids. When using stencil D2Q7(1, 1/4), we get almost the second order of accuracy in each norm, which is more stable than both algorithms with initial-value problem (IVP) and boundary-value problem (BVP) types of interpolation in [26], and similar to cell-vertex and cell-centered with linear explosion hybrid-recursive regularized (HRR) algorithms in [7]. When using stencil D2Q15(1, 25/38), we get the order of accuracy close to unity, which is resemble method with IVP interpolation in [26] and cell-centred with uniform explosion HRR algorithm in [7].

Let us note that the benchmark results can be improved by the use of advanced collision and boundary conditions [29,50]. Here we demonstrate only the impact of the proposed recalibration method, thus, the classical LBM with BGK collision is used.

In Figure 5, the difference between numerical and theoretical solution is plotted vs the x coordinate. Here $r=1$, $g=1\times {10}^{-6}$, $\nu =\sqrt{3}/24$. It can be seen that the grid transition introduced an error an leads to the non-physical current in the x direction, and the error, indeed, appears on the grid boundary, and is not caused by the bounce-back boundary.

4.2. Athermal vortex

In the second benchmark, the dynamics of the advected athermal vortex [51,52] is modelled. The initial conditions are

$${\left.u_{\phi }\left(r\right)\right|}_{t=0}=\frac{r{\beta }_0}{2\pi R}exp\left(\frac{1-r^2/R^2}{2}\right),{\left.\rho \left(r\right)\right|}_{t=0}=exp\left(-\frac{{\beta }_0^2}{8{\pi }^{2}T}exp\left(1-r^2/R^2\right)\right),$$

(47)

for the stationary vortex, and a constant $u_0$ is added to $u_x$ to model advection. Here ${\beta }_0$ is the amplitude of the vortex rotation, R is the vortex radius, T is the lattice temperature. According to the Navier-Stokes solution, the vortex relaxes to uniform flow, and the relaxation rate is proportional to viscosity. In our benchmarks, we set lower viscosity and study the preservation of the vortex shape under the influence of the grid boundaries.

We take the simulation region size $-512\le x\le 512$, $0\le y\le 512$, the left half of the region is refined. All boundaries are periodic. The vortex is initialized in the left half by setting the equilibrium populations on the fine grid with the flow parameters according to (47). The numerical parameters are $\nu =0.01$, ${\beta }_0=0.05$, $R=40$, $T=1/3$, $u_0=0.1$.

The D2Q15(1, 25/38) stencil is used on the grid transition. Figure 6 shows the moment of vortex interaction with the grid boundary. Its shape is preserved in the interaction, and even after it passes back onto the fine grid (Figure 7). The effect of the transition on the solution can be seen by magnifying the density range on the color map (Figure 6(b)). The effect is several orders below the solution.

5. Conclusion

As a conclusion, we have shown how the moment matching condition [11] can be used to perform conversions between LBM stencils on non-uniform grids. For LBM with BGK collision, we developed the method of recalibration of populations, which extends the one proposed in [9] to variable lattice temperatures, and supplemented it with the method of a transition to a different discrete velocity set.

We proved the recalibration to be valid with 2D benchmarks with one-dimensional grid transition interface. For this purpose, the algorithm that involves intermediate stencils on the transition interface was proposed.

We believe that this is an important contribution because it allows building non-uniform meshes while remaining in the LBM framework, without the use of interpolation, finite differences, series expansion, or other supplementary schemes. The recalibration method is based only on the correctness of quadrature rule in the velocity space, and the Chapman-Enskog analysis. This method can be used for different configuration of mesh interfaces, and even for dynamically refined or moving grids [53].

To apply the method further for more complex grid configurations, a combination of stencils for the grid interface has to be built. After that, the proposed method can be applied to perform the transitions.

The choice of the best configuration is yet to be found. In the following studies, we may search the source of errors that led to worse order of approximation in the Poiseuille flow benchmark.

Furthermore, it is interesting to formulate the recalibration in the ’Push’ framework based on the rescaling of the outgoing populations [32] and moment matching with more than two population sets. By combining a ’Pull’ and ’Push’ types of streaming, a configuration of stencils that lead to a grid transition without the mass loss may be found.

In the work we used the basic BGK collision and the basic bounce-back boundary to focus on the impact of the recalibration method. For other collision operators, the recalibration procedure, and the adaptation of the collision parameter to the stencil, may differ, but may be built on the same principles that were discussed in the current text.

While the classical LBM is known to be limited in the range of its applicability, it is still valid in the vicinity of the equilibrium function with fixed velocity (usually, zero), low Mach number, and on a sufficiently detailed space-time lattice. With the use of the moment matching condition, the LBM can be locally rebuilt with different background flow [33,54,55], with different compressibility assumptions [45]. Thus, the LBM stencil can change both in time and in space to adapt to the flow conditions. The contribution of our work is the demonstration of the fact that the moment matching can also be used to adapt the stencil to the local geometry, or to an increased Reynolds flows.