Integrating a Stabilized Radial Basis Function Method with Lattice Boltzmann Method

: The lattice Boltzmann method (LBM) has two key steps: collision and streaming. In a conventional LBM, the streaming is exact, where each distribution function is perfectly shifted to the neighbor node on the uniform mesh arrangement. This advantage may curtail the applicability of the method to problems with complex geometries. To overcome this issue, a high-order meshless interpolation-based approach is proposed to handle the streaming step. Owing to its high accuracy, the radial basis function (RBF) is one of the popular methods used for interpolation. In general, RBF-based approaches suffer from some stability issues, where their stability strongly depends on the shape parameter of the RBF. In the current work, a stabilized RBF approach is used to handle the streaming. The stabilized RBF approach has a weak dependency on the shape parameter, which improves the stability of the method and reduces the dependency of the shape parameter. Both the stabilized RBF method and the streaming of the LBM are used for solving three benchmark problems. The results of the stabilized method and the perfect streaming LBM are compared with analytical solutions or published results. Excellent agreements are observed, with a little advantage for the stabilized approach. Additionally, the computational cost is compared, where a marginal difference is observed in the favor of the streaming of the LBM. In conclusion, one could report that the stabilized method is a viable alternative to the streaming of the LBM in handling both simple and complex geometries.


Introduction
In recent years, the lattice Boltzmann method (LBM) has been used as a powerful numerical tool to simulate a wide range of fluid flow and heat transfer problems [1][2][3].The LBM has competitive advantages over conventional computational fluid dynamic (CFD) methods due to its simplicity, parallel computing nature, and the exactness of advection term (no dissipation error) [2,4].The solution in the LBM is obtained by two steps: collision (local) and streaming (non-local).The collision step relaxes the distribution functions to their equilibrium values.In the non-local streaming step, the distribution functions at the local points are streamed (advected) to the neighbor points using the discrete lattice speeds.Notably, this process leads to perfect shifting to the neighbor nodes because the time and space increments are equal.The perfect shifting justifies the exactness of the advection of the distribution functions and explains the high accuracy of the LBM.However, it constrains the applicability of the method to uniform meshes, which impedes the ability of the method to handle complex geometries [4].
The mesh refinement techniques require a second-order interpolation to achieve the correct Navier-Stokes equations.In addition, they need complex interfacing between different levels of refinement, which might adversely impact the conservation across the interface.Furthermore, they require recalculations of the collision for each refinement level.To avoid this, a second approach is used, which uncouples the time and space increments and uses conventional computational fluid dynamic (CFD) methods such as the finite-volume method (FVM) [17][18][19][20][21][22], the finite-difference method (FDM) [23][24][25][26][27][28][29][30][31][32][33][34], the finite-element method (FEM) [35][36][37][38][39][40][41], or the meshless method [42][43][44][45].This approach suffers from numerical diffusion, which leads to dissipation error in proportion to the mesh distribution.Additionally, complex mesh processing is required, thereby increasing the computational cost.In order to remediate this issue, the radial basis function (RBF) method is used because it has a large number of bases that equals the number of the considered neighbors.Lin et al. [46] used the local semi-Lagrangian RBF for the streaming step in LBM.Musavi et al. [47] handled the streaming of the LBM using the local Petrov-Galerkin RBF method, which is based on the weak form of the partial differential equation (PDE).Unfortunately, the methods depending on weak form solve for the average value over a finite element or volume, which have less accuracy than the perfect shifting of the standard LBM.Therefore, a very fine mesh should be used to achieve accuracy close to that of the perfect shifting (streaming).This mesh refinement comes with the price of high computational and storage costs.To reduce the computational cost, the interpolationbased approach (third approach) is used by researchers [48][49][50][51][52].In the third approach, the streaming in the LBM is treated as a Lagrangian, and then the shifted distribution functions are interpolated back to their original location.The interpolation-based approach usually generates numerical diffusion in the solution.This diffusion can be reduced by drastically increasing the number of grids or using a higher-order interpolation.Increasing the number of mesh points is fraught with many negative effects such as raising the computational cost.Although the higher-order interpolation is easy and applicable for the structured grid, it is difficult for non-uniform meshes.
The present work proposes using the interpolation-based approach to achieve a viable alternative to the streaming step, where the high-accuracy RBF interpolation is used to interpolate the post-collision distribution functions.The usage of high-order RBF interpolation enhances the accuracy and allows the usage of coarser mesh when compared to linear interpolation.However, the accuracy of RBF is predicated on the shape parameter, which needs to be selected carefully in order to achieve the desired stability and accuracy.Hence, it is necessary to stabilize RBF for a wide range of shape parameter values.The stability of the RBF depends strongly on shape parameter values that should be carefully selected.Using a novel stabilized RBF method proposed in our previous work [53] could beat this issue.The stabilized RBF approach has a weak dependency on the shape parameter, where any value of the shape parameter can be selected.In the current study, the stabilized approach is used to solve the streaming step of the LBM to enhance the accuracy and applicability of the method to problems with non-uniform grids.To evaluate the proposed method, three benchmark problems are solved using the stabilized RBF method and the streaming of the LBM.Moreover, the computational cost is compared for both used methods.
The rest of the paper is organized as follows.Section 2 describes the LBM.Section 3 explains the RBF interpolation method and how to use it for solving the streaming step of the LBM.Section 4 shows the results of one-and two-dimensional benchmark problems.Finally, Section 5 presents the conclusions.

Lattice Boltzmann Method
The LBM is derived from the Boltzmann transport equation after being discretized in velocity space using a finite velocity set   [54][55][56][57].
where ,    , , , , , , , and  denote the distribution function, the Maxwell-Boltzmann distribution function, the time, the position vector, the relaxation time, the force term, the density, the microscopic velocity vector, and the macroscopic velocity, respectively.Many approaches have been followed in the literature to discretize Equation (1).In the original LBM approach, the Lagrangian discretization is selected for the lefthand side [54][55][56][57].
Space and time increments are chosen to be equal to achieve the perfect shifting (streaming), where each distribution function is shifted perfectly to the neighbor node in the upstream direction.This approach is very efficient computationally and has no diffusion error.However, it works only on a structured uniform grid, which is very restrictive and might curb the applicability of the method to complex geometries.To relax this restriction, space and time increments are separated by Lee et al. [58] using the  method to discretize the Boltzmann transport equation (i.e., Equation ( 1)) along the characteristic, which then leads to the three-step method: pre-collision, streaming, and post-collision.

Pre-collision
Evidently, Equations ( 3) and ( 5) are local and denote the pre-and post-collision, respectively.However, Equation ( 4) is the non-local streaming step.The streaming step can be solved using conventional CFD methods such as FVM, FDM, or FEM.However, they are not recommended due to the complicated and pedantic mesh processing of 2D and 3D problems.Furthermore, the accuracy of those methods is considered low when compared to the RBF method for the same number of points.The RBF-based interpolation can be used to solve the streaming step (i.e., Equation ( 4)) even for an unstructured mesh.In the following, the RBF method and methodology of implementation into the LBM are presented in detail.

Radial Basis Function Method
The RBF method is used to interpolate a function by finding the coefficients of a basis function.The primary difference between the RBF and the polynomial approach is the basis (), which is a function of both distances between two points in space ((, ) = √( − ).( − )) and the shape parameter ().In the RBF method, the number of bases is the same as the number of interpolation points.For this reason, it generates a square interpolation matrix.To illustrate, the specified data (  ) at   locations can be interpolated by first defining the interpolant, where   represent the coefficients associated with each basis.The coefficients can be determined through the imposition of the interpolation constraints on the interpolant.This can be achieved by equating the interpolant (()) at points   to the specified data ((  )) at the same location.

𝐼(𝒙
Equation ( 7) is a linear system of equations that can be expressed in a matrix format The matrix [] is a symmetric and an invertible matrix if no repeated point is used.This property is found to be effective because it lowers the computational cost.Once the coefficients   are calculated, the data can be interpolated to any point  using Equation (6).Subsequently, this RBF-based interpolation can be used to solve the streaming step (i.e., Equation ( 4)) of the LBM.
Subsequently, the coefficients   are calculated from Equation ( 8) as follows: The coefficients are then used to evaluate the data at the original location .
where  () () is the interpolation matrix associated with  ℎ distribution function (  ).This calculation needs to be made once before starting the time iteration to improve the computational cost.Then, the interpolation operators ( () ()) can be used at each time step to obtain the solution of the streaming.However, it involves an unavoidable storage cost because there is a need to store huge interpolation operators ( () ()).This cost can be reduced by using sparse matrix techniques.
To explain the difference between the perfect shifting and the RBF method in handling the streaming step of the LBM, four cases for the D1Q2 model are presented in Figure 1.Case A represents the perfect streaming, where each distribution function is shifted in the original streaming step to its new position by magnitude and direction of    .The streaming step in the RBF method consists of two steps, imaginary and interpolation steps, as shown in Cases B, C, and D. The imaginary step is an unhappened step, mentioned here just to demonstrate the procedures of the RBF method.For example, we imagine that each distribution function in Case B is shifted to the imaginary position by magnitude and direction of   .Then, each one is repositioned to its correct location using the interpolation relation, Equation (11).Because the imaginary and the correct locations are the same in Case B, Cases C and D are presented.As shown in the imaginary step of Case C, each distribution function is shifted by magnitude and direction of    to the new imaginary position.To illustrate the interpolation step, the distribution function labeled 7 is resulted by using the interpolation between the distribution functions labeled 1 and 2. Case D shows that each distribution function generated in the interpolation in the direction   is produced from all neighbor distribution functions that are in the same direction weighted with appropriate weight (i.e., Equation (11)).For example, the distribution function labeled 13 is mainly produced from distribution functions labeled 1, 2, 3, 4, and 5.The weight function is monotonically decreasing.The far distribution functions (labeled 1 and 5) have less weight compared with the near distribution functions (labeled 2, 3, and 4).This makes the value of the interpolated distribution function labeled 13 closer to the values of distribution functions labeled 2, 3, and 4 than that of distribution functions labeled 1, 5. As mentioned previously, the interpolation is predicated on the evaluation of the basis function .The inverse of the matrix [], which is generated from the basis , is used to approximate a function via the interpolation operator.The shortcoming of this approach is its poor stability and hence occasionally yields inaccurate results, particularly when the shape parameter () approaches zero.A method to alleviate this problem is explained in the subsequent section.

Basis Function
The accuracy of the RBF-based method is known to be contingent on the shape parameter.In our previous work [53], the Hermite expansion with respect to the shape parameter was proposed to stabilize the RBF method and to achieve the same accuracy as the original RBF with weak or no dependency on the shape parameter.The RBF  can be expressed in terms of the Hermite polynomial by projecting the RBF on the Hermite polynomial space.
where   and   signify the  ℎ -order Hermite polynomial and its expansion coefficient, respectively.The proposed expansion of the RBF can be used to weaken the coupling between  and  by truncating the Hermite expansion at a particular order .
Notably, the Hermite polynomial is an orthogonal polynomial that satisfies the following recursion relation It is possible to obtain higher-order polynomials from Equation ( 13) by knowing the first two polynomials ( 0 () = 1 and  1 () = 2).
In the next section, these bases are used in the context of the RBF method with a view for solving the streaming step of the LBM.

Results and Discussions
In this work, 1D and 2D problems are used to evaluate the RBF method in handling the streaming step of the LBM.Furthermore, the results of the RBF method are compared with those of the regular streaming of the LBM.The analytical solutions or published results are presented to evaluate the accuracy of the RBF method and the streaming.Two mesh distributions for the RBF method are used: uniform and stretched distributions.In the uniform mesh distribution, the points are distributed uniformly along the x-and yaxis, while in the stretched mesh distribution, the points are stretched at the center and condensed at the edges of the domain to achieve a non-uniform mesh.The stretched mesh is presented to investigate the ability of the proposed method to solve problems associated with the irregular mesh.The following formula is used to map a uniform grid to a stretched grid: where erf denotes the error function.

One-Dimensional Problems
The solutions of 1D diffusion and 1D advection-diffusion equations are provided to exemplify the efficiency of the RBF combined with the modified bases for solving the streaming step.The Dirichlet boundary condition is used [2].Three mesh sizes of 100, 200, and 300 are considered for these problems.Notably, the D1Q3 lattice model can solve the diffusion and advection-diffusion equations, where it has the following lattice weights (), microscopic speeds (), speed of sound (), and equilibrium distribution (  ): , and ( 18) The advection velocity () can be substituted with zero in the case of pure diffusion.The diffusion coefficient () can be related to the relaxation time ( = 3 + /2 + ), where  denotes a parameter to control the implicitness of the method (see Equation ( 3)).

1D Diffusion
The 1D diffusion equation can be expressed in the following manner.
A step function centered at the origin is used as a benchmark for the diffusion equation within the range of (-50, 50).The analytical solution of this step function's diffusion can be written as a function of the complementary error function (erfc) [59] (, ) =  The results of the diffusion equation solved with the LBM and the conventional streaming are depicted in Figure 2. The results illustrate an excellent agreement with the analytical solution.Figure 3 presents the results of using the RBF interpolation to solve the streaming step.The figure compares the uniform and stretched mesh distributions of 100-, 200-, and 300-mesh sizes.The results of the 100-mesh size are shown to be inferior to the conventional streaming step, even for the stretched mesh distribution.This might be due to numerical diffusion, which leads to dissipation error.Increasing the mesh size to 200 or 300 reduces numerical error and raises the computational cost.The results of the uniformly distributed mesh are shown to deviate slightly from the analytical solution for the 200-mesh size, while the results of the stretched mesh reveal excellent agreement with the analytical solution.The results of the 300-mesh size improve the accuracy and achieve accuracy that is comparable to conventional streaming.Previous results are essentially qualitative, which prevents a fair comparison.Hence, it is necessary to make a quantitative comparison.Table 1 compares the error and the computational time of different mesh sizes and distributions.The error is evaluated with respect to the analytical solution.It reveals that the difference in the accuracy between the streaming and RBF interpolation decreases with an increase in the number of mesh points.To illustrate, the error per node at 300mesh is getting nearer to the error of the conventional streaming for both mesh distributions.However, the computational cost for the RBF interpolation approach is higher than that of the streaming approach.The difference is not very substantial, which demonstrates that the used approach serves as a viable alternative to the streaming step, especially when it is impossible to apply the conventional streaming as in the case of nonuniform meshes.

1D Advection-Diffusion
To investigate the proposed method furthermore, the 1D advection-diffusion equation is used to compare the RBF interpolation and the streaming.The 1D advectiondiffusion equation can be expressed as follows where  represents the advection velocity.The advection velocity represents a challenging numerical issue, where it plays the main role in the numerical stability of the method used.The same step function is used as a benchmark of the advection-diffusion equation.The analytical solution of the advection-diffusion of the step function can be written in the following manner.

𝜌(𝑥
The step function is centered at the origin within the range (-50, 50).The boundary conditions can be deduced from the analytical solution.The solution is obtained for  = 0.25  2 /,  = 0.1 /s, and different time steps up to  = 200 .
The LBM is used to solve the advection-diffusion equation that utilizes the conventional streaming, as illustrated in Figure 4.An outstanding agreement with the analytical solution is observed.Figure 5 illustrates the results of using the RBF interpolation as a feasible alternative for the streaming approach to solve the advectiondiffusion equation.The results of the interpolation are shown to be inferior to those of the conventional streaming approach.However, this gap can be reduced by increasing the number of mesh points.For example, the accuracy of the results of the 300-mesh size is comparable to that of conventional streaming with the 100-mesh size.Comparing the obtained results of the diffusion and advection-diffusion problems for the 100-mesh size, the diffusion results yielded by the RBF are more accurate than that of the advection-diffusion problem in both cases of the uniform and stretched mesh distributions.The authors believe that this huge error in advection-diffusion results for the 100-mesh size is due to the effect of the advection velocity, which causes dispersion error.However, both error sources (dispersion and dissipation) are reduced significantly by increasing the mesh size with a small computational overhead, as shown in Table 2.As a result, the proposed method is a powerful method to handle the streaming of the LBM with a comparable computational cost to that of the streaming.Table 2 compares the error and the computational time of the different mesh sizes and distributions.It demonstrates that the difference in the accuracy between the streaming and interpolation approaches can be reduced by increasing the number of mesh points.Moreover, a small difference is found in the computational cost between the interpolation approach and the streaming approach.As a result, it can be inferred that there is a low computational overhead in the utilization of the interpolation approach as an alternative to the streaming step.
In conclusion, one can say that though the present approach consumes a little bit more time compared to the perfect streaming of the LBM, it is still useful and effective in complex geometries.

2D Lid-Driven Cavity
The 2D lid-driven cavity problem is studied to show the applicability of the proposed method to higher dimensions.Notably, the D2Q9 lattice model can solve the Navier-Stokes equations, where it has the following lattice weights (  ), speeds (  ), and equilibrium distribution (  )  = [16,4,4,4,4 The bounce-back boundary condition is used for the left, right, and bottom boundaries, while the ZouHe boundary condition is used for the upper boundary (for more details, see [2]).In the present study, the Reynolds number () and lid velocity (  ) are set to be 100 and 0.1, respectively.In the current example, one uniform mesh size of 40 is used for the conventional streaming of the LBM, while 40-, 60-, and 80-mesh uniform or stretched distributions are used for the RBF approach.The obtained results are compared with that of Ghia et al. [60] for the horizontal velocity (  ) at the mid-x plane and vertical velocity (  ) at the mid-y plane.
The results of the 40-mesh size for the streaming and the RBF interpolation are shown in Figure 6.The results indicate that the RBF interpolation with a uniform mesh is better than the streaming, where the results are closer to the reference data.The error for horizontal velocity (  ) with respect to Ghia et al.'s result is 0.0410 for the streaming, 0.0353 for the RBF with a uniform mesh, and 0.1529 for the RBF with a stretched mesh, refer to Table 3.Moreover, the error for the vertical velocity (   ) is 0.0367 for the streaming, 0.03368 for the RBF with a uniform mesh, and 0.1028 for the RBF with a stretched mesh.Overall, the results show a marginal difference.However, the RBF consumes more time for convergence than the streaming.One may think that increasing the number of grid points could enhance the results; however, this does not occur in our case.For instance, increasing the number of grid points from 60 to 80 yields an error in   of 0.1214 and 0.0372, respectively, for the RBF interpolation with a uniform mesh, as shown in Figure 7 and Table 3.Moreover, the error in the stretched mesh in   is 0.1028 for the 40-mesh, 0.0448 for the 60-mesh, and 0.0532 for the 80-mesh sizes.This is due to the effect of the shape parameter on the accuracy of RBF interpolation, which does not appear for the mildly stretched mesh.In a very fine mesh, the shape parameter would require fine-tuning and updating of the matrix  () () at each time step to achieve the desired accuracy, which would increase the computational cost drastically.Hence, the accuracy of the stabilized RBF [53] should be further investigated and improved in the future for problems with very stretched mesh distribution to be less sensitive for changes in the shape parameter.This might be achieved using the local stabilized RBF approach, which uses a local support domain on each node instead of selecting all nodes.This results in a sparse matrix, which can be solved faster than a dense matrix.Hence, the computational cost would be improved.

Figure 1 .
Figure 1.Streaming step of LBM by perfect shifting and the RBF method.In all presented cases, the arrows with the BC label represent the distribution functions that are calculated based on the boundary conditions.

Figure 2 .Figure 3 .
Figure 2. Diffusion problem solved by LBM with regular streaming for 100 uniformly distributed points.

Figure 6 .
Figure 6.Lid-driven cavity solved using streaming with 40-mesh and RBF interpolation with 40-mesh uniform and stretched.

Figure 7 .
Figure 7. Lid-driven cavity solved using streaming with 40-mesh and RBF interpolation with 60-and 80-mesh uniform and stretched.

Table 1 .
Error and computational cost of the diffusion problem.

Table 2 .
Error and computational cost of the advection-diffusion problem.