Exploiting decoupled discretization in meshless lattice Boltzmann method for high-Reynolds number flows

Dawid Strzelczyk^∗, Maciej Matyka

Institute of Theoretical Physics, Faculty of Physics and Astronomy, University of Wrocław, pl. M. Borna 9, 50-204, Wrocław, Poland

^∗ dawid.strzelczyk@uwr.edu.pl

Abstract One of the most severe limitations of the Lattice Boltzmann Method in the context of simulating inertial flows is the coupling of the discretization of space to the velocity discretization. It causes the need to increase the size of computational meshes whenever one wants to increase the Reynolds number in the system for a fixed velocity and viscosity. This work aims at adopting the recently proposed meshless LBM formulation to the problem of high Reynolds number flows by using its internal property of decoupled space and velocity discretizations. In meshless formulation, one can change the reference length of the problem by scaling the distances between points that discretize the domain while the streaming distance remains fixed. By doing that, we increase the Reynolds number at no cost: additional nodes need not be inserted in the discretization. We measure the accuracy and efficiency of this approach in the flow around a circular obstacle at $Re=84$–$168$, the lid-driven cavity flow at $Re=100$–$5000$, and the flow through a porous sample at $Re=0.8$–$1500$. Additionally, we apply the meshless streaming step to the recently proposed fixed relaxation time $\tau=1$ LBM to extend its applicability to model inertial flows.

Keywords: Lattice Boltzmann Method, meshless methods, inertial flows, Reynolds number

1 Introduction

The transport of fluids through complex, porous structures is ubiquitous with many technological applications (CO2 storage, oil and shale gas extraction, electric battery design). It is of great interest for the study of human health, e.g. blood and cerebrospinal fluid flow, gas transport in the lungs[1, 2].

Experimental investigation of fluid flow phenomena in porous media is costly and requires sophisticated setups. Measuring velocity fields in creeping flows using Particle Image Velocimetry (PIV) methods has become popular [3, 4]. However, beyond Darcy Law (inertial regime) [5], PIV tracers used to track the flow face several problems. They may, for instance, damage the porous sample, and some of the flow regions may get excluded from the measurement due to the lack of tracers therein. Also, the tracers’ size must be chosen carefully to reconcile the accuracy of the scattered signal acquisition and the minimal interaction between the tracers, flow, and the pore space [6]. The remedy for this is to use numerical methods to solve fluid flow transport equations directly in pore space. Real samples, however, are relatively big compared to pore size, and using standard numerical methods based on computational grids may become impractical. This becomes most problematic at large Reynolds numbers, where fine flow structures appear.

Recently, the Lattice Boltzmann Method (LBM) became one of the most popular numerical tools used to solve the Navier-Stokes equations at the pore scale [7, 8]. It was applied to the simulation of fluid transport in various contexts of porous media [9], multiphase flows [10], semiclassical fluids flow [11], particulate suspensions [12] and especially inertial flows, and turbulent flows [13, 14]. Its popularity arises from the stable and local treatment of non-linear effects, a huge potential for parallelization of calculations, and simple implementation of basic no-slip boundary conditions on a regular grid [15, 16, 17]. The latter is especially attractive in complex porous media geometries, i.e. those arising from tomography scans in the form of binary images. Unfortunately, the memory overhead associated with storing the distribution function in LBM, and the need to regularly discretize the entire simulation space, limits the utility of this method for real-life, large samples with relatively small pore sizes.

One way to overcome the mentioned problem of excessive memory usage is to resort to LBM models with the fixed relaxation time $\tau=1$ (LBM1) [18, 19] where the distribution function is not explicitly stored in the memory. As shown in our previous work, LBM1 simplifies the main algorithm and reduces the memory requirement up to $86\%$ compared to the standard LBM solver [18]. Savings may also come from the use of the meshless discretization of space in LBM (meshless LBM, MLBM) [20]. This offers a significant advantage regarding the flexibility of space discretization, with LBM grid points distributed, i.e., in the pore-space only. We recently investigated its convergence in several scenarios, including low-Reynolds number flow through porous samples consisting of a network of two-dimensional channels and a three-dimensional periodic array of spheres [21].

By default, in both mentioned approaches, increasing the Reynolds number means either a decrease in viscosity, an increase in the flow velocity, or an increase in the size of the system. The first approach is limited by the minimum theoretical viscosity due to stability issues. Moreover, it becomes inaccessible in the $\tau=1$ formulation [18], where, by definition, we restrict the model to a fixed relaxation time, which involves fixing viscosity. Thus, considering the maximum velocity limit due to Mach number restrictions, increasing the Reynolds number by increasing the size of the system is the only way to obtain stable, high-Reynolds number fluid flow solutions in both MLBM and LBM1.

The aim of this work is to investigate the procedure of the domain reference length scaling in the meshless LBM (MLBM) to change the Reynolds number at zero computational and algorithmic cost [22] and to implement the meshless streaming step in the LBM1 model to achieve higher Reynolds numbers without additional memory overhead. We apply this scaling method to classic fluid flow problems: Kármán vortex street, lid-driven cavity and in the inertial regime of a flow through a porous medium.

2 Methods

2.1 Lattice Boltzmann Method

Lattice Boltzmann Method [15, 16] is a numerical tool for solving the discrete Boltzmann equation. During each timestep, two steps run - the streaming and the collision of the velocity distributions. The streaming is defined as:

$$f_{k}(t+1,\boldsymbol{x})=f_{k}^{\text{post}}(t,\boldsymbol{x}+\boldsymbol{e}_{k^{\prime}})$$ (1)

where $f_{k}$ is the $k$-th distribution function with its streaming vector $\boldsymbol{e}_{k}$, $k^{\prime}$ denotes the vector antiparallel to $\boldsymbol{e}_{k}$ ($\boldsymbol{e}_{k}=-\boldsymbol{e}_{k^{\prime}}$) and the superscript ’post’ denotes the post-collision distribution function. According to this notation, $\boldsymbol{x}+\boldsymbol{e}_{k^{\prime}}$ denotes the departure node of the $k$-th population being streamed to $\boldsymbol{x}$. Eq. (1) is written in a non-dimensionalized form with the timestep length and the lattice spacing equal to 1.

In this work, we use D2Q9 BGK model with 9 streaming directions ($k=0,1,...,8$):

$$\boldsymbol{e}_{k}\in\left((0,0),(1,0),(0,1),(-1,0),(0,-1),(1,1),(-1,1),(-1,-1),(1,-1)\right).$$ (2)

Due to the Lagrangian nature of the streaming step, such a choice of the streaming vectors imposes the use of a square lattice (refer to Fig. 1 for a graphical representation of the LBM lattice and streaming directions). The collision step calculates $f^{\text{post}}_{k}$ from the current distributions and macroscopic fields. This work uses two relaxation time collision term [23]:

$$f^{\text{post}}_{k}=f_{k}(t,\boldsymbol{x})-\frac{1}{\tau^{+}}\left[f_{k}^{+}(t,\boldsymbol{x})-f^{\text{eq}+}_{k}(t,\boldsymbol{x})\right]-\frac{1}{\tau^{-}}\left[f_{k}^{-}(t,\boldsymbol{x})-f^{\text{eq}-}_{k}(t,\boldsymbol{x})\right]$$ (3)

where $\tau^{\pm}$, are relaxation times for the symmetric and anti-symmetric parts of the distributions, which are defined as:

$$\begin{array}[]{c}\begin{aligned} f_{k}^{+}=\frac{f_{k}+f_{k^{\prime}}}{2},&\quad&f_{k}^{-}=\frac{f_{k}-f_{k^{\prime}}}{2}\\[2.15277pt] f_{k}^{\text{eq}+}=\frac{f_{k}^{\text{eq}}+f_{k^{\prime}}^{\text{eq}}}{2},&\quad&f_{k}^{\text{eq}-}=\frac{f_{k}^{\text{eq}}-f_{k^{\prime}}^{\text{eq}}}{2}.\\ \end{aligned}\end{array}$$ (4)

The relaxation times are related to each other by:

$$\Lambda=\left(\tau^{+}-\frac{1}{2}\right)\left(\tau^{-}-\frac{1}{2}\right).$$ (5)

For the non-fixed relaxation time LBM and MLBM we use the value of $\Lambda=3/16$ [7]. The kinematic viscosity in the lattice units is then:

$$\nu_{lb}=c_{s}^{2}\left(\tau^{+}-\frac{1}{2}\right)$$ (6)

where $c_{s}=1/\sqrt{3}$ denotes the lattice speed of sound. The equilibrium distributions $f_{k}^{\text{eq}}$ are expressed as:

$$f^{\text{eq}}_{k}(t,\boldsymbol{x})=\omega_{k}\rho_{lb}\left[1+\frac{\boldsymbol{e}_{k}\cdot\boldsymbol{u}_{lb}}{c_{s}^{2}}+\frac{(\boldsymbol{e}_{k}\cdot\boldsymbol{u}_{lb})^{2}}{2c_{s}^{4}}-\frac{\boldsymbol{u}_{lb}^{2}}{2c_{s}^{2}}\right]$$ (7)

where $\omega_{k}$ is the weight specific to the $k$-th streaming direction:

$$\omega=\left(\frac{4}{9},\frac{1}{9},\frac{1}{9},\frac{1}{9},\frac{1}{9},\frac{1}{36},\frac{1}{36},\frac{1}{36},\frac{1}{36}\right).$$ (8)

The macroscopic density and velocity in lattice units, $\rho_{lb}\!=\!\rho_{lb}(t,\boldsymbol{x})$ and $\boldsymbol{u}_{lb}\!=\!\boldsymbol{u}_{lb}(t,\boldsymbol{x})$ respectively, at time $t$ and point $\boldsymbol{x}$, are obtained from the discrete populations:

$$\rho_{lb}= \sum\limits_{k=1}^{q}f_{k}$$ (9)

$$\boldsymbol{u}_{lb}= \frac{1}{\rho_{lb}}\sum\limits_{k=1}^{q}f_{k}\boldsymbol{e}_{k}$$

The pressure is related to the density via the lattice speed of sound: $p_{lb}=\rho_{lb}c_{s}^{2}$.

In the numerical implementations, at each time step, Eq.  (3) is first calculated to obtain the values of the post-collision distribution function (collision step). Then,  (1) advects post-collision distributions to neighboring nodes (streaming step). Because lattice nodes $\boldsymbol{x}$ coincide with the departure/arrival nodes of the streaming step, transport is purely Lagrangian and amounts to an index shift in the distribution function array.

The fixed-relaxation time variant of LBM assumes $\tau^{+}=1$ and $\Lambda=1/4$ [18, 19], which by Eq. (5) sets $\tau^{-}=1$. In such case, Eq.(4) simplifies the streaming and collision procedures from Eqs. (1) and (3) to the same form as if Bhatnagar-Gross-Krook collision [24] was used with the unit relaxation time:

$$f_{k}(t+1,\boldsymbol{x})=f^{\text{eq}}_{k}(t,\boldsymbol{x}+\boldsymbol{e}_{k^{\prime}}).$$ (10)

This substitution alleviates the need to store the distribution function explicitly. In each timestep, the equilibrium distributions are calculated directly from the current macroscopic fields according to Eq. (7) and summed as in Eq. (9) to obtain the macroscopic fields in the next time step.

2.2 Formulation of meshless LBM

The meshless LBM algorithm considered in this work belongs to the family of off-grid LBM [25, 26, 27, 28, 20]. Those methods trade the Lagrangian approach to solving the streaming step for the semi-Lagrangian or Eulerian method, such that the space and velocity discretizations get decoupled, and nodes no longer need to be arranged in a regular grid. In the implementations taking the semi-Lagrangian approach to solve the streaming step, like the one discussed here, the discretization points (so-called Eulerian nodes) need not coincide with the departure points of the distribution function (so-called Lagrangian nodes). Instead, with each Eulerian point $\boldsymbol{x}_{i}$ we identify $q$ Lagrangian points to which the corresponding distributions are interpolated and from which they are advected during the streaming step (see Fig. 1). To relate the Eulerian and Lagrangian nodes’ positions to one another, in MLBM the streaming vectors need to be expressed in physical units: $\delta\boldsymbol{x}_{k}=\boldsymbol{e}_{k}\delta x$ where $\delta x$ is the streaming distance. In this manner, the positions of the departure nodes of the $k$-th population are related to the positions of their Eulerian nodes as $\boldsymbol{x}+\delta\boldsymbol{x}_{k^{\prime}}$. We introduce the interpolation step between the collision and the streaming to approximate the post-collision distributions in the appropriate Lagrangian points. Once the approximation is done, the streaming and the collision step are performed in the same way as in LBM (see. Fig. 2). We note that in the models with zero-velocity population discretized (i.e. $\boldsymbol{e}_{0}$ in Eq. (2)), one does not perform the streaming for this population explicitly.

Among many interpolation methods available to solve the streaming step on irregular discretizations, meshless methods [29, 30, 31] are less computationally demanding and error-prone in the discretization process, they achieve boundary-compliant discretizations, offer high order approximations as well as feasible local refinement and adaptivity. Such approaches have been frequently used to solve transport partial differential equations [32, 33, 34, 35, 36, 37].

Figure 1: Comparison of space discretizations in the standard (left) and meshless (right) LBM. The standard LBM requires a square grid; the meshless LBM can operate on scattered, boundary-compliant node sets. Note that in MLBM, velocity is still discretized using a lattice, resulting in several departure points (so-called Lagrangian nodes, filled squares) assigned to each discretization point (so-called Eulerian nodes, filled circles).

Figure 2: The procedure performed at a single timestep of the meshless LBM algorithm. Circles represent Eulerian points; squares denote Lagrangian points. A dashed loop encloses the stencil of a Lagrangian node, which consists of the nine closest neighbors of the Lagrangian node (blue circles). The interpolated distribution function is streamed to the Eulerian point at the center of the presented square lattice.

The semi-Lagrangian streaming step in the discussed MLBM uses meshless interpolation based on radial basis functions-generated finite differences (RBF-FD)[29]. This allows one to approximate the value of a linear operator $\mathcal{L}$ applied to the underlying function in an arbitrary point as a linear combination of the known function values stored in the neighboring nodes:

$$\mathcal{L}(f)|_{\boldsymbol{x}_{0}}\approx\sum\limits_{i=1}^{N_{L}}f(\boldsymbol{x}_{i})w_{i}$$ (11)

where $\boldsymbol{x}_{0}$ denotes the approximation point, $\boldsymbol{x}_{i}$ denote the $N_{L}$ approximation stencil members and $w_{i}$ are interpolation weights. In particular one can speak of the unit operator $\mathcal{L}$ ($\mathcal{L}(f)\equiv f$) which results in the interpolation formula in Eq. (11). We use local approximation supports for each Lagrangian node consisting of its $N_{L}=25$ closest Eulerian neighbors, unless otherwise stated. To achieve high-order approximation we use polynomial augmentation [38] of order two, unless otherwise stated. For a detailed description of the algorithm, we refer the reader to our previous work [20].

The application of the above formalism to LBM1 is the following. The macroscopic density and velocities are interpolated to each of the Lagrangian points. Then, the equilibrium distributions are calculated in the Lagrangian nodes, Eq.(7), and their values are summed according to Eq.(9) to obtain the next-timestep macroscopic fields in the Eulerian nodes. The interpolation of the macroscopic fields, rather than the equilibrium distributions, frees one from the need to store $f_{k}^{\text{eq}}$ in the memory.

2.3 Increasing Reynolds number in meshless LBM

To assess the ratio between inertial and viscous forces in a flow, the dimensionless Reynolds number is used:

$$Re=\frac{UL}{\nu}$$ (12)

arising naturally from non-dimensionalization of incompressible Navier-Stokes momentum equation:

$$\begin{array}[]{l}\dfrac{\partial\boldsymbol{u}}{\partial t}+(\boldsymbol{u}\cdot\nabla)\boldsymbol{u}=-\dfrac{\nabla p}{\rho}+\boldsymbol{g}+\nu\nabla^{2}\boldsymbol{u}\\ \dfrac{\partial\boldsymbol{u^{*}}}{\partial t^{*}}+(\boldsymbol{u^{*}}\cdot\nabla^{*})\boldsymbol{u^{*}}-\nabla^{*}p+\boldsymbol{g}^{*}+\dfrac{1}{Re}\nabla^{*2}\boldsymbol{u^{*}}\\ \text{where }\boldsymbol{x}^{*}=\boldsymbol{x}/L,\quad\boldsymbol{u}^{*}=\boldsymbol{u}/U,\quad t^{*}=t\dfrac{U}{L}=tT,\quad p^{*}=p/\rho U^{2},\quad\boldsymbol{g}^{*}=\boldsymbol{g}\dfrac{U^{2}}{L}\\ \end{array}$$ (13)

where asterisk denotes non-dimensionalized quantities, $\nu$ is the fluid’s kinematic viscosity, $\boldsymbol{g}$ is the acceleration, and $L$, $U$, $T$ denotes the reference length, velocity, and time, respectively. Conventionally, when inertial forces dominate the system, the Reynolds number is larger than unity $Re>1$. The practical utility of the Reynolds number lies in its ability to relate flows of different length and time scales to each other, e.g., flow past a real object versus its scaled model used in an experiment. The dependence of hydrodynamic coefficients, like drag, lift, or permeability, is often presented in relation to $Re$, thus making the results more universal. Reynolds number also allows to estimate the relative importance of non-linearity, turbulence, and momentum diffusion phenomena in a flow when choosing an appropriate physical model to describe the system.

In MLBM, the streaming distance length $\delta x$ can be varied independently from the positions of the nodes discretizing the domain. It means that the whole space spanned by the streaming distance $\delta x$ and the Eulerian discretization parameter $h$ is available in MLBM, in contrast to traveling along the $\delta x=h$ line in the standard LBM (see Fig. 3). For instance, consider a flow domain of size $L=100$. Let us assume that the streaming distance is $\delta x=0.1$ (we do this by specifying the positions of Lagrangian points to differ from their Eulerian lattice centers by $\delta\boldsymbol{x}_{k}$ vectors). With such choice of $L$ and $\delta x$ we arrive at the domain size in lattice units equal to $L_{lb}=L/\delta x=1000$, no matter what Eulerian discretization parameter $h$ we choose. Now, increasing the size of the domain $L$, e.g., by multiplying the position vectors of the Eulerian points $\boldsymbol{x}_{i}$ by some constant $\alpha_{E}$ (with $\delta x$ unaltered) changes the domain size in lattice units $\alpha_{E}$ times. The same effect can be achieved by multiplying the streaming distance by $\alpha_{\delta}=1/\alpha_{E}$ and keeping the positions of Eulerian points constant (we use the notation $\alpha_{\delta}$ to distinguish the factor for the scaling of the streaming distance from the one for the scaling of the Eulerian points positions). This is because the approximation coefficients $w_{i}$ from Eq. (11) are invariant under translation, rotation, and scaling of the stencil nodes and the interpolation point, and the presented MLBM algorithm uses the non-dimensionalized formulation of LBM ($\delta x=1$ and $\delta t=1$ in lattice units).

Both strategies (depending on scaling the domain size in physical units or the streaming distance length) alter the flow’s Reynolds number in the following ways. Considering the Reynolds number calculated from the quantities in lattice units, both strategies change the size of the domain, changing the $L_{lb}$ value from the previous paragraph. The same happens with $Re$ calculated using the quantities in physical units when the $\alpha_{E}$ scaling is used. Finally, when the scaling is applied to the streaming distance (i.e. $\delta x\rightarrow\alpha_{\delta}\delta x$), to keep the physical viscosity unchanged, timestep length must scale with the square of $\alpha_{\delta}$, since $\nu=c_{s}^{2}(\tau-1/2)\delta x^{2}/\delta t$. When this happens, the physical velocity $U$ changes into its counterpart in the scaled setup $U^{\delta}$ as

$$U=U_{lb}\frac{\delta x}{\delta t}\>\rightarrow\>U^{\delta}=U_{lb}\frac{\delta x\alpha_{\delta}}{\delta t\alpha_{\delta}^{2}}=U\>\frac{1}{\alpha_{\delta}}$$ (14)

which in the case of decreasing the streaming distance ($\alpha_{\delta}<1$) means increasing the velocity $U$ in physical units.

Figure 3: The $\delta x-h$ discretization parameters space and its three discussed directions (denoted by gray dashed arrows). Blue circles denote Eulerian points, and small dark gray squares denote Lagrangian nodes. Moving along each of the shown directions results in different scaling of the streaming length $\delta x$ and Eulerian discretization parameter $h$. In particular, in LBM, moving only along the constant $\delta x/h$ direction is possible.

In general, the scaling of the streaming distance is followed by the scaling of the conversion factors from lattice units to physical units:

$$\begin{array}[]{lcr}C_{L}=\delta x&\rightarrow&C_{L}^{\delta}=\delta x\alpha_{\delta}=C_{L}\alpha_{\delta}\\ C_{t}=\delta t&\rightarrow&C_{t}^{\delta}=\delta t\alpha_{\delta}^{2}=C_{t}\alpha_{\delta}^{2}\\ C_{u}=\delta x/\delta t&\rightarrow&C_{u}^{\delta}=\delta x\alpha_{\delta}/\delta t\alpha_{\delta}^{2}=C_{u}/\alpha_{\delta}\\ C_{\nu}=\delta x^{2}/\delta t&\rightarrow&C_{\nu}^{\delta}=C_{\nu}\\ C_{g}=\delta x/\delta t^{2}&\rightarrow&C_{g}^{\delta}=\delta x\alpha_{\delta}/\delta t^{2}\alpha_{\delta}^{4}=C_{g}/\alpha_{\delta}^{3}\end{array}$$ (15)

We note that the discussed method of increasing the Reynolds number in MLBM should apply when Eulerian schemes are used for solving the streaming step[39]. The scaled streaming distance would then be equivalent to the scaled streaming velocity in the term $\boldsymbol{c}\cdot\nabla$ of the Boltzmann equation, where $\boldsymbol{c}$ is the microscopic population velocity.

3 Results

3.1 Validation and execution time analysis

We investigate the flow over a cylindrical obstacle in a plane channel, Fig. 4. The length of the channel is $3$, the height of the channel is $1$ and the cylindrical obstacle of diameter $d=0.14$ is placed on the channel’s axis, at a distance $0.5$ from its inlet.

For MLBM simulations we use the scattered discretization near the cylinder and in its wake and the regular discretization elsewhere. We use the meshless RBF-FD approximation for the Lagrangian points whose closest Eulerian neighbors do not form a regular lattice. In the remaining points we use biquadratic Lagrange interpolation. The distance between the nodes vary from $h_{\text{min}}=1/200$ on the cylinder’s surface to $h_{\text{max}}=1/25$ away from the cylinder. The set of points used for the discretization is symmetric about the channel’s axis to achieve a symmetric lift force profile in time. The boundary conditions on the obstacle’s surface is the no-slip wall. It is implemented using a simple bounceback rule for the populations for which the Lagrangian points lie inside the cylinder. The inlet and the outlet has the prescribed parabolic velocity profile with the maximum value of $U_{0,lb}=10^{-2}$ and unit density $\rho=1$. Zero macroscopic velocity $\boldsymbol{u}=\boldsymbol{0}$ and unit density $\rho=1$ is imposed on the top and bottom boundaries. The former four boundaries are implemented by overwriting all populations in their nodes with the equilibrium distributions parametrized with the desired macroscopic density and velocity. The initial condition are the equilibrium populations parametrized with the unit density $\rho=1$ and zero macroscopic velocity $\boldsymbol{u}=\boldsymbol{0}$.

Figure 4: The hybrid regular-irregular discretization of the flow over a cylindrical obstacle used in the simulations (left) and a zoom in at the discretization around the obstacle (right).

The Reynolds number of each case was calculated based on the inlet velocity $U_{0}$ and the cylinder’s diameter $d$:

$$Re=\frac{U_{0}d}{\nu}.$$ (16)

The values considered in MLBM simulations are $Re=105,126,147,168$. We increase $Re$ by decreasing the streaming distance length from $\delta x=1.33\cdot 10^{-3}$ for $Re=105$ to $\delta x=8.33\cdot 10^{-4}$ for $Re=168$. The relaxation time is $\tau^{+}=0.53$ and the diffusive scaling of the timestep is used.

For the reference in timings analysis we perform a series of standard LBM simulations with the same initial and boundary conditions as in the MLBM setup. The Reynolds numbers considered in LBM simulations are $Re=84,105,126,147,168$. The lattice parameter for $Re=84$ is $\delta x=1/200$ and the increase in $Re$ is achieved by the appropriate decrease in $\delta x$, such that $Re=168$ yields the value $\delta x=1/400$. The relaxation time is $\tau^{+}=0.51$ which is close to the stability limit. Also here we use the diffusive scaling of the timestep. We note that MLBM setups use higher value of $\tau^{+}$ due to stability issues.

We investigate the onset of inertial effects in the flow qualitatively by observing the appearance of the Kármán vortex street in the obstacle’s wake and quantitatively via the increase of the Strouhal number $St$:

$$St=\frac{fd}{U_{0}},$$ (17)

where $f$ is the frequency of the wake vortices shedding.

Figure 5: Streamlines of the velocity fields of the flow over a cylindrical obstacle for Reynolds numbers 84 (LBM results) and 168 (MLBM result).

Fig. 5 shows velocity streamlines of the channel flow for two extreme Reynolds numbers: $Re=84$ obtained with LBM and $Re=168$ obtained with MLBM. The onset of inertial effects is visible as the distance between the repeating velocity field patterns in the obstacle’s wake decreases.

Figure 6: Strouhal number values achieved in flows around a cylindrical obstacle as a function of the flow’s Reynolds number. The dashed part of the reference line for Norberg et al. [40] indicates the region where irregular shedding took place.

Fig. 6 shows the values of the Strouhal number $St$ against the Reynolds number of the flows. The present results are in good agreement with a fit to the experimental data for the cylinder’s diameter-to-length ratio equal to 2000 [40] and earlier experimental works [41, 42, 43].

Let us now move to the timings analysis. For this comparison, LBM and MLBM codes are executed single-threaded on 2x Intel Xeon E5520 with 12GB RAM machine. Fig. 7 shows the CPU time needed to simulate 1s of vortex shedding. For both methods, the execution time increases with $Re$. This is caused by the diffusive scaling of the timestep with the streaming distance $\delta x$ and the decreasing value of $\delta x$ with the increasing $Re$. This means that more timesteps is needed to reach 1s in the simulation for larger $Re$. Nevertheless, for MLBM the execution time rises slower than for LBM ($\mathcal{O}(Re^{2})$ vs. $\mathcal{O}(Re^{4})$), since in the former it is only the timestep number that increases with $Re$. For LBM, on the other hand, not only the timestep number increases with $Re$, but also the total number of nodes in the domain does. For the highest value of $Re$ considered here, MLBM performs $\sim 3$ times faster than LBM.

Figure 7: CPU time needed to execute 1s of the vortex shedding for LBM and MLBM in a single threaded run for various Reynolds numbers. Both axes use logarithmic scale. The dashed lines denote the 2nd and the 4th order slopes.

To analyze the sources of MLBM speedup, we consider the CPU time needed to perform $10^{3}$ iterations in a single-threaded run on the same machine as before, shown in Fig. 8. The data are presented for LBM and for three scenarios of MLBM. The first scenario uses fully scattered discretization with RBF-FD approximation at each Lagrangian node, with uniform in space internodal distance equal to $1/200$. In the second scenario we exchange the scattered discretization away from the cylinder and its wake for a regular points arrangement. We also use biquadratic Lagrange interpolation for the Lagrangian nodes whose closest Eulerian neighbors form a regular lattice. The third scenario additionally features coarser discretization away from the cylinder and is the one used in the calculation of the Strouhal number. The number of nodes in the LBM and MLBM setups is summarized in Table 1. The presented execution time increases with $Re$ for LBM because of the increasing number of nodes in the domain when the streaming distance $\delta x$ is decreased. This is in contrast to MLBM, when the increase of $Re$ is achieved only by decreasing the value of $\delta x$, with no additional nodes inserted in the domain. Due to this, the execution time stays constant for each scenario across the considered Reynolds numbers. Further, the execution times differ between the MLBM scenarios. About $40\%$ speedup between scenario 1 and 2 is achieved by switching from RBF-FD approximation with $N_{L}=25$ nodes to the computationally cheaper biquadratic Lagrange interpolation away from the cylinder. It is achieved in spite of the increase of the total number of nodes by about $10\%$. Decreasing the discretization density away from the cylinder accounts for the subsequent $86\%$ speedup between scenario 2 and 3.

Figure 8: CPU time needed to execute $10^{3}$ timesteps of LBM and MLBM in a single threaded run for various Reynolds numbers. The visualizations of the discretizations of the three MLBM scenarios are shown on the right hand side. Note that two uniform visualizations show the discretizations coarser than $h=1/200$ for clarity.

Table 1: Number of nodes $N$ discretizing the domain in the LBM and MLBM simulations of the Kármán vortex street.

LBM case	$N$
$Re=84$	$119988$
$Re=105$	$187290$
$Re=126$	$269522$
$Re=147$	$366668$
$Re=168$	$478740$

MLBM case	$N$
scattered uniform	$104825$
hybrid uniform	$115185$
hybrid coarsened	$10117$

3.2 Application of MLBM1 in the lid-driven cavity

The flow domain is a $[0,1]\times[0,1]$ square (Fig. 9). We use irregular discretizations with the Eulerian discretization parameter $h=1/160$ uniform in space. We set the top lid velocity to $u_{lb}(x,1)=U_{0,{lb}}(1-(2x^{*}-1)^{24})$ where $U_{0,{lb}}=0.1$ is the reference velocity in LB units and $x^{*}=\min(\max(0,x),1)$ is the truncated $x$-coordinate. The no-slip boundary condition on the left, right, and bottom walls is implemented via multireflection bounceback [44]. The unknown populations at the top boundary are equilibrium distributions parametrized with $\rho_{lb}=1$ and velocities: $u_{lb}(x,1)$ and $v_{lb}(x,1)=0$. We introduce $x^{*}$ to be able to calculate $f^{\text{eq}}_{k}$ streamed from $x<0$ or $x>1$. To increase the Reynolds number, we scale the positions of Eulerian points by the factor of $\alpha_{E}=1,4,10,32,50$. Along with the streaming distance $\delta x=6\cdot 10^{-3}$ and physical viscosity $\nu=1$, this gives Reynolds numbers based on the square side length and top lid velocity equal to $100,400,1000,3200$ and $5000$. We initialize the simulation with equilibrium distributions parametrized with zero velocity and unit density and iterate the simulation until the relative change of the $x$-component of the velocity between two timesteps defined as:

$$\Delta u=\frac{\max\limits_{i}\left(\left|u_{i}^{n+1}-u_{i}^{n+1}\right|\right)}{U_{0}}$$ (18)

falls below $10^{-12}$ or up to $6\cdot 10^{6}$ (for $\alpha_{E}=32$) or $12\cdot 10^{6}$ (for $\alpha_{E}=50$) iterations. By $u_{i}^{n}$ we mean the value of the $x$-component of the velocity at the $i$-th Eulerian node at the $n$-th timestep.

Figure 9: The irregular discretization used in the simulations of the flow in a driven cavity.

Figure 10: Profiles of $u$ and $v$ components of velocity along the lines $x=0.5$ and $y=0.5$ for the driven cavity flow. Small symbols denote MLBM1 results, larger squares denote the reference solution by Ghia [45]. Colors indicate the Reynolds number.

Figure 11: Positions of the centers of the vortices in the driven cavity. Left: the primary vortex, center: the secondary bottom left vortex and right: the secondary bottom right vortex. The dotted line denotes the sequence of increasing $Re$: $100,400,1000,3200$. The MLBM1 data are compared with the results of Ghia [45].

Figure 12: Paths traversed by the centers of the primary vortex (left) and the secondary bottom left vortex (right), $Re=5000$. The MLBM1 data are compared with the vortices centers reported by Ghia [45].

Figure 10 shows the normalized $u$ and $v$ component of velocity profiles along the centerlines of the cavity ($x=0.5$, $y=0.5$ after rescaling the domain by $1/\alpha_{E}$). Small circles denote MLBM1 results and are compared with the results of Ghia [45] denoted by larger squares. A good compliance between the two methods is visible, especially in the middle of the primary vortex. The difference in the velocity profiles increases closer to the boundaries, possibly due to large velocity gradients and shear stresses. Fig. 11 shows the centers of the primary vortex and the secondary vortices in the bottom left and bottom right corners of the cavity. The compliance between MLBM1 results and Ghia’s solution [45] is within about 1% of the domain side length. As the $Re=5000$ case exhibited an oscillatory behavior of slowly decreasing amplitude, we plot the path of the primary vortex and the secondary vortex in the bottom left corner of the cavity in Fig. 12. After $12\cdot 10^{6}$ iterations, corresponding to $t=72$, the centers of the vortices oscillate around points close to those reported by Ghia [45]. We note that to obtain a stable solution in $Re=5000$ case, we needed to increase the order of augmentation in the meshless approximation to the fourth and the number of stencil members to $N_{L}=70$. Despite the stability, spurious shockwaves propagating from the right boundary could be seen throughout the simulation.

3.3 Inertial effects in a flow through a porous sample

Figure 13: The point cloud used for the simulations of a flow through a porous sample: the whole discretized domain and a zoom at the discretization detail (right).

We investigate the onset of inertial effects in a flow through a porous medium (see Fig. 13) of porosity $\phi=0.64$ with MLBM1. We use irregular discretization with the Eulerian discretization parameter $h=5\cdot 10^{-3}$. The no-slip boundary condition is imposed on the top and bottom boundary and the circular obstacles’ surfaces using the multireflection bounce back rule [44]. The periodic boundary condition is imposed on the left and right boundary by periodic search of the stencil members. The flow is forced by a constant acceleration of value $g_{lb}=10^{-7}$ along the $x$ direction. The meshless discretization of the domain corresponds to its size of $1\times 1$. The investigated values of the scaling parameter $\alpha_{E}$ range from $1$ to $17$. We note that for $\alpha_{E}=18$, we could not obtain a stable solution despite increasing the stencil size to $N_{L}=70$ and using fourth-order polynomial augmentation. The streaming distance length is chosen to be $\delta x=2.5\cdot 10^{-3}$, timestep length is $\delta t=1.04\cdot 10^{-6}$ and the kinematic viscosity is $\nu=1$. The Reynolds number for each case was calculated based on the mean flow velocity in the $x$-direction $\langle u\rangle$:

$$\langle u\rangle=\frac{1}{\phi}\int\limits_{D}d\Omega\>u(\boldsymbol{x}),$$ (19)

sample’s side length equal to $\alpha_{E}$ and kinematic viscosity $\nu$.

We show the visualization of the velocity field for three chosen values of $\alpha_{E}$ in Fig. 14. As the Reynolds number increases, recirculation zones grow causing the streamlines to separate from the grains’ wakes.

Figure 14: Streamlines of the velocity field of flows through the porous domain scaled with the factors $\alpha_{E}=1,8,17$ ($Re=0.76$, $305$, and $1511$, respectively). The points coordinates were rescaled by a factor of $1/\alpha_{E}$ prior to plotting.

The inertial effects appearance can be quantified in terms of kinetic energy distribution in the system [5]. The left plot in Fig. 15 shows the fraction of the kinetic energy defined as $e_{k}(\boldsymbol{x})=\boldsymbol{u}^{2}(\boldsymbol{x})$ confined in vortices to the total kinetic energy in the system:

$$\frac{\displaystyle\int\limits_{\boldsymbol{x}\>:\>Q<0}e_{k}(\boldsymbol{x})\>d\Omega}{\displaystyle\int e_{k}(\boldsymbol{x})\>d\Omega}\>.$$ (20)

We use the Q-criterion to define a vortex as a volume where:

$$Q=\frac{1}{2}\left(|S|^{2}-|\Omega|^{2}\right)<0,\text{ where }S=\frac{1}{2}\left(\frac{\partial u_{i}}{\partial x_{j}}+\frac{\partial u_{j}}{\partial x_{i}}\right),\>\Omega=\frac{1}{2}\left(\frac{\partial u_{i}}{\partial x_{j}}-\frac{\partial u_{j}}{\partial x_{i}}\right)$$ (21)

and $|A|^{2}=\sum\limits_{i,j}A_{ij}^{2}$. By $u_{i}$ we mean the velocity components (i.e. $u_{0}\equiv u$ and $u_{1}\equiv v$). In the range of $Re=100$–$700$ a rapid increase in the value of the fraction from Eq. (20) occurs.

We can relate this to the transition in the behavior of the $\langle u\rangle$. In the Darcy regime (inertial effects negligible), $\langle u\rangle$ is proportional to the pressure loss between the inlet and the outlet of the sample:

$$\langle u\rangle=-k\frac{\Delta P}{L\mu}=-k\frac{g}{\nu}$$ (22)

where $k$ is the sample’s permeability, $L$ is the sample length, and the last equality holds if we relate the acceleration $g$ forcing the flow to the mean pressure drop via $\rho g=\Delta P/L$. Next, Eq. (22) needs to be related to the changes in the size of the sample since this is the way of increasing $Re$ in our study. If in Eq. (13) we omit the pressure gradient term and change the reference length from $L$ to $\alpha_{E}L$, we find that the only difference between the two systems (the one rescaled by $\alpha_{E}$ and the non-rescaled one) is in the term $Re=UL/\nu=L^{2}/T\nu$. This means that a system scaled by a factor of $\alpha_{E}$ behaves as if it was not scaled, but the viscosity was changed by a factor of $\alpha_{E}^{-2}$ (due to the term $L^{2}$ in the denominator of $Re$). Inserting this relation into Eq. (22) gives:

$$\langle u\rangle\propto k\frac{g}{\nu\alpha_{E}^{-2}}\propto\alpha_{E}^{2}$$ (23)

A similar procedure applied in the Forchheimer regime (leaving only the second order term $\langle u\rangle^{2}$) yields:

$$\langle u\rangle^{2}\propto\alpha_{E}^{2}.$$ (24)

Presenting Eqs. (23) and (24) in logarithmic scale gives:

$$\log\langle u\rangle\propto\begin{cases}2\log\alpha_{E},&\text{Darcy regime}\\ \log\alpha_{E},&\text{Forchheimer regime}\\ \end{cases}$$ (25)

which is visible in the right plot of Fig. 15. The transition between the second and the first order relations, denoted by dashed grey lines in Fig. 15, occurs at the scaling factor $\alpha_{E}\approx 7$ for which $Re\approx 25$ (see the inset in the left subplot of Fig. 15) and corresponds well with the predictions of the onset of inertial effects from the vortices’ kinetic energy approach.

Figure 15: Left: the fraction of the kinetic energy confined in the $Q<0$ volumes as a function of the Reynolds number. The inset shows the relation between the Reynolds number and the scaling factor $\alpha_{E}$. Right: mean $x$-component velocity defined by Eq.(19) as a function of the scaling factor. The dashed lines denote the 1st and the 2nd order slopes.

Computational complexity and memory demands of the MLBM1 algorithm

Matyka and Dzikowski [18] derived the expression for the memory demands of their algorithm, which yields:

$$M_{\text{LBM1}}=N_{\text{LBM1}}\cdot 2N_{f}c$$ (26)

where $N_{\text{LBM1}}$ is the number of nodes in the discretization, $N_{f}$ is the number of macroscopic fields (in our case 3 - $\rho$,$u$ and $v$), and $c$ is the number of bytes per single number (e.g. 8 for double). MLBM1 uses additional memory to store the interpolation weight vectors for each of the $N(q-1)$ Lagrangian nodes (the term $-1$ appears because no interpolation is needed for the populations with zero microscopic velocity). Assuming that each interpolation stencil has the same number of nodes $N_{L}$, this results in the memory demands for MLBM1:

$$M_{\text{MLBM1}}=N_{\text{MLBM1}}\cdot(2N_{f}+N_{L}(q-1))c$$ (27)

where $N_{\text{MLBM1}}$ is the number of Eulerian nodes. Taking the ratio of the two gives:

$$r_{M}=\frac{M_{\text{MLBM1}}}{M_{\text{LBM1}}}=\frac{N_{\text{MLBM1}}}{N_{\text{LBM1}}}\cdot\frac{2N_{f}+N_{L}(q-1)}{2N_{f}}=\frac{N_{\text{MLBM1}}}{N_{\text{LBM1}}}\cdot\left(1+\frac{N_{L}(q-1)}{2N_{f}}\right)$$ (28)

In our case, $N_{f}=3$, $q=9$, and for most of the presented cases, $N_{L}=25$. Considering that sparser discretizations can be used in the meshless algorithm, it turns out that MLBM1 can outperform LBM1 regarding memory demands for the same streaming distance value. For example, in the driven cavity flow, we used $\delta x=6\cdot 10^{-3}$, which would give the number of nodes in the LBM1 setup ranging from $2.8\cdot 10^{4}$ for $\alpha_{E}=1$ to $7\cdot 10^{7}$ for $\alpha_{E}=50$. At the same time the term in brackets in Eq. (28) is equal to $34$ for $\alpha_{E}<50$ and $94$ for $\alpha_{E}=50$ (due to the change of $N_{L}$ for $\alpha_{E}=50$). This means that for the MLBM1 domain with $N=22615$ Eulerian nodes, the ratio $r_{M}$ ranges from $27$ times more memory in the case of MLBM1 for $\alpha_{E}=1$ to $32$ times more memory in the case of LBM1 for $\alpha_{E}=50$. This shows the possibility of considerably reducing the number of nodes in LBM1 models with an off-grid (especially meshless) streaming procedure. When a proper local refinement strategy is used and less computationally intensive approximation scheme is employed e.g. in the bulk of the flow, as shown in Section 3.1, the memory savings may come at no cost to the solution’s accuracy.

The computational complexity of MLBM1 carries the overhead of interpolation during the streaming step. To exploit the absence of explicit storage of the distribution function in the memory, the collision (i.e., the calculation of $f^{\text{eq}}_{k}$) must take place in Lagrangian nodes, and the macroscopic fields must be interpolated to those points. Otherwise (if the collision took place in Eulerian nodes), $f^{\text{eq}}_{k}$ would have to be stored in memory between its calculation and interpolation. This means that for each Eulerian point, one needs $qO_{\text{col}}+(q-1)N_{f}O_{\text{int}}$ operations per timestep, where $O_{\text{col}}$ is the number of operations needed to calculate $f^{\text{eq}}_{k}$ and $O_{\text{int}}$ is the number of operations needed to interpolate a scalar within a stencil (note that for $k=0$ one does not need to interpolate, but still the calculation of $f^{\text{eq}}_{0}$ is necessary). This shows that the savings in memory for $q$ distributions come at the cost of $N_{f}$ times more interpolations which have to be done for each Lagrangian node, compared to MLBM (in MLBM, one interpolates only one population to each Lagrangian point).

4 Discussion

The decoupled space and velocity discretizations in MLBM change the paradigm of increasing the model’s size typical for LBM and make it more similar to the one of the Navier-Stokes solvers. In LBM, the information about the characteristic size of the problem is encoded in the number of lattice links/nodes that discretize this size. So is the accuracy of the approximation in space. In MLBM, the information about the system’s size is stored in the ratio between the streaming distance and the distance between Eulerian points representing the characteristic length. The approximation accuracy in space depends on the distances between the approximation stencil members. It is invariant under their translation, rotation, and scaling by a constant as long as the approximated function is devoid of discontinuities. Consequently, in MLBM, the velocity discretization parameters are determined by the desired compressibility errors. In contrast, the space discretization parameters are determined by the flow domain’s geometry and the desired approximation accuracy in space. This is interpreted as hypervolumes in the space of those two groups of discretization parameters accessible to MLBM, where each point might be assigned its Reynolds number, in contrast to hypersurfaces in the case of LBM.

We chose the presented way of increasing the Reynolds number in MLBM to make the procedure as similar to the LBM approach as possible and highlight how the two methods relate to each other. This allows us to state the following. First, the qualitative results show that the described approach causes the Reynolds number to increase, in line with the findings of previous authors [22]. It also suggests that the propagation of information in the model based on the approximation, rather than the exact streaming, does not impair the physical mechanism of mass and momentum transfer in the discrete Boltzmann equation. Large streaming distances compared to the spacings between the Eulerian points correspond to fast momentum diffusion, small $\delta x/h$ ratios - on the contrary. Second, quantitative analysis underlines the importance of the approximation accuracy for the solution. At the same time, the presented discretization parameters may serve as a crude estimation of the method’s stability limits along with the aids to mitigate the divergence.

The accuracy and convergence of MLBM and other off-grid LBM methods were investigated in numerous works [27, 46, 20], especially in application to inertial flows in works by Lin and others [47] or Musavi and Ashrafizaadeh [48, 28]. According to the findings by, e.g., Krämer and others [46], the decreasing ratio of $\delta x/h$ introduces larger approximation errors in OLBM. When semi-Lagrangian streaming is used, this has directly to do with the error term $\mathcal{O}\left(h^{p+1}/\delta x\right)$ characteristic for semi-Lagrangian methods, as suggested in our previous work [20]. To prevent the divergence of this term, for the growing $h$, one would need to either increase $\delta x$, which would again decrease the size of the domain, or increase the accuracy of the approximation. We presented this procedure in two tests. First, we used regular space discretization in the far-field subdomain of the flow past a cylinder, in contrast to an irregular discretization around the cylinder. Second, in the driven cavity test for the highest Reynolds number ($Re=5000$), the augmentation order and the number of stencil members were increased compared to the cases with lower $Re$ to mitigate the divergence of the simulation. However, the mentioned flow through a porous medium for $\alpha_{E}=18$ shows that increasing the approximation order or the stencil size does not always guarantee stability retrieval. In such cases, the Eulerian discretization refinement should be used. On the other hand, from the LBM error point of view, decreasing the $\delta x/h$ ratio acts in favor of the compressibility error, which would increase when higher $Re$ would be achieved by larger inlet velocity or acceleration. Finally, the error of the solution obtained with MLBM on rescaled domains may also come from the under-resolution of physical phenomena, such as small vortices typical for inertial flows. However, this under-resolution might be exploited to the method’s advantage, as suggested by Chen [49].

The solution’s accuracy considerations must always be counterbalanced by preserving the simulation’s efficiency and robustness. In MLBM, one can look for memory and computational savings in discretizing space and velocity independently. First, the feasibility of the local static and dynamic refinement in meshless approximation methods offers great promise to reduce the number of Eulerian points in MLBM compared to LBM. Apart from an arbitrary refinement, as shown in the present work in the case of the flow around a cylinder, the choice of a proper error indicator, such as quantities based on vorticity suggested by Fakhari and Lee [50] or on the velocity field divergence [15], can allow for dynamic $h$- or $p$-refinement during the simulation run. In turn, the error of the solution becomes more uniform in space granting a better accuracy at a little increase of the computational expense. Second, switching to regular discretizations with the approximation schemes computationally less intensive than the scattered ones, can considerably decrease the computational burden in the parts of the domain devoid of geometrical intricacies. The research on this approach is given attention in the field of meshless analysis [51, 52, 53]. Finally, LBM1 models [18, 19] reduce the memory overhead by explicitly storing just a few macroscopic fields rather than a large number of velocity distributions. We show that by introducing meshless space discretization, it is possible to achieve a several-fold memory advantage in comparison to the lattice-based LBM1, which is already estimated to use about $76\%$ less memory than its arbitrary relaxation time counterpart [18].

5 Conclusions

We investigate the procedure to increase the Reynolds number of a flow from laminar to the transitional regime in meshless LBM simulations by increasing the system’s domain size. We achieve this solely by multiplying the coordinates of Eulerian points by some constant $\alpha_{E}>1$, or equivalently - the streaming distance by its inverse $1/\alpha_{E}$, thus no additional nodes are inserted into the domain. We apply this method of the Reynolds number increase to a general meshless LBM scheme and to the meshless implementation of a memory-efficient LBM model with the relaxation time $\tau=1$ (LBM1). We show that the procedure works by numerically investigating three test cases in a range of Reynolds numbers - von Kârmn̂ vortex street behind a cylindrical obstacle of infinite length, driven cavity flow, and a flow through a porous sample. The qualitative observation of the velocity fields indicates the onset of inertial effects in the flows. Also, the obtained values of hydrodynamic parameters agree well with the numerical and experimental references in a wide range of Reynolds numbers. We highlight the computational speedup of MLBM compared to the standard LBM when a proper choice of the nodes arrangement and approximation scheme is made. From the performed simulations, we estimate the possible memory savings of the meshless LBM1 compared to the standard (lattice-based) LBM1 to be as high as 32-fold and pinpoint the source of the computational overhead of the meshless LBM1.

Our findings open new perspectives for further studies on the development and applications of MLBM, and its competitiveness with LBM. Questions yet to be answered address the issues of the possibility of simultaneous scaling of space and velocity discretizations, or the detailed mechanism of the information propagation when the non-purely Lagrangian streaming step is used.

6 Acknowledgments

Funded by National Science Centre, Poland under the OPUS call in the Weave programme 2021/43/I/ST3/00228. This research was funded in whole or in part by National Science Centre (2021/43/I/ST3/00228). For the purpose of Open Access, the author has applied a CC-BY public copyright licence to any Author Accepted Manuscript (AAM) version arising from this submission.