Enhanced single-node boundary condition for the Lattice Boltzmann Method
Francesco Marson, Yann Thorimbert, Jonas Latt, Bastien Chopard
EEnhanced single-node boundary condition for the Lattice Boltzmann Method
Francesco Marson, ∗ Yann Thorimbert, Jonas Lätt, and Bastien Chopard
Department of Computer Science, University of Geneva, 1204 Geneva, Switzerland
We propose a new way to implement Dirichlet boundary conditions for complex shapes using datafrom a single node only, in the context of the lattice Boltzmann method. The resulting novel methodexhibits second-order convergence for the velocity field and shows similar or better accuracy than thewell established Bouzidi et al. [1] boundary condition for curved walls, despite its local nature. Themethod also proves to be suitable to simulate moving rigid objects or immersed surfaces either withor without prescribed motion. The core idea of the new approach is to generalize the descriptionof boundary conditions that combine bounce-back rule with interpolations and to enhance them bylimiting the information involved in the interpolation to a close proximity of the boundary.
I. INTRODUCTION
In recent years, the rising interest in complex flowsin numerous applications such as particulate suspen-sions [2], porous media [3], blood flow [4] and multiphaseflow [5] gave a new impulse to research on local boundaryconditions for the Lattice Boltzmann Method (LBM). Lo-cal boundary methods for curved geometries can delivera precise flow description, needing to access the flow vari-ables only on a single node located next to the surface.Thanks to these characteristics, it is possible to improvethe geometry description, yet maintain an efficient mem-ory access pattern and limit the communications betweenthreads in parallel simulations.Since the standard lattice Boltzmann method is in-herently bounded to its regular and structured lattice,boundary conditions that aim to recover realistic shapesare often of off-lattice nature. This implies that someamount of information located outside the concernedmesh node needs to be integrated into the mathemati-cal model, generally using interpolations. If this is notdone, the accuracy of the boundary condition degradesto first-order in space, due to the inability to follow thecurved shape of the wall. A first-order representation de-teriorates the overall accuracy of the simulation and mayrequire an increase of mesh resolution. For example, thecommon half way bounce-back rule [6], is an up to third-order accurate scheme [7–9], but the solution degeneratesto a first-order “stair-cased” representation when appliedto curved boundaries. In most cases, the use of interpo-lations or extrapolations causes the loss of locality of themethod. Roughly speaking, if we call boundary nodes thenodes located next to the surface, the boundary conditionwill need to access the second layer of nodes, that herewe call secondary nodes. Nevertheless, in the last threedecades some local curved boundary conditions have beenproposed [7, 10–15].To cope with the large number of different approaches,it is useful to split boundary techniques into two groups.The first is based on the computation of unknown popula-tions through a unique operation applied to the current ∗ [email protected]; [email protected] node. The second sequentially resolves the unknownsthrough independent operations in each lattice direction.These last, are often referred to as link-wise because theyoperate on “links” that connect the boundary nodes withthe wall along the discrete lattice directions. For this rea-son, they do not require any information from the otherpopulations to reconstruct the pre-collision value of thepopulation on a given link. On the contrary, the formerare named node-based .A subset of the link-wise group consists of techniquesinspired by the half-way bounce-back (HWBB) rule [6]and commonly referred to as interpolated bounceback inrecent literature [11, 16–18]. The most common HWBBextension to treat curved boundary conditions is theBouzidi, Firdaouss and Lallemand method (BFL) [1]. Inthe (linear) BFL, the populations in an adjacent layer ofnodes ( secondary nodes ) are additionally used to carryout interpolations dependent on the wall position. Thewall position impacts not only the interpolation coeffi-cients but also the choice of nodes and populations in-volved in the interpolation scheme. In the present article,we call fragmented methods with this propriety. On thecontrary, we use the term “unified” to refer to algorithmsin which the wall position only determines the coefficientsof the interpolation, but does not modify the expressioninterpolation scheme. Two years after the proposition ofthe BFL, Yu [19] proposed a scheme that can be consid-ered as a unified version of the linear BFL, given thatit uses the same populations as the BFL to perform theinterpolation. The BFL and Yu methods share the draw-back of not being local, needing a second layer of nodesto operate.In the last two decades, few attempts to create local in-terpolated bounce-back schemes have been proposed [11–14]. Among those, the Zhao and Yong (ZY) [12] inge-niously chose the LBM populations used to interpolatethe unknown mixing pre-collision and post-collision val-ues to build a single-node boundary condition, withoutthe need to introduce further elements to the model. Thisboundary condition, further developed in references [20–22], is second-order accurate in space under diffusive scal-ing hypothesis and first-order accurate in time [20]. TheZY method has been tested by Peng et al. in [16, 23].Just like the original paper, they report second-order con- a r X i v : . [ phy s i c s . c o m p - ph ] S e p vergence and accuracy similar to the BFL. The remaininglocal interpolated bounce-back methods [11, 13, 14] followa different approach introduced by Chun and Ladd [11]:first they reconstruct the boundary populations using thewall velocity and a constant approximation for the den-sity, and then they use an approximated version of thenon-equilibrium bounce-back of Zou and He (see [24]).However, the method of Chun and Ladd (CL) is notstrictly speaking local, as it can require information fromother nodes in some situations. In the wake of Chun etLadd, in recent years two new local boundary conditionshave been proposed: the Tao et al. [13] and the Liu etal. [14].Among the first local link-wise techniques, the Filip-pova and Hanel (FH) [25, 26] can be singled out. The FHbelongs to the family of ghost methods [9], also known as extrapolation methods [27] or fictitious equilibrium meth-ods [28], that uses additional fictitious nodes on the solidside of the boundary together with extrapolations to re-construct the unknown populations. In FH the ghostnode used for the interpolation is built guessing a veloc-ity beyond the wall with an extrapolation of the bound-ary velocity. Unfortunately, the FH method has knownstability issues, solved by the Mei, Li and Shyy, but sac-rificing the locality of the method [9, 29].Junk and Yang (JY) [10] proposed a single node boundary condition based on a correction of the half-waybounce-back scheme. To perform the correction, it is nec-essary to solve a linear system on each node to ensure thecompliance of the numerical result with the expected so-lution at Navier-Stokes level. For this reason, it cannotbe considered link-wise . Besides, the solution of the lin-ear system for each node adds a layer of complexity tothe implementation. This approach leads to an almostlocal mass conservative boundary condition that showssecond-order convergence for the velocity. However, theoriginal paper [10] reports that the JY method is slightlyless accurate than BFL for the velocity and pressure fieldsin the case of the flow inside a cylinder. The method hasbeen extensively tested by Yang [30] that concluded thatthe method is almost as stable as the BFL method, andit has comparable or better accuracy of the BFL and FHmethods. Nevertheless, the JY method has been testedby Nash et al. who reported in [31] poor stability prop-erties.It is worth mentioning two other local non- link-wise methodologies. The partially saturated bounce-back(PSBB) [32, 33] is a local method based on full-waybounce back. The full-way bounce-back rule is based on amodified LBM collision step, while the half-way bounce-back scheme modifies the streaming step. In the PSBB,the underlying idea is to use the knowledge of the frac-tion of fluid in the boundaries cells to operate a mixedfluid-solid collision [9, 33]. The PSBB is exactly mass-conservative and allows for a smooth transition betweensolid and fluid nodes in the case of moving objects. More-over, it does not require the exact knowledge of the shapeof the surface: this is particularly suitable for porous me- dia applications. However, this can turn into a disadvan-tage when it is necessary to guarantee an exact no-slipcondition at the surface, because the fluid fraction is nota sufficient piece of information for the method to “know”the boundary position and orientation. It also requiresan additional computation step if the fluid fraction ofboundary nodes needs to be recovered from the geomet-rical shape of the wall. Furthermore, Chen et al. compar-ing different boundary conditions in [34] reported a lowaccuracy in the computation of the cylinder drag whenusing the PSBB. From the algorithmic point of view, thePSBB cannot be used to represent thin shells becausethis method constrains the user to allocate solid nodes inthe simulation.We finally mention the Local Second-Order Bound-ary (LSOB) method of Ginzburg [7]. The technique isbased on a precise computation of the boundary nodesaccording to the Chapman-Enskog expansion to relatemacroscopic fields such as density and velocity with themesoscopic populations of the LBM. The LSOB is a high-fidelity third-order accurate local method, but it is lim-ited to laminar flows [15] and its implementation is latticeand problem-dependent.Despite the existence of the single-node boundary con-ditions that we have reviewed, further research is neededto make local boundary conditions appealing in terms ofaccuracy and simplicity of the implementation. Further-more, the relation between the existing local boundarymethod should be clarified. To this end, in the presentarticle, we develop a framework to generalize interpolatedbounce-back schemes, including the CL, Tao, Liu, ZY lo-cal methods, and the well-established BFL and Yu meth-ods. Within this framework, we also develop a family ofnovel boundary conditions to improve the compactnessof the interpolation range and the accuracy of the non-equilibrium approximation adopted in the CL, Tao, andLiu methods.This article is structured in the following way. Af-ter briefly presenting the LBM in section II A and theinterpolated bounce-back methods in section II B, a gen-eral description of the local ELIBB is presented in sec-tion II C. In section II D some specific variants of thegeneral ELIBB scheme are proposed. Finally, the imple-mentation in the open-source software PALABOS [35]of the ELIBB is tested for three configurations whoseanalytical solution is known. Namely, the impulsively-started unsteady Couette flow in section III A, the steadycylindrical Couette flow in section III B and the Jeffery’sorbit in section III C. II. NUMERICAL METHODOLOGYA. The lattice Boltzmann Method
The Boltzmann equation (BE) ∂ t f + ( ξ · ∇ ) f = Q BE ( f, f ) (1)describes the space and time evolution of the probabilitydistribution function f ( x , ξ , t ) of finding a particle withvelocity ξ at position x and time t . The latter is sub-ject to advection in the velocity space ξ , as well as col-lision as illustrated by the Boltzmann’s collision integral Q BE ( f, f ) [36–39]. The complexity of the Boltzmann’scollision integral is the major obstacle to the solutionand analysis of the equation. This is why Q BE ( f, f ) iscommonly approximated, with relaxation towards equi-librium models. One of the oldest and most successfulrelaxation models is the BGK operator that was formu-lated independently by Bhatnagar, Gross, Krook [40] andby Welander [41].In the LB context, the BE (1) is first decomposed ina finite set of equations resulting from the velocity spacediscretization [42, 43]. Those equations are known as theDiscrete Velocity Boltzmann Equation (DVBE) and readas ∂ t f i + ( ξ · ∇ ) f i = Q DVBE ( f i , f i ) (2)Similarly to the BE, the DVBE expresses the timeevolution of discrete probability distribution functions f i ( x , ξ i , t ) . But contrarily to the BE, the latter “pop-ulations” now propagate at constant velocity ξ i with i ∈ { , , . . . , Q } . Applying the method of the character-istics and the trapezoidal integration rule –to the LHSand RHS terms of the DVBE (2) respectively– one endsup with the Lattice Boltzmann Equation (LBE) that isthe cornerstone of LBMs [9, 44] f i ( x + c i , t + 1) = f i ( x , t ) − Ω f neq i (cid:124) (cid:123)(cid:122) (cid:125) Q LBM (3)where lattice units are implied. Q LBM consists in aheuristic approximation of Q BE , and Ω is the relaxationparameter. f eq i and f neq i = f i − f eq i are the equilib-rium and non-equilibrium populations respectively. Inpractice, the LBE (3) is solved through two successivesteps. The local collision step (RHS of equation (3)) andthe non-local Lagrangian streaming step (LHS of equa-tion (3)). Through the normalization of discrete veloci-ties ( c i = ξ i c s with c s the lattice constant), the stream-ing satisfies the on-grid condition that leads to an exact advection of populations from one grid node to anotherone [9]. To compute the collision term, one needs to de-rive equilibrium and non-equilibrium populations. f eq i isthe discrete counterpart of the Maxwell-Boltzmann dis-tribution [39] f eq = ρ ( x , t )(2 πθ ( x , t )) D/ exp (cid:18) − ( u ( x , t ) − ξ ) θ ( x , t ) (cid:19) , (4)where D is the dimensionality of the problem. ρ , u and θ are the macroscopic density, velocity and reduced tem-perature respectively. One way to derive f eq i from theMaxwellian (4) is to rely on the Gauss-Hermite quadra-ture [42, 43]. In the context of isothermal LBMs, thereduced temperature θ = T is constant and directly re- lated to the lattice constant T = c s , where c s = 1 / forstandard velocity discretizations (D2Q9, D3Q27, etc).Regarding collision models, one of the most popular isthe BGK approximation that leads to [45] Ω BGK = 1¯ τ + 1 / , f neq i = (cid:0) f i − f eq i, (cid:1) (5)where ¯ τ is the relaxation time, even though, in the LBcommunity, it is also common to refer to τ = ¯ τ + 1 / asrelaxation time. f eq i, = w i ρ (cid:20) c iα u α c + u α u α ( c iα c iα − c δ α α )2 c (cid:21) (6)is the discrete equilibrium up to the second-order, where α i ∈ { x , . . . , x d } ∀ i and Einstein’s summation rule isassumed for the sake of compactness. The BGK collisionmodel is easy to understand and to implement, but it hasseveral drawbacks. Concerning the bounce-back method,the coupling with the BGK collision model gives rise toa second-order error (related to the exact location of thewall) that is commonly referred to as viscosity-dependenterror . Among the other limitations, it can lead to nu-merical instabilities in the case of under resolved meshes(typically at high Reynolds numbers) and non-vanishingMach numbers [9, 27, 46].To overcome these issues, several extended collisionmodels were proposed [47] with varying degrees of suc-cess [46]. As a first study, hereafter, we will restrictourselves to Multi-Relaxation-Time (MRT) [48–50], Two-Relaxation-Time (TRT) [51–54], and regularization [55–59] approaches, which all have a different impact on theboundary condition performances [60]. a. MRT models: The core idea of the MRT is tocarry out the relaxation collision process in the momentspace. In other terms, from Q populations, Q momentsare computed. The latter is (independently) relaxed to-wards their equilibrium values using different relaxationparameters. Referring to equation (3), we can mathe-matically describe the MRT collision model as follows Ω MRT = M − SM , f neq i = (cid:0) f i − f eq i, (cid:1) (7)where, in the original formulations, S is a diagonal re-laxation matrix, while M and M − are two orthogo-nal matrices allowing to move from the Q-dimensionalvelocity space to the Q-dimensional orthogonal momentspace and vice versa . In this context, the LBE (3) can berewritten for the MRT case using the bra-ket notation | f i ( x + c i , t + 1) (cid:105) = | f eq i (cid:105) + (cid:0) I − M − SM (cid:1) | f neq i (cid:105) (8)where the symbol | . . . (cid:105) stands for the vectorial notationpopulations. This method depends on many free param-eters, though, which need to be handpicked through anadequate procedure. b. TRT models: To reduce the number of free pa-rameters, Ginzburg et al. [51–54] proposed a two-relaxation formulation of the MRT collision model (8).These TRT-LBMs independently relax even (symmetric)and odd (anti-symmetric) moments, each of which areassigned an individual collision frequency: ω + and ω − respectively. Interestingly, the TRT behavior can be de-scribed by the magic parameter Λ that relates the twocollision frequencies in the following way Λ = (cid:18) ω + − (cid:19) (cid:18) ω − − (cid:19) . (9)TRT-based collision models are particularly interestingfor the bounce-back boundary condition. It has beenproven that for Λ = / the bounce-back conditionis viscosity-independent meaning, i.e., the wall locationdoes not depend on the viscosity anymore. The latterfeature is particularly critical for low-viscosity flows, andit further leads to third-order accuracy for the half-waybounce back [9]. c. Standard and recursive regularized models: Theinstabilities of the BGK model can be addressed in away that is an alternative to the MRT and TRT ap-proach (even if can be reformulated in the same formal-ism for particular moment spaces and relaxation frequen-cies [46, 47]). The basic idea of the regularized (RBGK)and the recursive regularized (RRBGK) models is to fil-ter out non-hydrodynamic modes in the BGK solutionarising from the numerical discretization of the veloc-ity space. The latter spurious contributions are hiddenin non-equilibrium populations and can be deleted by“manually” imposing a particular form –compliant withthe macroscopic behavior of interest– of f neq i . More pre-cisely, starting from a Knudsen number expansion of pop-ulations (as done, e.g., for any Chapman-Enskog expan-sion), one obtains f i = ∞ (cid:88) l =0 f ( l ) i = f (0) i + f (1) i + f (2) i + . . . (10)with f ( l ) i ∼ O (Kn l f (0) ) the l -order contribution (withrespect to the Knudsen number Kn) to populations f i ,and f (0) i = f eq i [61]. Each regularization step aims atmanually discarding terms that are not compliant withthe Navier-Stokes-level of physics, i.e., f neq i = f i − f (0) i ≈ f (1) i (11)where f (1) i can have several forms depending on the as-sumption of the considered regularized approach [55–59].An in-depth discussion regarding that point can be foundin appendix B and appendix C. B. Interpolated bounce-back methods
The purpose of boundary algorithms is to reconstructmissing populations on nodes next to the wall ( bound-ary nodes ), after the streaming step. On the boundarynodes
F, the unknown population are those associatedwith discrete velocities c i that “leave” the wall (see fig-ure 1). We call such populations “ incoming to fluid ”and denote them with index i , while we call “ outgoingfrom fluid ” and denote with the index ¯ i the others [27].The main idea of the Interpolated Bounce-Back methods(IBB) is to perform a one dimensional polynomial inter-polation of the known population close to the boundariesalong the discrete directions of the lattice (links direc-tions) to recover the unknown incoming populations atthe boundary nodes . We can formalize this general ideain mathematical terms writing f i ( x F , t + 1) = (cid:88) j a j f ∗ i ( x j , t ) + (cid:88) k a k f ∗ ¯ i ( x k , t ) + K (12a)where f ∗ = f − Ω f neq is the post-collision population, K is a hypothetical correction factor, the symbol a de-notes the interpolation coefficients, x is an interpolationpoint and t is the current iteration. In practice, for thelinear case and referring to figure 1, the previous formulagenerally reduces to f i ( x F , t + 1) = a f ∗ ¯ i ( x F F , t ) + a f ∗ ¯ i ( x F , t )+ a f ∗ i ( x F , t ) + K (12b)equations (12) represent a generic formulation of theinterpolated-bounceback approach. To derive a specificFigure 1: 2D representation of the boundary nodes,normalized distance q , discrete lattice velocities (cid:126)c i , (cid:126)c ¯ i ,links (dashed segments) and locations of boundaries atthe intersections with links ( • ). ∆ x is the space step.method, it is necessary to specify the expressions of theinterpolation coefficients and points. To this end, twoviable solutions exist. The first one consists in writingclosure relations by exploiting the macroscopic no-slip condition [8, 12, 62, 63]. In practice, the no-slip con-dition is expanded using a formal mathematical expan-sion and subsequently equations (12) are injected ontoit. The second solution relies on a mesoscopic, geometri-cal approach and was proposed by Bouzidi et al. amongother authors [1, 11, 13, 19]. The idea in this case is touse the bounce-back rule, understood as a modificationto the streaming step, to compute the interpolation co-efficients. Roughly speaking, the bounce-back operatormodifies the streaming operator from a simple translationin space to a translation-reflection-translation. The pop-ulations subjected to the bounce-back, during the trans-lation, are reflected when they encounter the wall. Owingto this bounce-back rule, the interpolation coefficients arethose that allow to geometrically compute the unknowneither at time step t + 1 or at its virtual off-lattice post-collision state at time t .To illustrate this concept we consider the BFL algo-rithm [1]. For the linear BFL algorithm, equation (12b)becomes f t +1 i ( x F ) = a (cid:122)(cid:123) q f ∗ ¯ i ( x F ) + a (cid:122) (cid:123) (1 − q ) f ∗ ¯ i ( x F F ) q < . (13a) f t +1 i ( x F ) = 12 q (cid:124) (cid:125) a f ∗ ¯ i ( x F ) + 2 q − q (cid:124) (cid:125) a f ∗ i ( x F ) q ≥ . (13b)where q denotes the distance of the boundary node Ffrom the wall, normalized by the norm of the discretevelocity (cid:107) c i (cid:107) = (cid:107) x F F − x F (cid:107) . As showed by figure 2, inthe IBB methods, populations can be thought as lumpedmass elements moving according to their discrete latticevelocities c i . During the streaming step, each populationundergoes either to the free streaming (along a straightline) or to the bounce-back streaming (in case of a wallencounter). In the latter case, the populations reverttheir streaming direction before completing their trajec-tory of the length of (cid:107) c i (cid:107) . When q < / (figure 2a), theinterpolation is carried out when all populations are atthe time-step t in their post-collision state. In this case,the target location of interpolation is the former positionof f t +1 i ( x F ) before the bounce-back streaming step. Onthe contrary, when q ≥ / (figure 2b), the interpolationfactors must be computed after the streaming procedure.Consequently, either the streaming or the bounce-backrule is applied to each population at time-step t (post-collision state). After that, the interpolation is carriedout at time-step t + 1 . Figure 2: 1D link-wise representation of bouncing-backprocedure in the BFL method: (left) q < / and(right) q ≥ / . The coefficients a , a , a are defined inequations (13). The arrows represent the populationsappearing in equations (13). The dashed curved arrowsrepresent the bounce-back rule. Left: the streaming andthe bounce-back rule are applied after the interpolationof the off-lattice population. Right: the streaming andthe bounce-back rule are applied before theinterpolation, that is carried out at timestep t + 1 . C. Single-node interpolated bounce-back methodsand the new ELIBB scheme
In order to design local link-wise BCs, one must dis-card the non-local contribution f ∗ ( x F F ) that appears inequations (12b) and (13). From the algorithmic pointof view, one could simply apply the boundary methodafter streaming and save the other population at nodeF before the streaming. Unfortunately, this algorithmiclocality leads to some issues. In fact, it does not allowto describe corners or narrow gaps (where a second wallis located between the nodes F and FF without intro-ducing a special treatment of these cases [11]). To getrid of the unwanted f ∗ ( x F F ) population, yet maintain-ing the link-wise nature of the method, two approacheshave been proposed. The first one comes from Zhao etal. [12, 20] and consists on the following first-order intime approximation f ∗ ¯ i ( x F F ) = f t +1¯ i ( x F ) ≈ f t ¯ i ( x F ) . (14)The second one, was introduced by Chun and Ladd [11],further developed by Tao, Liu et al. [13] and also im-proved by the present article. It consist building virtual(approximated) population located at the wall position x W . In the proposition of Tao et al. [13], only the popu-lation ˜ f t +1 i ( x W ) was introduced in equation (12b), wherethe ∼ cap indicates an approximated value. Hereafter,we propose a novel Enhanced Local Interpolated Bounce-Back method (ELIBB) which additionally accounts forthe population ˜ f ∗ i ( x W ) hence, extending (12b) in the fol-lowing way E L I BB (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) f i ( x F , t + 1) = (cid:24)(cid:24)(cid:24)(cid:24)(cid:24)(cid:24) a f ∗ ¯ i ( x F F ) + a f ∗ ¯ i ( x F ) (cid:124) (cid:123)(cid:122) (cid:125) Ladd’s BB + a f ∗ i ( x F ) + a ˜ f ¯ i ( x W , t + 1) (cid:124) (cid:123)(cid:122) (cid:125) by Tao et al. + a ˜ f ∗ ¯ i ( x W ) (cid:124) (cid:123)(cid:122) (cid:125) new ELIBB . (15)The previous equation represents the scheme underly-ing the new ELIBB method. However, it is still neces-sary to specify the expression of the populations and theinterpolation coefficients. It will be done in the follow-ing section II C 2. Changing the interpolation coefficientsin equation (15), different variants will be developed (sec-tion II D). From this perspective, the model proposedby Tao et al. [13] directly flows from the ELIBB (equa-tion (15)).
1. Approximation of the wall populations
In the approach introduced by Chun in [11], the fun-damental idea to realize a single-node IBB is to exploitthe knowledge of the boundary velocity to rebuild some virtual populations at the boundary location. This ap-proach consists of using the boundary velocity and anestimation of the density at the boundary to rebuild theequilibrium population. In [11], it has been proved thatfor a “slow” flow the approximation ˜ f eq ,t +1 i ( x W ) ≈ f eq ( ρ ( x F , t ) , u ( x W , t + 1)) (16)is second order accurate. Regarding the approximation ofthe non-equilibrium component at the wall ˜ f eq ,t +1 i ( x W ) we can use an approximated non-equilibrium bounce-back ([11, 13, 14]). This is a first-order approximation of the non-equilibrium bounce-back method of Zou and He [24].We discuss it in appendix C [see in particular equa-tion (C10)]. This leads to the following second orderaccurate approximation ˜ f neq ,t +1 i ( x W ) ≈ f neq¯ i ( x F ) (17)As demonstrated in [11] the reason why the approxi-mated non-equilibrium bounce-back leads to a second-order accurate boundary condition resides in the factthat the non-equilibrium component is a second-ordercorrection over the equilibrium. This fact allows for asecond-order approximation with an only first-order ap-proximation of the non-equilibrium part.In the present article, we additionally propose a moregeneral approach to estimate the non-equilibrium com-ponent. The idea is to use Malaspinas’ [56] recursive for-mulas to recompute the wall non-equilibrium componentusing the Hermite basis expansion truncated up to thefourth-order (refer to appendix C for the details of this procedure). This approach allows for higher flexibility inthe modeling of the wall non-equilibrium population.
2. Generalized computation of interpolated coefficients inIBB methods
To obtain a more uniform picture, we propose to ex-press the populations in all these IBB methods in their pre-collision state at time t + 1 in a similar fashion ofthe case q ≥ / of the BFL method (figure 2). To beable to extend this description to all methods and for anyvalue of q , we introduce signed normalized distance fromthe wall s at time t + 1 . At the time t + 1 some popu-lation has been streamed following the free stream rule,whilst, others near the wall have been streamed using thebounce-back rule described in section II B and figure 2b.In this condition the generalized coordinate s reads s ( f I ( x , t + 1)) def = ( x − x W ) · c I (cid:107) c i (cid:107) ∀ I ∈ { i, ¯ i } , (18)where x is the coordinate of the population f I after thestreaming/bounce-back step. The coordinate s turns outto be a simple yet effective tool to describe and comparelink-wise boundary conditions. Using equation (18), wecan define a set of simple rules to move from the x coor-dinate metric to the s coordinate metric: f t +1 i ( x F ) ≡ f ( s = q ) f t, ∗ i ( x F ) = f t +1 i ( x F F ) ≡ f ( s = q + 1) f t, ∗ ¯ i ( x F ) ≡ f ( s = − q + 1) f t, ∗ ¯ i ( x F F ) = f t +1¯ i ( x F ) ≡ f ( s = − q ) f t, ∗ ¯ i ( x W ) ≡ f ( s = 1 + u w · c i / (cid:107) c i (cid:107) ) f t +1¯ i ( x W ) ≡ f ( s = 0) , (19a)(19b)(19c)(19d)(19e)(19f)where the index of x indicates either a node (F or FF)or the virtual node W located at the intersection of thelink with the boundary surface (as a consequence f ( x W ) are some virtual population which will be useful later)and u w is the wall velocity (in lattice units).With the help of the coordinate s , it is now possibleto generalize the formulas used in the BFL and the otherIBB methods with the Sylvester-Lagrange [64, 65] poly-nomial interpolation formula: f ( s ref , t + 1) = n (cid:88) j =0 a j ( s j ) f ( s j , t +1) (20)where n is the interpolation order, j is the index of theinterpolation point, and a j are the interpolation coeffi-cients given by a i = (cid:89) ≤ j ≤ n j (cid:54) = i s ref − s j s i − s j . (21)In the linear case the interpolation coefficients a j ∈{ a α ( s α ) , a β ( s β ) } can be easily recovered from the valuesof s j ∈ { s α , s β } in the following way a α ( s α ( q )) = 1 − s ref − s α s β − s α (22a) a β ( s β ( q )) = s ref − s α s β − s α (22b) s α < s ref < s β (22c)where in our case s ref = q . equations (22) will be usedto develop variants of the general ELIBB formula (equa-tion (15)). D. Enhanced Single-Node boundary conditionvariants
The novel generalized coordinate introduced in the pre-vious section (equation (18)) is a useful tool to developvariants of the general scheme proposed in equation (15).The first method that we propose is a unified method,and it can be written as:ELIBB-U | f t +1 i ( x F ) = a ˜ f t +1 i ( x W )+ a ˜ f ∗ i ( x W ) (23a)where ˜ f is the approximation of population f consist-ing of a separate evaluation of the equilibrium ˜ f eq andnon-equilibrium ˜ f neq parts that was discussed in sec-tion II C 1. As anticipated in the introduction, the ad-jective unified refers to the fact that the interpolationscheme does not depends on the value of q . The in-terpolation factors a , in equation (23a) can recoveredconverting the x coordinates in s coordinates with equa-tions (19) and then using equations (22a) and (22b) a = 1 − q (23b) a = q . (23c)For the wall populations we adopted the following ap-proximations ˜ f ∗ i ( x W ) def = ˜ f eq i ( x W ) + ˜ f neq ∗ i ( x W )˜ f t +1 i ( x W ) def = ˜ f eq ,t +1 i ( x W ) + ˜ f neq ,t +1 i ( x W )˜ f eq ,t +1 i ( x W ) ≈ f eq ( ρ F ( t ) , u ( x W , t + 1))˜ f eq i ( x W ) ≈ f eq ( ρ F ( t ) , u ( x W , t ))˜ f neq i ( x W ) ≈ use (17) or (C8) , (24a)(24b)(24c)(24d)(24e)where in (24c) and (24d) we used (16). In this re-gard, one may note that in general x W ( t ) (cid:54) = x W ( t + 1) and u ( x W ( t )) (cid:54) = u ( x W ( t + 1)) , but their actual expres-sions depend on the time advance scheme of the wall( e.g. explicit euler, implicit euler, etc. ). Nonetheless, if ∆ u w (cid:28) u w one can set u ( x W ( t )) ≈ u ( x W ( t + 1)) . Thisfirst variant of the ELIBB method is different from theother because it does not use any nodal population in the interpolation. In particular, if the non-equilibriumcomponent at the wall is computed from a regularizationprocedure ( e.g. equation (C8)) instead of using equa-tion (17), a good name for this method could be localregularized link-wise method . This name emphasizes thefact that this specific variant only relies on the recon-struction of the wall populations just like the regularizedboundary condition presented in the articles [66, 67].The physical locality of (23a) can be improved adopt-ing a fragmented interpolation scheme, taking advantageof the knowledge of the population f ∗ ¯ i ( x F ) , E L I BB - F (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) f t +1 i ( x F ) = a ˜ f t +1 i ( x W )+ a f ∗ ¯ i ( x F ) q < . f t +1 i ( x F ) = a f ∗ ¯ i ( x F )+ a ˜ f ∗ i ( x W ) q ≥ . (25a)(25b)where: a = (1 − q ) / (1 − q ) q < . (25c) a = q/ (1 − q ) q < . (25d) a = (2 − q ) /q q ≥ . (25e) a = (2 q − /q q ≥ . (25f)Finally, the combination of all the information ofthe two previous methods, it is possible to perform aquadratic interpolation. That leads to a quadratic inter-polation unified method E L I BB - U Q (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) f t +1 i ( x F ) = a ˜ f t +1 i ( x W )+ a f ∗ ¯ i ( x F )+ a ˜ f ∗ i ( x W ) (26)where the interpolation coefficients a i are computed us-ing equations (21) (after using (19) to compute the s coordinates).Adding the knowledge of ˜ f ∗ i ( x W ) we can also formalizea fragmented quadratic method E L I BB - F Q (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) f t +1 i ( x F ) = a ˜ f t +1 i ( x W ) q < . a f ∗ ¯ i ( x F )+ a ˜ f ∗ i ( x W ) f t +1 i ( x F ) = a f ∗ ¯ i ( x F ) q > . a ˜ f ∗ i ( x W )+ a f ∗ i ( x F ) (27a)(27b)where the interpolation coefficients are given again byequation (21). The above quadratic interpolation for-mulas defined by equations (27) and (26) improve theaccuracy of the ELIBB-F. Nevertheless, it is necessaryto specify a cutoff minimal value q mincutoff below which theELIBB-F is used in place of the ELIBB-FQ. Thereby, weavoid instabilities due to the overlapping of interpolationpoints when q approaches zero. For the results presentedin this article we use q mincutoff = 0 . .We remark that the equations seen in this section arevalid for resting boundary with respect to the lattice. Inthe case of a moving boundary, it is necessary to adda momentum correction. We discuss this in the nextsection. E. Moving Boundaries
In the case of a moving boundary, it is necessary tointroduce a moving boundary correction at wall nodes,as proposed for the HWBB in in [6, 68]. This correctionterm is used to generalize to moving boundaries the IBBseen in the previous paragraphs modifying each nodal (non-wall) population in the rhs of the equations withthe following simple rule, written in terms of generalizedcoordinate s if s < then f ( s ) → f ( s ) + 2 w i ρ c i · u w c s , (28)where u w is the velocity of the wall at the link intersec-tion x W . F. Force computation
There exist two main ways to compute the force andtorque acting on the surfaces: the stress tensor integra-tion [69] and the Momentum Exchange Algorithm [6]. Asrecommended by Mei et al. [70], the momentum exchangealgorithm is used for the computation of the force in thefollowing section.There are different numerical approaches to computethe momentum exchange from the fluid to the surface [6,70–75] that have been summarized by Tao et al. in [76].For the experiments in this paper we have decided toadopt the method described in [75]. In [75] Wen et al. propose the following receipt to compute the force on theboundary of on lattice node F = (cid:40)(cid:88) i − (( c i − u w ) f i ( x f ) − ( c ¯ i − u w ) f ∗ i ( x f )) (cid:41) out (29)where F is the force acting on the boundary surface dueto one fluid lattice and the index “out” means that whenthe boundary is described as a closed surface, the forcecomputation has to be carried out only for the “external”fluid. III. RESULT AND DISCUSSION
The variants of the novel boundary condition methodhave been implemented for the D2Q9 and D3Q27 lat-tices, in the open-source library PALABOS [35]. In the results, we mainly used the non-equilibrium bounce-back approach to estimate the non-equilibrium compo-nent (equation (17)) because it is simpler to implementand gives similar results of the regularized approach (ap-pendix C). Nevertheless, in the second test-case, we in-cluded a brief comparison between the results obtainedwith the two different methods. All the physical quan-tities appearing in the next sections are either in latticeunits , or in non-dimensional units specific to each prob-lem.
A. Impulsively-started unsteady Couette flow
The impulsively started Couette flow is the fluid config-uration obtained abruptly moving one of the two parallelwalls containing a quiet fluid, from the rest position tothe constant velocity U . In this specific case, we considerthe upper wall moving along the x direction and locatedat y = h and the bottom one resting at y = 0 .In the context of low Reynolds number flows, con-vective phenomena are negligible for this configuration,whose evolution is then governed by [77] ∂u∂t = v ∂ u∂y (30)with the following boundary and initial conditions u (0 , t ) = 0 for all tu ( h, t ) = U for t > u ( y,
0) = 0 for (cid:54) y < h . (31)The problem defined by equations (30) and (31) has asolution in the form of slow converging series [77] u ∗ th = u ∗ th U = yh + 2 π (cid:88) n =1 ( − n n e − n π νt/h sin nπyh . (32)In [77] Erdoˇgan mentioned that at time t ∗ = νt / h equa-tion (32) truncated at the 45th term is sufficient to obtaina numerical solution compatible with a double floatingpoint precision numerics simulation (given that the eval-uation is computationally very light, we used 100 termsin the present test case).We used the ELIBB variants, the HWBB, the BFL,and the Tao’s methods to simulate the impulsively-started Couette flow. The numerical domain is squared,bounded in the y direction by the walls, and by the peri-odicity condition in the x direction. The LBM simulationis carried out with the D2Q9 lattice and the BGK col-lision model. The simulations are symmetric in the twodimensions and the top and bottom layer of nodes arelocated at y = q and y = h − q (lattice units).The relative Root Squared Error rRSE (also known asFigure 3: Convergence of rRSE for impulsively started Couette flow at dimensionless time t ∗ = νt / h . Results fordifferent Reynolds number Re , relaxation time τ , and normalized distance of the walls from the first layer of nodes q . The LBM topology is D2Q9, the collision model is BGK.0Figure 4: Representation of computational domain. Ingreen the inter cylinders Couette flow. In the externalyellow region, the fluid is at rest because the externalcylinder of radius r is at rest, in the inner yellow regionthe flow is rotation, moved by the inner cylinder ofradius r that is moving with angular velocity ω . Boxeddots represent the nodes where the error along thex-axis is computed. L -error function), defined as rRSE ( u ∗ ) = (cid:118)(cid:117)(cid:117)(cid:116) (cid:80) Ni ( u ∗ ( x i ) − u ∗ th ( x i )) (cid:80) Ni u ∗ ( x i ) , (33)is used to evaluate the convergence of computational er-ror in the cylindrical Couette flow region, where x i isthe coordinate of a lattice node, u ∗ is the computed non-dimensional macroscopic velocity norm and u ∗ th is thetheoretical velocity norm given by equation (32).The results for the rRSE for increasing height h , ex-pressed in lattice units are shown in figure 3. For smallrelaxation numbers ( Re = 330 , τ = 5 . ) all methodsshow a second-order convergence and almost identical ac-curacy but the HWBB because q (cid:54) = 1 / . For τ in therange [0 . , . the following observations can be made:1. quadratic variants of ELIBB (ELIBB-UQ, ELIBB-FQ) are generally more accurate than other meth-ods;2. ELIBB-UF and ELIBB-FQ performes identicallyfor q < / because they have the same interpo-lation scheme in this range;3. when q ≥ / ELIBB-FQ appear to be more accu-rate than ELIBB-UQ for τ (cid:46) , but vice-versa for τ (cid:38) ;4. ELIBB variants are generally more accurate thanthe Tao’s method in this range. x L E ( u y ( x c )) (a) x × − L E ( u y ( x c )) (b) Figure 5: Plot of the linear error LE = (cid:12)(cid:12)(cid:12)(cid:12) u y ( x c ) − u theory y ( x c ) u ref y (cid:12)(cid:12)(cid:12)(cid:12) along the x-axis for thecylindrical Couette flow characterized by β = / . (a) r − r = 10 rmlu , (b) r − r = 40 rmlu The verticallines serve as markers of the x coordinate ofintersections with the cylinders: x < . and x > . rest flow region, . < x < . and . < x < . cylindrical couette flow region, . < x < . innerregion (linear flow). B. Steady state cylindrical Couette flow
The cylindrical Couette flow is a common benchmarkto test the accuracy of curved boundary conditions. Weimplemented this test case using a D3Q19 lattice. In thistest case, two coaxial cylinders are placed in the centerof the simulation domain. The cylinders axis is parallelto the z direction, along which the periodicity conditionhas been imposed. The inner cylinder of radius r rotateswith angular velocity ω = u θ /r and u θ = ν Re/ ( r − r ) (34)1Figure 6: Cylindrical Couette flow, convergence of the relative root squared error rRSE ( i.e. L error funciton)for the fragmented variants of the present schemes(ELIBB-FL, ELIBB-FQ) and for the Bouzidi et al. (BFL) method. Parameters: TRT collision model with Λ = / , β = 0 . .tangential velocity, while the outer of radius r is at rest.The inter-cylinder distance r − r can be expressed usingthe cylinder ratio β = r /r parameter as r (1 /β − .The velocity flow inside the inner cylinder is linear ifthe flow is laminar, with a maximum tangential velocity u θ = 1 close to the inner cylinder and u θ = 0 in thecenter. Between the two cylinders, the solution for thefluid velocity at steady state is given by the cylindrical-Couette flow tangential velocity u θ = (cid:0) r − β r (cid:1) ω (1 − β ) r (35)where r is the radial distance from the axis of the cylin-ders.The relative Root Squared Error rRSE, defined byequation (33) is used to evaluate the convergence of com-putational error in the cylindrical Couette flow region(see figure 4), wherein this case u ∗ th = u θ is the theoreti-cal velocity norm. The rRSE is also used to evaluate theerror of the tangential velocity along the x-axis, in this Figure 7: Cylindrical Couette flow, convergence of the relative root squared error rRSE ( i.e. L error function)for the unified variant of the present schemes(ELIBB-UL) and for the Yu et al. method. Parameters:TRT collision model with Λ = / , β = 0 . .case, the expression reads rRSE ( u θ ( y c )) = (cid:118)(cid:117)(cid:117)(cid:116) (cid:80) Ni ( u θ ( x i , y c ) − u th ( x i , y c )) (cid:80) Ni u ( x i , y c ) , (36)where y c is the y − coordinate of the cylinders axis andthe index i refers to the index of the node along the x -direction. To visualize the numerical error along the x -axis, we also use the linear error (LE), that we defineas: LE ( x i ) = (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) u y ( x i , y c ) − u th y ( x i , y c ) u th y (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) , (37)where x i is the x coordinate of a node located at y c ,which is the y coordinate of the center of the cylinders.
1. Space distribution of the linear error figure 5 show the trend of LE for two resolutions r − r ( lu and lu) and for seven different boundary condi-tions: the BFL method [1] (linear version), the unifiedscheme of Yu et al. [19], the single node method of Tao2Figure 8: Comparison of the rRSE (L-2 error) betweenthe BFL and ELIBB-Q methods simulated for thispaper and other results from the literature. The error isrelative to the cylindrical Couette region (blue region infigure 4). For the results of the Liu et al. method datafrom [14] is used. For the ZY methods data from [12] isused. Other parameters: TRT Λ = / , lattice topologyfor the cases simulated = D3Q27, τ = 2 . et al. [13] and the three variants of the ELIBB, plus theHWBB in the 10lu case. The simulated experiments werecarried out for Re = 10 in the case of the TRT collisionmodel using a D Q lattice layout and for q = 0 (seefigure 4). From these figures is possible to see that,1. fragmented methods BFL, ELIBB-FL and ELIBB-FQ generally more accurate than the unified meth-ods ELIBB-U and Tao;2. the space error distribution is similar for all themethods;3. the only local methods that performs similarly tothe BFL are the ELIBB-F and ELIBB-FQ;4. for low resolutions the present fragmented schemes(ELIBB-F and ELIBB-FQ) show smaller error thanthe BFL.
2. rRSE convergence
All the previously mentioned methods show approxi-mately second order convergence as displayed in figures 6and 7. The results appearing in the figures concern a lowreynolds fluid flow Re = 10 between cylinders with di-ameters ratio β = r a r b / . The choosen collision model isthe TRT with magic parameter Λ = / , viscous relax-ation numbers τ = 0 . and τ = 2 . Those figures also confirm the observations discussed above: the three frag-mented methods tested (BFL, ELIBB-FL, ELIBB-FQ)are generally slightly more accurate than the three uni-fied techniques (Yu, Tao, ELIBB-UL).We finally compare the results of the ELIBB-FQmethod with some results from the literature, namelythe Liu et al. [14] method and the ZY method [12]. TheLiu et al. can be seen as an extension of the Tao et al. ,while in the ZY method the single-node characteristic isachieved setting f ( f ∗ ¯ i ( x F F ) , t +1) ≈ f ( f ∗ ¯ i ( x F F ) , t ) . Evenif figure 8 cannot be considered a direct comparison, itsuggests that the ELIBB-FQ method results in a moreaccurate solution for the considered set of parameters.Note that in the steady-state test described by figure 8the ZY and BFL method should deliver the same resultsgiven that the ZY method is a first-order in time ap-proximation of the BFL. Therefore, the different resultsbetween the green and the red points in figure 8 are mostprobably caused by the different lattice topology, errorcomputation, numerical precision, or position relative tothe lattice of the cylinders axis.
3. rRSE as a function of q Some authors [34] suggested that unified methods canhelp to improve the stability of moving boundaries be-cause they do not need to change interpolation pointscrossing the midpoint between two nodes. We tried toassess this behavior, modifying the parameter q for dif-ferent resolutions keeping the other parameters constant.The results are shown in figure 9. From the results of thisspecific experiment, it is not possible to confirm this in-sight, since the distance q has a small impact on errorlevels for all methods but the HWBB. This kind of testhas also been recently carried out for BFL and Yu meth-ods in [17], for a Poiseuille flow with moving walls.
4. Non-equilibrium computation at the wall
Section II C 1 provided two ways of computing the non-equilibrium component at the wall position with a first-order approximation. In our experience, the two methodslead to similar results for low Re numbers. In figure 11,we present an example of a comparison of the two meth-ods. The figure shows that the technique used for non-equilibrium computation does not seem to have any im-pact on the error, at least for this range of parameters.However, there is evidence that regularized approachesare more stable for higher values of the Reynolds number.Hence, as future work, it will be interesting to furthercompare both approaches in the low viscosity regime.It is somehow important to mention that the non-equilibrium bounce-back (equation (17)) is generally eas-ier to implement than the regularized approach, becausethe latter can be used to impose both macroscopic val-ues and their gradient similarly to Robin boundary con-3 (a) (b)
Figure 9: Comparison of rRSE (36) for two different resolutions (a) r − r = 20 + 2 q , (b) r − r = 30 + 2 q as afunction of q for a set of methods. Parameters used for the simulation: Re = 10 , TRT collision model, Λ = / , τ ν = 1 . (a) (b) (c) Figure 10: Sensitivity to τ of the relative root squared error rRSE (36), i.e. L2 error function, in the case ofcylindrical Couette flow. Parameters of the simulation: Re = 10 , inter-cylinder distance L = 20 lu , β = 0 . . For theTRT collision model and the magic parameter Λ is fixed to / .ditions. This can be useful for advanced boundary mod-eling, like for turbulent wall modeling [78].
5. Evaluation of the viscosity dependence of the error
The sensitivity of the results to the viscosity relaxationtime τ , under diffusive scaling, has been investigated forBGK, TRT, and RRBGK collision models. In the case ofBGK (figure 10a) all models, and especially the HWBB,in accordance with the work by Ginzburg et al. [8, 9] Forthe lower values of τ , the beneficial displacement of the boundary caused by the viscosity dependence is higherthan the accuracy deterioration owed to the increasedtime-step (consequent to increased τ ). For higher valuesof τ , the error rise for all methods. This confirms a non-negligible viscosity dependence effect in these methods ifcoupled with the simple BGK. It is anyhow interestingto notice how, in this experiment and for higher valuesof τ , the simple half-way bounce-back method performsbetter than any interpolated version in the BGK case. Inthe case of the TRT collision model, the τ dependence ofthe HWBB is almost wiped out and its rRSE becomesalmost a horizontal line in figure 10. As a final comment4Figure 11: Comparison of the convergence of the rRSEfor the proposed ELIBB methods variants in the case ofthe two different non-equilibrium componentcomputation. The neqbb acronym stands for methoddescribed by equations (17) and the acronym reg standsfor the method described by equation (C8)).on figure 10, one may notice that the present ELIBBschemes show a resilient behavior for high values of τ :this fact combined with the good accuracy at low reso-lutions make the ELIBB particularly robust methods forcoarse space and time resolution simulations.
6. Mass conservation violation
A common issue with interpolated bounce-back isthe violation of mass conservation. We investigatedthis concern computing the average density fluctuation (cid:12)(cid:12)(cid:12)(cid:80) Ni ( ρ i − ρ ) /N (cid:12)(cid:12)(cid:12) in the Taylor-Couette region at non-dimensional time t ∗ = t/t ref = 3 . where t ref = ( r − r ) /u and N is the number of nodes in the Taylor-Couetteregion; the results are presented in figures 12 and 13. Theconclusion is that, even though the ELIBB are pretty ac-curate at low time and space resolutions, in this regimethey show a higher mass violation. Nonetheless, for finermeshes they exhibit similar values of the average densityfluctuation. C. Jeffery’s orbit: Ellipsoidal Cylinder andEllipsoid Rotation
Jeffery’s orbit is a common benchmark test for curvedboundary conditions. It describes the rotation of el-lipsoidal objects induced by a shear flow in Stoke’sregime [79]. Therefore, it is well suited to verify the ca-pability of a numerical method to describe a fluid-solidinteraction problem. In our experiment, the ellipsoid is Figure 12: Average density fluctuation | ρ LB − | in thecylindrical Couette flow region at t ∗ = t/t ref = 3 . where t ref = ( r − r ) /u . Parameters: TRT collisionmodel Λ = / , τ = 1 .Figure 13: Average density fluctuation | ρ LB − | in thecylindrical Couette flow region at at t ∗ = t/t ref = 3 . where t ref = ( r − r ) /u . Parameters: TRT collisionmodel Λ = / , τ = 2 .located at the center of a channel. The channel is de-limited in the y direction by two horizontal walls, thatare impulsively moved along the x direction at the be-ginning of the simulation with the fluid at rest. In theinitial condition, the prolated ellipsoid lies at the centerof the channel, in vertical position, with its longer di-ameter aligned with the y direction. The computationaldomain is periodic in the x and z directions. At the time t = 0 the upper wall is abruptly accelerated to its ter-minal velocity u lid , such that it generates the following5Figure 14: Jeffery’s orbit described by a 3D prolateellipsoid characterized by a ratio of radii β = r a / r b = 2 .Plots of the angular velocity obtained with differentmethods compared with the analytical solution.Parameters of the simulations τ = 2 , BGK collisionmodel, Re = 4 , H = 10 r e (channel height), r e = (cid:112) r a + 2 · r b = 30 (equivalent radius).strain rate in the channel ˙ γ = u lid H (38)where H is the channel height. The ellipsoid start ac-celerating until it reaches its steady state Jeffery’s orbit,that reads [79] ˙ θ = ˙ γr e + 1 (cid:0) r e cos θ + sin θ (cid:1) (39)where θ is the inclination of the ellipsoid axis correspond-ing to the major radius r a with respect to the y -axis (thevertical one), ˙ θ is the corresponding angular velocity and r e = (cid:112) r a + 2 · r b is the equivalent radius. The Reynoldsnumber in this scenario is redefined with the shear stress Re = ˙ γr ν = u lidr e / H ν . (40)In this experiment, we describe the ellipsoid with athin-shell surface whose dynamics are computed follow-ing the rigid body motion equations. The thin-shell rep-resenting the ellipsoid is filled with fluid. Nevertheless,the internal fluid is virtual and has no impact on the dy-namics. This because the forces from the fluid to the rigidbody are computed only considering the external fluid,using equation (29). All the investigated methods lead tosimilar values of the ellipsoid angular velocity evolutionin time in the case of a prolate ellipsoid with a diame-ter ratio β = r a / r b = 2 . In particular in figure 14, wecompared the results of a 3D ellipsoid angular velocity Figure 15: Jeffery’s orbit described by a 3D prolateellipsoid characterized by a ratio of radii β = r a / r b = 2 .Plots of the squared fluctuation of the torque obtainedwith equation (D3).Figure 16: Jeffery’s orbit described by a 3D prolateellipsoid characterized by a ratio of radii λ = r a / r b = 2 .Plots of the squared fluctuation of the torque obtainedwith equation (D3).evolution for the HWBB, BFL, and the present ELIBB-FQ, but the results are similar for all the ELIBB variantsand the other methods investigated in the previous sec-tion.
1. Numerical noise
When dealing with moving boundaries, LBM bound-ary methods give rise to spurious pressure oscillations.One of the main sources of these oscillations is the re-6sults of some nodes, named fresh nodes , changing theside of the boundary surface. In the upwind part of thesurface fresh nodes appear as uninitialized nodes withwrong populations’ values: this causes the triggering ofpressure waves. The process of recomputing the valuesfor the distribution functions in the fresh nodes is called refilling and can reduce the magnitude of pressure os-cillations. The refill methods have recently been com-pared in [76]. The authors indicate the local iterationrefill (LIR) [34] as the most effective in reducing oscil-lations. For this reason, we choose the LIR to correctthe spurious pressure oscillations, implementing it in aslightly modified version to make it consistent for thecase of thin-shell boundaries (see appendix D 2). In ourexperiments, using a thin-shell approach, the simulationsare stable even if no refilling algorithm is used. There-fore, we tested different boundary conditions before andafter the implementation of the LIR.The pressure oscillations due to the boundary motionare transferred to the rigid body through the momen-tum coupling, leading to a noisy torque time evolution.The measure of the oscillations in the torque (or resultingforce) acting on the rigid body, is a common way to esti-mate the magnitude of pressure oscillations [76, 80, 81].It is anyhow important to notice that perturbations in thetorque acting on the body are only an indirect measureof the effect of pressure waves and can be also influencedby the techniques used for the force computation.We decided to use a qualitative approach to computethe TSF that allows producing smooth graphs that areeasy to compare. The detail of the computations is shownin appendix D 1. Here, we only point out that the pre-sented values of the TSF are interpolated values usingbest fitting polynomials. Therefore, they should be in-terpreted as qualitative measures that do not aim at ac-curate measurements.Using the momentum exchange algorithm and withoutrefilling techniques, local methods show higher fluctua-tions (figure 15). Nevertheless, after the implementationof the LIR, both local and non-local interpolated meth-ods show similar performance in terms of torque squaredfluctuation (figure 16).In figure 17 we also computed the average torquesquared fluctuation
T SF on < θ (cid:46) π / for a ellipsoidalcylinder and plotted it for an increasing resolution of theellipsoids. The results show convergence with approxi-mate slope x − (this is expected for second-order schemessince we are working on the squared fluctuations).To summarize the results of this section, the ELIBBshows good stability propriety in the simulation of a fluid-rigid body problem interaction. The ELIBB leads to re-sults comparable with the BFL for the ellipsoid dynam-ics, both in terms of angular velocity evolution and interms of torque fluctuation. Nevertheless, for the localmethods tested the refill algorithm is important to reducethe effects on the torque due to the pressure oscillation. Figure 17: Average root squared fluctuation of thetorque for an ellipsoidal cylinder (with ratio betweenaxes d a / d b = 2 ) computed from the initial condition ofresting fluid, ω = 0 and ellipsoid in verical position( θ = 0 ), to θ ≈ π / . Local label stands for ELIBB andTao methods, LIR means local iteration refilling . IV. CONCLUSIONS
In this paper, we presented a novel class of enhancedsingle-node boundary conditions (ELIBB). The enhance-ment theoretically derives from the introduction of wallpopulations and the optimization of coefficients corre-lated with the quality of the interpolation in equa-tion (15). The resulting boundary methods are single-node , which means that they are suitable to simulatecomplex shapes with narrow gaps containing only one lat-tice node, without introducing special conditions. This physical locality feature of the ELIBB, facilitates themanagement of complex geometries. Thus, the ELIBBcan be implemented in algorithmic local way without anyneed of introducing special cases near corners or nar-row gaps. In parallel simulations no data will need tobe recovered from neighboring cells, improving in thisway parallel computing performances. Therefore, ELIBBboundary conditions are attractive for GPU based imple-mentations.The novel boundary conditions show, in the investi-gated experiments, a robust behavior in terms of accu-racy of the velocity field at low resolutions and largevalues of the relaxation parameters τ , and, in thisrange, they appear to be more accurate than the well-established Bouzidi-Firdaouss-Lallemand [1] boundarycondition. The ELIBB variants shares with other in-terpolated bounce-back methods some limitations: massviolation and viscosity dependence. In the performed ex-periment, mass violation of the novel methods is highercompared to the BFL at low resolutions. At higher res-olutions, this gap is bridged and the mass violation be-comes similar to the one of BFL. Further investigationsare needed to quantify these preliminary findings.7We showed that the method recently proposed by Tao et al. [13], can be interpreted as a variant of the ELIBBfamily. Nonetheless, the variants proposed in the presentpaper show a non-negligible improvement with respect tothe Tao et al. [13] method, at least for the low Reynoldsnumber regime.The novel class of boundary conditions is suitable todescribe moving boundaries immersed in the fluid. Inparticular, it proved to be stable and well-behaving whendescribing the dynamics of a rigid body in a shear Be-sides, local iteration refill algorithm has proved to be agood companion for the ELIBB. After its adoption toreduce pressure oscillations, single-node techniques andthe BFL method showed the same level of noise on thetorque acting on the rigid body.We additionally provided an alternative non-equilibrium estimation scheme that can be valuable, inperspective when a more advanced boundary modeling isneeded in place of a simple Dirichlet boundary condition( e.g. wall models for turbulent flows). This aspect needsfurther investigation in the future.From our first initial results, the ELIBB display en-couraging improvements over some earlier single-nodesboundary conditions, yet maintaining an easy implemen-tation characteristic. Therefore it is a good general-purpose candidate to replace the BFL in many applica-tions where a single-node method is required. Neverthe-less, more systematic studies are necessary to properlycompare it with other techniques, especially in the rangeof higher Reynolds numbers. ACKNOWLEDGMENTS
Fruitful discussions with Christophe Coreixas andChristos Kotsalos are gratefully acknowledged. Thisproject has received funding from the European Union’sHorizon 2020 research and innovation programme undergrant agreement No 823712 (CompBioMed2 project).
Appendix A: Hermite Polynomials
The Hermite functions are the solution solutions to theHermite ordinary differential equation [82]: H (cid:48)(cid:48) n ( x ) − x H (cid:48) n ( x ) + 2 n H n ( x ) = 0 . (A1)The functions H n are of integer degree n solution of (A1)are the physicist Hermite polynomials [83]. Their gen-eral expression can be evaluated using the Rodriguez for-mula [82] H n ( x ) = ( − n e x d n dx n e − x . (A2) or, equivalently, using the generating function [82, 84] g ( x, t ) = e − t +2 tx = ∞ (cid:88) n =0 H n ( x ) t n n ! (A3)In the LB context is more common to refer to the proba-bilistic Hermite polynomials, which are a rescaled versionof the H n ( x ) . In one dimension their expression is [83] H n ( x ) = 2 − n H n (cid:18) x √ (cid:19) = ( − n e x d n dx n e − x . (A4)For the polynomial expansion of the probability densityfunction in the velocity space it is necessary to performa multivariate extension of the classic Hermite polyno-mial [85]. In this case the Rodrigues’ formulas reads H n ( x ) = H α ...α n ( x ) = ( − n e x ∂ n ∂ α ...α n e − x . (A5)where x = ( x , . . . , x d ) , d is the dimensionality and α i ∈{ x , . . . , x d } . In the velocity space the expression of thefirst polynomials reads: H = 1 H α = ξ α H α α = ξ α ξ α − δ α α H α α α = ξ α ξ α ξ α − ( δ α α ξ α + δ α α ξ α + δ α α ξ α ) H α α α α = ξ α ξ α ξ α ξ α − ( δ α α ξ α ξ α + δ α α ξ α ξ α + δ α α ξ α ξ α + δ α α ξ α ξ α + δ α α ξ α ξ α + δ α α ξ α ξ α )+ ( δ α α δ α α + δ α α δ α α + δ α α δ α α ) (A6)In the discrete case, the Hermite polynomial needs tobe rescaled as a consequence of the rescaling of the dis-crete velocities c i = ξ i c s H i = 1 (A7) H i,α = c i,α (A8) H i,α α = c i,α c i,α − c s δ α α (A9) H i,α α α = c i,α c i,α c i,α (A10) − c s ( δ α α c i,α + δ α α c i,α + δ α α c i,α ) (A11) H i,α α α α = c i,α c i,α c i,α c i,α (A12) − c s ( δ α α c i,α c i,α + δ α α c i,α c i,α (A13) + δ α α c i,α c i,α + δ α α c i,α c i,α (A14) + δ α α c i,α c i,α + δ α α c i,α c i,α ) (A15) + c s ( δ α α δ α α + δ α α δ α α + δ α α δ α α ) (A16)where H i,α ..αn = H α ..αn ( c i ) and consequently H i,n = H n ( c i ) .8 Appendix B: Regularization procedures
The concept of regularized populations was originallyintroduced by Skordos [86] and further developed byLatt, Malaspinas et al. [55, 56, 87–89]. The cornerstoneof this approach is to reconstruct populations keepingonly the minimal information to ensure that the cor-rect macroscopic behavior is recovered. The purposeis to discard the information contained in the popula-tion that does not refer to the Navier-Stokes physics,to filter numerical errors. To this effect, populationsare reconstructed from a truncated Hermite polynomialexpansion of f (0) i and f (1) i at order M and N respec-tively [42, 43, 85]. It is possible to exemplify the regular-ized LB considering the following steps:1. Consider the initial condition, at iteration t − ,in which the non-equilibrium populations consistssolely of hydrodynamic components f neq i = f (1) i ∼O (Kn f eq ) with null high order Hermite polyno-mials components. Where higher order means –not necessary to recover the Navier-Stokes levelphysics–. We have, f i = f (0) i,M + f (1) i,N , (B1a) f (0) i = f (0) i,M , (B1b) f (1) i = f (1) i,N , (B1c)where M and N are the polynomial orders. And,after collision, the populations reads f ∗ i = f (0) i,M + (cid:18) − τ (cid:19) f (1) i,N (B2)where f ∗ i denotes a post-collision population.2. After streaming, at time t , non-hydrodynamic errorcomponents ε i caused by the numerical discretiza-tion appear ˜ f i = ˜ f eq (cid:122) (cid:125)(cid:124) (cid:123) f (0) i,M (cid:8)(cid:8)(cid:8) + ε eq i (cid:124) (cid:123)(cid:122) (cid:125) to be filtered + ˜ f neq (cid:122) (cid:125)(cid:124) (cid:123) f (1) i,N (cid:8)(cid:8)(cid:8) + ε neq i (cid:124) (cid:123)(cid:122) (cid:125) to be filtered (B3)where ˜ f i , ˜ f eq i , and ˜ f neq i are respectively the non-filtered population, equilibrium component andnon-equilibrium component.3. Projective regularization step: it is a filtering pro-cess consisting in estimating f (0) i,M and f (1) i,N from theHermite moments of ˜ f eq i and ˜ f neq i discarding highorder moments contributions. The projective regu- larization step reads f (1) i,N ≈ f neq i,N = w i N (cid:88) n =0 n ! c ns ˜ a neq n : H i,n (B4) ˜ a neq n ≈ Q − (cid:88) i =0 w i (cid:16) ˜ f i − f eq i,M (cid:17) H i,n (B5) f eq i,M = f (0) i,M = w i M (cid:88) n =0 n ! c ns a (0) n : H i,n (B6)where H i,n = H n ( c i = c s ξ i ) is a Hermite polyno-mial of order n (see appendix A), a n is its corre-sponding tensor of coefficients, “:” is the Frobeniusinner product [90], w i are the quadrature weightsand c s is the lattice constant.4. finally, filtered population f i is reconstructed as f i ≈ f (0) i,M + f neq i,N . (B7)Even if the RBGK was successful, it does not filtersout all the non-hydrodynamic components, because the f neq i computed from the projection still contains numer-ical errors due to higher perturbation order (in Knud-sen number) components ˜ f (2) i + ˜ f (3) i + . . . . For thisreason, it has been extended by the recursive regular-ized model (RRBGK). In this case, the Hermite coeffi-cients appearing in equations (B5) are recomputed usingthe recursive formulas derived in references [56, 88, 89](see equations (C5) and (C6) in appendix C). The re-cursive formulation improves the filtering of numericalnon-hydrodynamics errors for the Hermite componentsof order higher than two. Appendix C: Recursive estimation of thenon-equilibrium populations
In this section, we present the formal techniques usedto enhance the boundary condition with an estimationof the wall populations. To achieve this, it is necessaryto formally unfold the relation that elapses between themesoscopic representation of the Navier-Stokes equationsand the mesoscopic perspective of the BE. Although di-verse mathematical approaches exist, the most commonin the LB community is the Chapman-Enskog perturba-tive expansion (CE) [61]. This method has been extendedwith the Grad-Hermite (GH) expansion [85] leading to anelegant and systematic approach recovering the hydrody-namic limit of both the BE and LBE [43, 56, 87–89].In the Chapman-Enskog-Hermite expansion [56, 61,88, 89] the populations and the LBE are expanded ina perturbative formulation and decomposed in a set of9Hermite-basis moments equations. The perturbative ex-pansion of the velocity distribution function is givenby equation (10). The equilibrium f (0) and first non-equilibrium component f (1) components of equation (10)can be projected in the Hermite polynomial basis f eq i,M = f (0) i,M = w i M (cid:88) n =0 n ! c ns a (0) n : H i,n (C1) f (1) i,N = w i N (cid:88) n =0 n ! c ns a (1) n : H i,n (C2)where a (0) n and a (1) n are n -th order Hermite tensor coef-ficients of respectively the zero and the first order of theCE expansion (equation (10)) a (0) n = Q − (cid:88) i =0 w i f (0) i H i,n (C3) a (1) n = Q − (cid:88) i =0 w i f (1) i H i,n . (C4)Performing the CE of the LBE give rise an hierarchy ofequations. Taking the moments of the first two equa-tions (order zero and one of the CE) and comparingthem with the Navier-Stokes equations, one can derivethe Malaspinas’ recursive formulas [56, 88, 89], that inthe isothermal case read respectively ρa (0) α ..α n = a (0) α n a (0) α ..α n − ∀ n ≥ , (C5)and ρa (1) α ..α n = a (0) α n a (1) α ..α n − + n − (cid:88) l =1 a (0) α ..α l − α l +1 ..α n − a (1) α l α n (C6) ∀ n ≥ . The comparison between the Hermite momentsof the Chapman-Enskog expanded LBE and the Navier-Stokes equation allows also to bind the relaxation timeand the fluid viscosity through the following relation a (1) α α = τ ρ Λ α α (C7)where Λ α α = ∂ α u α + ∂ α u α is the macroscopicisothermal-incompressible stress tensor and D is theproblem dimensionality (one, two or three).In the RRBGK the recursive formulas are used inequations (C1) and (C2) to rebuild the populations upto the first oreder of the CE expansion, filtering non-hydrodynamic components. To correctly recover a Her-mite moment of order n at the Knunsen order k , M should be greater or equal to n + k +1 and N ≥ n + k [43].Thus, even in the isothermal case, the optimal order oftruncation for the equilibrium is the fourth. Note that,even if in the case of isothermal lattices only momentsup to the second polynomial order are correctly recov-ered, equation (C5) allows to recompute a (0) > from a (0)2 . Regarding the order of truncation N , the minimal valuein the isothermal case is N = 3 , but there is evidencethat the optimal value is N = M = 4 [91]. Here weshow the procedure for the order three truncation forsake of conciseness. Using the Malaspinas’ formula (C6)in (C2) truncated at the third order and recalling thatequations (C7) one gets f (1) i, ≈ − w i τ ρc s (cid:20) c s H i,α α Λ α α + 16 c s H i,α α α (Λ α α u α + Λ α α u α + Λ α α u α ) (cid:21) , (C8)where we used the Einstein notation, therefore summa-tion is implicit for repeated indexes, α i ∈ { x , . . . , x d } , d is the dimensionality of the velocity space and x , . . . , x d are the names of the space axes.Thanks to the symmetry proprieties of the Hermitepolynomials, from the previous equation follows the sym-metry of non-equilibrium populations f neq i ( x α , t α ) = f neq¯ i ( x α , t α ) . (C9)The latter formula is the starting point of the wet-nodeZou-He [24, 92] boundary condition, also know as non-equilibrium bounce-back. Some authors [13, 14] usethis idea to compute off-lattice unknown non-equilibriumcomponents in an approximate way f neq i ( x α , t α ) ≈ f neq¯ i ( x β ≈ x α + δ x x α , t β ≈ t α + δ t t α ) (C10)where x α is the coordinate of the off-lattice position of in-terest, x β is the coordinate of a lattice node, δ x is a smallspace distance, δ t denotes a small time interval and t β isa timestep that can differ from t α . The approximatednon-equilibrium bounce-back is a first-order approxima-tion for the non-equilibrium component, nevertheless asshown by Chun et al. [11] the non-equilibrium componentis a second-order correction over the equilibrium compo-nent. Thanks to this fact, equation (C10) can still beused to develop a second-order accurate boundary con-dition.To be more general, in this paper we propose an addi-tional approach to explicitly estimate the off-lattice val-ues of f (1) i . Namely, using the expression (C8) to ap-proximate the non-equilibrium component at the walllocation. In equations (C8) all the variables appearingshould be interpreted in our case as computed at walllocation. For the stresses we use a first-order approxi-mation Λ α α ( x W ) ≈ Λ α α ( x F ) . Adopting (C8) can re-quire a slightly higher computational cost, but it allowsfor more flexibility on the boundary modeling if needed.For example, when is necessary to model Λ α α ( x W ) asin references [78, 93].0 Fresh Nodes?Yeslocal C&S on a local sub-lattice copyupdate fresh nodes on the original latticeconditionYesEND LI NoNoCollide & SreamBEGIN LIFresh Nodes detection
Figure 18: The modified LIR algorithm.
Appendix D: Technical details of Jeffery’s orbittest-case1. Qualitative approach to evaluate torqueoscillations
The theoretical value of the torque acting on the el-lipsoid (or ellipsoidal cylinder) changes over time. If theellipsoid inertia is small, the magnitude of the torque os-cillation exceeds the value of the theoretical torque. Un-fortunately, the analytical value of the torque is not avail-able and this makes it impossible to compute the oscilla-tion because the baseline solution is unknown. To give aqualitative representation of the evolution of the torque squared fluctuation, we decided to compute a smoothnumerical baseline solution using interpolating polyno-mials. This baseline regular numerical solution is thenused to compute the squared fluctuation of the torque(TSF) over time. The TSF of the torque computed inthis way is noisy and not enough precise to obtain quan-titative results. Nevertheless, we can use it to perform asecond interpolation with a polynomial of the same orderto get a qualitative estimation of the TSF evolution overtime and use it to compare visually different methods.In detail, the first polynomial least squares interpola-tion is of order n = 30 and reads T ip,n ( t j ) = argmin k (cid:88) j =0 | p n ( t j ) − T j ( t j ) | . (D1)After the computation of the torque squared fluctuation(TSF) T SF = ( T − T ip ,n ) , (D2)a second interpolating polynomial of the same order iscomputed for the fluctuations T SF ip,n ( t j ) = argmin k (cid:88) j =0 | p n ( t j ) − T ip ,n,j ( t j ) | (D3)where t j is the j − th iteration. This procedure allowsto produce regular and qualitative graphs of the torquesquared fluctuations in time.
2. Local Iteration Refilling
The Local Iteration refilling (LIR) proposed in [34] ismodified to make it consistent with a thin shell two-dimensional representation of the boundary. To do so,the LIR is applied before the global collide and stream,and not after . This detail guarantees that the streamingstep does not move wrong populations outside the freshnodes before the application of the LIR. The modifiedalgorithm is represented in figure 18. [1] M. Bouzidi, M. Firdaouss, and P. Lallemand, Phys. Flu-ids , 3452 (2001).[2] Y. Thorimbert, F. Marson, A. Parmigiani, B. Chopard,and J. Latt, Computers & Fluids , 286 (2018).[3] C. Huber, A. Parmigiani, J. Latt, and J. Dufek, WaterResour. Res. , 6371 (2013).[4] C. Kotsalos, J. Latt, and B. Chopard, J. Comput. Phys. , 108905 (2019).[5] S. Leclaire, A. Parmigiani, O. Malaspinas, B. Chopard,and J. Latt, Physical Review E , 033306 (2017).[6] A. J. C. Ladd, J. Fluid Mech. , 285 (1994).[7] I. Ginzbourg and D. d’Humières, J. Stat. Phys. , 927(1996). [8] I. Ginzburg and D. d’Humières, Physical Review E ,066614 (2003).[9] Kruger, Timm and Kusumaatmaja, Halim and Kuzmin,Alexandr and Shardt, Orest and Silva, Goncalo, The lat-tice Boltzmann method (a, New York, NY, 2016).[10] M. Junk and Z. Yang, Physical Review E , 066701(2005).[11] B. Chun and A. J. C. Ladd, Physical Review E ,066705 (2007).[12] W. Zhao and W.-A. Yong, J. Comput. Phys. , 1(2017).[13] S. Tao, Q. He, B. Chen, X. Yang, and S. Huang, Comput-ers & Mathematics with Applications , 1593 (2018). [14] S. Liu, T. Zhou, S. Tao, Z. Wu, and G. Yang, Interna-tional Journal of Modern Physics C , 1950041 (2019).[15] G. Silva and I. Ginzburg, Philosophical Transactions ofthe Royal Society A: Mathematical, Physical and Engi-neering Sciences , 20190404 (2020).[16] C. Peng, O. M. Ayala, and L.-P. Wang, Computers &Fluids 10.1016/j.compfluid.2019.06.032 (2019).[17] P. Zhang, S. A. Galindo-Torres, H. Tang, G. Jin, andL. Li, content , 24.[18] A. De Rosis, S. Ubertini, and F. Ubertini, J. Sci. Comput. , 477 (2014).[19] D. Yu, R. Mei, and W. Shyy, in (American Institute of Aeronauticsand Astronautics, 2003).[20] W. Zhao, J. Huang, and W.-A. Yong, Multiscale Model-ing & Simulation , 854 (2019).[21] L. Wang, S. Tao, X. Meng, K. Zhang, and G. Lu, PhysicalReview E , 063307 (2020).[22] X. Meng, L. Wang, W. Zhao, and X. Yang, Advances inWater Resources , 103583 (2020).[23] C. Peng, O. M. Ayala, J. C. Brändle de Motta, and L.-P.Wang, Computers & Fluids , 104251 (2019).[24] Q. Zou and X. He, Phys. Fluids , 1591 (1997).[25] O. Filippova and D. Hänel, J. Comput. Phys. , 219(1998).[26] O. Filippova and D. Hänel, Computers & Fluids , 697(1997).[27] Succi, S., The Lattice Boltzmann equation: for complexstates of flowing matter , First edition ed. (Oxford Uni-versity Press, Oxford, 2018).[28] Guo, Zhaoli and Shu, Chang,
Lattice Boltzmann methodand its applications in engineering , Advances in compu-tational fluid dynamics (World Scientfic, 2013).[29] R. Mei, L.-S. Luo, and W. Shyy, J. Comput. Phys. ,307 (1999).[30] Z. Yang,
Analysis of Lattice Boltzmann Boundary Con-ditions , Ph.D. thesis (2007).[31] R. W. Nash, H. B. Carver, M. O. Bernabeu, J. Hether-ington, D. Groen, T. Krüger, and P. V. Coveney, PhysicalReview E , 023303 (2014).[32] D. R. Noble and J. R. Torczynski, International Jour-nal of Modern Physics C 10.1142/S0129183198001084(1998).[33] D. R. Noble and J. R. Torczynski, International Jour-nal of Modern Physics C 10.1142/S0129183198001084(2011).[34] L. Chen, Y. Yu, J. Lu, and G. Hou, Int. J. Numer. Meth-ods Fluids , 439 (2014).[35] J. Latt, O. Malaspinas, D. Kontaxakis, A. Parmigiani,D. Lagrava, F. Brogi, M. B. Belgacem, Y. Thorimbert,S. Leclaire, S. Li, F. Marson, J. Lemus, C. Kotsalos,R. Conradin, C. Coreixas, R. Petkantchin, F. Raynaud,J. Beny, and B. Chopard, Computers & Mathematicswith Applications 10.1016/j.camwa.2020.03.022 (2020).[36] Cercignani, Carlo, The Boltzmann Equation and Its Ap-plications , edited by John, F. and Marsden, J. E. andSirovich, L., Applied Mathematical Sciences, Vol. 67(Springer New York, New York, NY, 1988).[37] Cercignani, Carlo and Illner, Reinhard and Pulvirenti,M,
The mathematical theory of dilute gases (1994).[38] C. Cercignani and E. Gabetta, eds.,
Transport phenom-ena and kinetic theory: applications to gases, semi-conductors, photons, and biological systems , Modelingand simulation in science, engineering and technology (Birkhauser, Boston, 2007).[39] G. M. Kremer,
An introduction to the Boltzmann equa-tion and transport processes in gases , Interaction of me-chanics and mathematics (Springer, Berlin ; New York,2010).[40] P. L. Bhatnagar, E. P. Gross, and M. Krook, Phys. Rev. , 511 (1954).[41] P. Welander, Arkiv Fysik Vol: 7 (1954).[42] X. Shan and X. He, Phys. Rev. Lett. , 65 (1998).[43] X. Shan, X.-F. Yuan, and H. Chen, J. Fluid Mech. ,413 (2006).[44] P. J. Dellar, Comput. Math. Appl. , 129 (2013).[45] Y. H. Qian, D. D\textquotesingleHumiÚres, andP. Lallemand, Europhysics Letters (EPL) , 479 (1992).[46] C. Coreixas, G. Wissocq, B. Chopard, and J. Latt, Phil.Trans. R. Soc. A , 20190397 (2020).[47] C. Coreixas, B. Chopard, and J. Latt, Physical ReviewE , 033305 (2019).[48] D. d’Humieres, Rarefied gas dynamics (1992).[49] P. Lallemand and L.-S. Luo, Phys. Rev. E , 6546(2000).[50] D. d’Humières, I. Ginzburg, M. Krafczyk, P. Lallemand,and L.-S. Luo, Philosophical Transactions of the RoyalSociety of London. Series A: Mathematical, Physical andEngineering Sciences , 437 (2002).[51] I. Ginzburg and F. Verhaeghe, Commun. Comput. Phys., 53 (2008).[52] I. Ginzburg and F. Verhaeghe, Commun. Comput. Phys., 64 (2008).[53] D. d’Humières and I. Ginzburg, Comput. Math. Appl. , 823 (2009).[54] I. Ginzburg, D. d’Humières, and A. Kuzmin, J. Stat.Phys. , 1090 (2010).[55] J. Latt and B. Chopard, Math. Comput. Simul DiscreteSimulation of Fluid Dynamics in Complex Systems, ,165 (2006).[56] O. Malaspinas, Arxiv (2015).[57] C. Coreixas, G. Wissocq, G. Puigt, J.-F. Boussuge, andP. Sagaut, Physical Review E , 033306 (2017).[58] J. Jacob, O. Malaspinas, and P. Sagaut, J. Turb. , 1051(2018).[59] G. Wissocq, C. Coreixas, and J. Boussuge, arXivpreprint arXiv:2006.07353 (2020), arXiv:2006.07353[physics.comp-ph].[60] G. Silva and V. Semiao, Phys. Rev. E , 013311 (2017).[61] S. Chapman, T. Cowling, D. Burnett, and C. Cercig-nani, The Mathematical Theory of Non-uniform Gases:An Account of the Kinetic Theory of Viscosity, ThermalConduction and Diffusion in Gases , Cambridge Mathe-matical Library (Cambridge University Press, 1953).[62] I. Ginzbourg and P. M. Adler, Journal de Physique II ,191 (1994).[63] I. Ginzburg and F. Verhaeghe, Commun. Comput. Phys., 53 (2008).[64] J. L. de Lagrange, Leçons élémentaires sur les Mathéma-tiques, données à l’École normale, en 1795 (1812).[65] J. Sylvester, The London, Edinburgh, and Dublin Philo-sophical Magazine and Journal of Science , 267 (1883).[66] J. Latt, B. Chopard, O. Malaspinas, M. Deville, andA. Michler, Physical Review E , 056703 (2008).[67] O. Malaspinas, B. Chopard, and J. Latt, Computers &Fluids , 29 (2011).[68] A. J. C. Ladd, in Computer Simulation in Materials Sci-ence: Interatomic Potentials, Simulation Techniques and Applications , NATO ASI Series, edited by M. Meyer andV. Pontikis (Springer Netherlands, Dordrecht, 1991) pp.481–504.[69] T. Inamuro, K. Maeba, and F. Ogino, International Jour-nal of Multiphase Flow , 1981 (2000).[70] R. Mei, D. Yu, W. Shyy, and L.-S. Luo, Physical ReviewE , 041203 (2002).[71] Y. Hu, D. Li, S. Shu, and X. Niu, Physical Review E ,033301 (2015).[72] E.-J. Ding and C. K. Aidun, J. Stat. Phys. , 685(2003).[73] Y. Chen, Q. Cai, Z. Xia, M. Wang, and S. Chen, PhysicalReview E , 013303 (2013).[74] H. Li, X. Lu, H. Fang, and Y. Qian, Physical Review E , 026701 (2004).[75] B. Wen, C. Zhang, Y. Tu, C. Wang, and H. Fang, J.Comput. Phys. , 161 (2014).[76] S. Tao, J. Hu, and Z. Guo, Computers & Fluids , 1(2016).[77] M. Emin Erdoˇgan, International Journal of Non-LinearMechanics , 1091 (2002).[78] O. Malaspinas and P. Sagaut, Journal of ComputationalPhysics , 25 (2014).[79] G. B. Jeffery and L. N. G. Filon, Proceedings of theRoyal Society of London. Series A, Containing Papers ofa Mathematical and Physical Character , 161 (1922).[80] A. Caiazzo and M. Junk, Computers & Mathematics withApplications , 1415 (2008).[81] X. Li, F. Jiang, and C. Hu, Computers & Fluids , 33(2016). [82] G. B. Arfken, H. J. Weber, and F. E. Harris, in Math-ematical Methods for Physicists (Seventh Edition) (Aca-demic Press, Boston, 2013) pp. 871–933.[83] M. Abramowitz and I. Stegun,
Handbook of MathematicalFunctions: With Formulas, Graphs, and MathematicalTables , Applied mathematics series (Dover Publications,1965).[84] G. B. Arfken, H. J. Weber, and F. E. Harris, in
Mathe-matical Methods for Physicists (Seventh Edition) , editedby G. B. Arfken, H. J. Weber, and F. E. Harris (AcademicPress, Boston, 2013) pp. 551–598.[85] H. Grad, Commun. Pure Appl. Math. , 331 (1949).[86] P. A. Skordos, Phys. Rev. E , 4823 (1993).[87] J. Latt, Hydrodynamic Limit of Lattice Boltzmann Equa-tions , Ph.D. thesis, University of Geneva (2007).[88] Malaspinas, Orestis Pileas,
Lattice Boltzmann method forthe simulation of viscoelastic fluid flows , Ph.D. thesis(2009).[89] O. Malaspinas and P. Sagaut, J. Fluid Mech. , 514(2012).[90] A. R. Amir-Moéz and C. Davis, Mathematische Annalen , 107 (1960).[91] C. Coreixas,
High-order extension of the recursive reg-ularized lattice Boltzmann method , Ph.D. thesis, INPToulouse (2018).[92] X. He and Q. Zou, (1995), arxiv:comp-gas/9507002[comp-gas].[93] G. I. Park and P. Moin, Physics of Fluids26