Lincheng Xu
email: lincheng.xu@univ-amu.fr.

Eric Serre

Pierre Sagaut

A theoretical analysis of mass leakage at boundaries within the lattice Boltzmann method

Mass leakage at boundaries can be a critical issue for reliability of the lattice Boltzmann (LB) method based on Cartesian grids. Despite numerous work based on the LB method, the intrinsic macroscopic mechanisms causing mass leakage are still not fully characterised, but are essential to improve the mass conservation of LB simulations. In this paper, an original theoretical investigation of mass leakage at boundaries is proposed within the general LB framework. It is demonstrated that the mass leakage originates from the intrinsic deficiency of the wall-cut LB links at boundary nodes in recovering macroscopic momenta. From a mesoscopic-level definition, i.e. the net loss of distribution functions during the streaming process, the local mass leakage at individual boundary nodes and its averaged value along smooth boundaries are mathematically expressed using macroscopic variables. The local mass leakage is shown to be dominated by terms proportional to the tangential momentum component. In contrast, the averaged mass leakage is shown to be contributed from various terms including the boundary curvature, the tangential momentum, and the gradients of density, momentum and momentum flux. Meanwhile, amplitude of the averaged mass leakage is theoretically estimated to be proportional to the local grid spacing, based on which a first-order accurate correction scheme is proposed. In addition, both the local and averaged mass leakage are demonstrated to be significantly dependent on boundary orientation with respect to the grid. The proposed theoretical analysis is assessed by performing numerical experiments. Two-dimensional weakly compressible flows through straight and curved moving channels are considered to estimate each term appearing in the theoretical analysis. The numerical results are in very good agreement with the proposed analysis, and the proposed mass correction scheme based on the averaged mass leakage effectively cures the mass leakage problems in the considered test cases.

I. INTRODUCTION

During the last three decades, the lattice Boltzmann (LB) method has received a growing interest as an alternative computational fluid dynamic (CFD) tool, and it is nowadays widely applied in various problems [START_REF] Aidun | Lattice-Boltzmann method for complex flows[END_REF][START_REF] Zarghami | Assessment of interaction potential in simulating nonisothermal multiphase systems by means of lattice Boltzmann modeling[END_REF][START_REF] Xu | An immersed boundary-lattice Boltzmann method for gaseous slip flow[END_REF] . Compared to the conventional CFD techniques the LB method is praised for its linear convective term, simple parallel implementation and high computational efficiency for unsteady flow simulations. However, since the LB method is usually implemented on Cartesian grids, the stair-wise Cartesian grid boundaries (defined by "boundary nodes") near curved and/or inclined solid walls raise challenges for accurate boundary treatments. Especially, the mass conservation (or inversely the mass leakage minimization) at boundary nodes is of fundamental importance to guarantee reliable and accurate solutions, in particular in internal flows, but it is still not fully understood and requires to be theoretically addressed.

Despite the stair-wise distribution of boundary nodes, various boundary treatments have been developed within the LB framework. The simplest and most commonly used boundary method is the bounce-back (BB) scheme in which the distribution functions streaming towards boundaries are directly returned to the boundaries node along the reversed directions to mimic particles collision dynamics on non-slip interfaces [START_REF] Ziegler | Boundary conditions for lattice Boltzmann simulations[END_REF][START_REF] Ginzbourg | Boundary flow condition analysis for the three-dimensional lattice Boltzmann model[END_REF][START_REF] Peng | A comparative study of immersed boundary method and interpolated bounce-back scheme for no-slip boundary treatment in the lattice Boltzmann method: Part i, laminar flows[END_REF] . Although this method is simple to implement for arbitrary geometries and satisfies the mass conservation well, it has a zero-order accuracy for pressure and first-order accuracy for velocity only when non-aligned boundaries involved [START_REF] Ginzbourg | Local second-order boundary methods for lattice Boltzmann models[END_REF][START_REF] Junk | One-point boundary condition for the lattice Boltzmann method[END_REF][START_REF] Ginzburg | Spurious interface and boundary behaviour beyond physical solutions in lattice Boltzmann schemes[END_REF] , while the LB method itself is of second-order accurate [START_REF] Reider | Accuracy of discrete-velocity BGK models for the simulation of the incompressible Navier-Stokes equations[END_REF][START_REF] Oulaid | On the origin of numerical errors in the bounce-back boundary treatment of the lattice Boltzmann method: A remedy for artificial boundary slip and mass leakage[END_REF] . To develop more accurate boundary treatments, Chen et al. [START_REF] Chen | On boundary conditions in lattice Boltzmann methods[END_REF] proposed an extrapolation scheme to reconstruct the unknown distribution functions directly from their neighbouring counterparts. Filipova and Hanel [START_REF] Filippova | Grid refinement for lattice-BGK models[END_REF] proposed a scheme to reconstruct the unknown distribution functions from macroscopic variables and their gradients, which are extrapolated from those known in the fluid domain. These two schemes are expected to exhibit a second-order accuracy, but suffer from numerical instability in some situations due to the extrapolation operations. The extrapolation-caused instability was cured by Mei et al. [START_REF] Mei | An accurate curved boundary treatment in the lattice Boltzmann method[END_REF] by adjusting the extrapolation stencil properly. Later, Bouzidi et al. [START_REF] Bouzidi | Momentum transfer of a Boltzmann-lattice fluid with boundaries[END_REF] proposed a second-order accurate method which avoids extrapolation by combining the bounce-back concept with interpolation. This is achieved by extending the interpolation stencil for a prescribed LB link reversely to the fluid domain when a node in the non-fluid region is required. The scheme was further improved by Ginzburg and d'Humieres [START_REF] Ginzburg | Multireflection boundary conditions for lattice Boltzmann models[END_REF] who proposed a more general framework from which a third-order accurate multi-reflection scheme was derived. These schemes have been demonstrated to be ap-parently more accurate than the BB method. However, except in some special cases (e.g. simple steady Stokes flow simulations where the LB process is purely linear [START_REF] Ginzbourg | Local second-order boundary methods for lattice Boltzmann models[END_REF][START_REF] Ginzburg | Steady-state two-relaxation-time lattice Boltzmann formulation for transport and flow, closed with the compact multi-reflection boundary and interface-conjugate schemes[END_REF][START_REF] Ginzburg | Mass-balance and locality versus accuracy with the new boundary and interface-conjugate approaches in advection-diffusion lattice Boltzmann method[END_REF] ), the mass conservation is generally compromised by the applied interpolation or extrapolation operations [START_REF] Lallemand | Lattice Boltzmann method for moving boundaries[END_REF][START_REF] Krüger | The Lattice Boltzmann Method -Principles and Practice[END_REF] .

Apart from the above mentioned schemes focusing on the reconstruction of the unknown distribution functions, there is another kind of boundary treatment relying on the use of an external body force tuned to enforce some constrains on the macroscopic quantities. Following the seminal work on immersed boundary method by Peskin [START_REF] Peskin | Flow patterns around heart valves: a numerical method[END_REF] , many variants and improvements have been proposed and extensively applied within the LB framework, e.g. [START_REF] Xu | A novel geometry-adaptive Cartesian grid based immersed boundary-lattice Boltzmann method for fluid-structure interactions at moderate and high Reynolds numbers[END_REF][START_REF] Wang | Recent progress of immersed boundary method and its applications in compressible fluid flow[END_REF] . Usually, this boundary strategy uses lattice nodes in non-fluid regions so that the streaming process can be implemented without sensing the boundaries. On the one hand this feature is very favourable for moving and/or deformable boundary interfaces, but on the other hand it allows for local non-physical information leakage, e.g. mass and momentum, to the non-fluid regions.

Many research efforts have been devoted to studying the mass leakage in LB simulations. Conceptually, the LB method relies on the preservation of conservativity, i.e. preserving the zeroth and first-order moments of the distribution function to recover density and momentum, along with the associated collision invariants to recover the correct macroscopic mass and momentum conservation equations.

Ginzburg & d'Humières [START_REF] Ginzbourg | Local second-order boundary methods for lattice Boltzmann models[END_REF] proposed a local second order accurate boundary (LSOB) treatment reconstructing unknown distribution functions from macroscopic variables and their gradients according to Chapman-Enskog analysis up to second order. In the analysis of simple Stokes flows (Reynolds number much less than unit), including the Couette and Poiseuille flow, using linear LB simulations, the LSOB method gave exact results without mass leakage. Meanwhile, the local mass leakage along aligned boundaries was analysed to be linearly dependant on the normal momentum as well its high order (≥2) normal gradients. In a more general background beyond Stokes flow, Lallemand and Luo [START_REF] Lallemand | Lattice Boltzmann method for moving boundaries[END_REF] pointed out that interpolation during the reconstruction of unknown distribution functions could lead to a breakdown of mass conservation, thus yielding mass leakage. To avoid interpolation-related issues, Kao and Yang [START_REF] Kao | An investigation into curved and moving boundary treatments in the lattice Boltzmann method[END_REF] proposed an interpolation-free scheme based on the grid refinement concept [START_REF] Filippova | Grid refinement for lattice-BGK models[END_REF] without refining the grid in practice. However, this scheme still does not guarantee mass conservation [START_REF] Sanjeevi | Choice of no-slip curved boundary condition for lattice Boltzmann simulations of high-reynolds-number flows[END_REF] . Rohde et al. [START_REF] Rohde | Volumetric method for calculating the flow around moving objects in lattice-Boltzmann schemes[END_REF] observed mass leakage in LB simulations applying finite-volume flux techniques at moving boundaries. Bao et al. [START_REF] Bao | A mass conserving boundary condition for the lattice Boltzmann equation method[END_REF] observed apparent mass leakage in simulating flow through a U-shaped tube using interpolationbased boundary methods, and proposed a scheme to enforce mass conservation by directly ad- to the method proposed by Filipova and Hanel [START_REF] Filippova | Grid refinement for lattice-BGK models[END_REF] , "MLS" refers to the boundary treatment proposed by Mei et al. [START_REF] Mei | An accurate curved boundary treatment in the lattice Boltzmann method[END_REF] , "IPF" refers to the interpolation-free method proposed by Kao and Yang [START_REF] Kao | An investigation into curved and moving boundary treatments in the lattice Boltzmann method[END_REF] , "NEE" is short for "Non-equilibrium extrapolation", "Bouzidi" refers to the schemes proposed by Bouzidi et al. [START_REF] Bouzidi | Momentum transfer of a Boltzmann-lattice fluid with boundaries[END_REF] , "ZY" refers to the scheme proposed by Zhao and Yong [START_REF] Zhao | Single-node second-order boundary schemes for the lattice Boltzmann method[END_REF] , "LSOB" refers to the scheme proposed by Ginzburg 7 , and "MR" refers to the scheme proposed by Ginzburg and d'Humieres [START_REF] Ginzburg | Multireflection boundary conditions for lattice Boltzmann models[END_REF] .

Ref. Boundary methods Configurations

Bao et al. justing local density values. Similar mass correction schemes are also proposed and tested by other researchers [START_REF] Coupanec | A mass conserving boundary condition for the lattice Boltzmann method for tangentially moving walls[END_REF][START_REF] Yu | Modified curved boundary scheme for two-phase lattice Boltzmann simulations[END_REF] . There are also some mass-conservative strategies specifically valid for flat walls aligned with coordinate axes [START_REF] Oulaid | On the origin of numerical errors in the bounce-back boundary treatment of the lattice Boltzmann method: A remedy for artificial boundary slip and mass leakage[END_REF][START_REF] Coupanec | A mass conserving boundary condition for the lattice Boltzmann method for tangentially moving walls[END_REF][START_REF] Feng | Mass-conserved wall treatment of the non-equilibrium extrapolation boundary condition in lattice Boltzmann method[END_REF] . Especially, Feng and Lim 30 reported that mass leakage around sharp boundary corners is much more significant than along aligned walls in their LB simulations. In addition, it was reported that mass leakage in LB simulations of Stokes flow is so low that local mass correction is not necessary and even harmful in some situations [START_REF] Chun | Interpolated boundary condition for lattice Boltzmann simulations of flows in narrow gaps[END_REF][START_REF] Yin | Mass and momentum transfer across solid-fluid boundaries in the lattice-Boltzmann method[END_REF] , which is consistent with the observation that mass leakage is more significant at high Reynolds numbers and/or at moving boundaries [START_REF] Sanjeevi | Choice of no-slip curved boundary condition for lattice Boltzmann simulations of high-reynolds-number flows[END_REF] . More recently, Ginzburg [START_REF] Ginzburg | Steady-state two-relaxation-time lattice Boltzmann formulation for transport and flow, closed with the compact multi-reflection boundary and interface-conjugate schemes[END_REF] proposed to include mass sources in the multi-reflection boundary scheme within a two-relaxation-time LB framework, and it was successfully applied to achieve mass balance across interfaces in advection-diffusion and Stokes flow simulations. Besides, within the same LB framework, a uniform normal mass flux was demonstrated to outperform a uniform mass-source in advection-diffusion simulations [START_REF] Ginzburg | Mass-balance and locality versus accuracy with the new boundary and interface-conjugate approaches in advection-diffusion lattice Boltzmann method[END_REF] . Some examples of mass leakage observation are summarised in Table I. Apparently, the observed mass leakage is supposed to depends on many factors, such as the boundary orientation, the Reynolds number value and specific flow characteristics. However, mechanisms responsible for the mass leakage at boundaries within the LB framework are still not clearly and rigorously characterised, which is the motivation of this paper.

This paper is focused on macroscopic mechanisms of mass leakage in LB simulations. Firstly, some basic characteristics of the LB method are presented in §II. Secondly, the mass leakage is defined as the net loss of distribution functions from a mesoscopic view in §III. Thirdly, local mass leakage at individual boundary nodes is mathematically expressed using macroscopic variables in §IV. Based on the local mass leakage analysis, averaged mass leakage over smooth boundaries is theoretically investigated in §V. After that, a mass correction scheme based on the averaged mass leakage is proposed in §VI. To validate the proposed theoretical analysis and assess the mass correction scheme, numerical experiments dealing with two-dimensional flows through straight and curved channels with moving boundaries are simulated using an in-house LB solver in §VII.

Finally, some conclusions are drawn in §VIII.

II. FUNDAMENTALS OF THE LB METHOD

As a fundamental basis to study mass leakage in LB simulations, some key features of the LB method are reminded with an emphasis on the mass conservation issue in this section. The LB method models fluid dynamics by considering a finite set of discrete velocities e i and distribution functions f i defined at a mesoscopic level. Specifically, the distribution functions are assumed to obey a set of coupled advection-relaxation partial differential equations, which can be interpreted as a discretisation in space, time and velocity space of the original continuous Boltzmann equation in kinetic theory. The relaxation term is a model of the exact collision term that appears in the Boltzmann equation, stemming from the idea that the net statistical effect of collisions is to drive the flow toward thermodynamic equilibrium.

Generally, the LB method can be formulated as follows:

f i (x,t) -f i (x -e i ∆t,t -∆t) = i (f -f eq ) ( 1
)
where ∆t is the time step, f and f eq are vectors containing f i (x -e i ∆t,t -∆t) and the equilibrium distribution functions f eq i , respectively, and i is a scalar collision function of ff eq for f i . Each component of e i is equal to either 0 or ±c with c = ∆x/∆t and ∆x being the grid spacing.

The left-hand side of Eq. ( 1) originates in a Lagrangian implementation of the linear convective streaming process, while the right-hand side is generally a nonlinear collision process (might also be linear when the physics is linear, e.g. advection-diffusion Stokes flow). By selecting adequate collision operator i , Eq. ( 1) is able to recover the classic single-relaxation-time (SRT) [START_REF] Bhatnagar | A model for collision processes in gases. I. Small amplitude processes in charged and neutral one-component systems[END_REF] , the multiple-relaxation-times (MRT) [START_REF] Lallemand | Theory of the lattice Boltzmann method: Dispersion, dissipation, isotropy, Galilean invariance, and stability[END_REF] , and the regularised (RLB) [START_REF] Latt | Lattice Boltzmann method with regularized pre-collision distribution functions[END_REF] LB methods. The most widely applied definition of f eq i is a truncated polynomial expansion of the original Maxwell-Boltzmann equilibrium distribution function, i.e.

f eq i = ω i ρ 1 + e i • u c 2 s + (e i • u) 2 2c 4 s - u • u 2c 2 s + O(Ma 3 ) . ( 2
)
where ω i is the i th weighting factor, c s is the lattice sound speed, ρ = ∑ f i , u = ∑ f i e i /ρ, and , Ma is the local Mach number u /c s . It should be noticed that c s /c is a tunable parameter depending on the lattice configuration as well as the adopted discrete velocities (e.g., c s /c = 1 √ 3 for regular lattices using the D2Q9 model) [START_REF] Ginzburg | Multireflection boundary conditions for lattice Boltzmann models[END_REF][START_REF] Qian | Lattice BGK models for Navier-Stokes equation[END_REF] . By using terms up to the first-order in f eq i (i.e. without non-linear terms in Eq. ( 2)) the LB method could restore to linear advection-diffusion equations.

Considering that mass leakage in this kind of linear LB processes has been well analysed by Ginzburg 7 , and demonstrated to be fully avoidable by choosing proper high order accurate boundary treatments in simple Stokes flow simulations [START_REF] Ginzbourg | Local second-order boundary methods for lattice Boltzmann models[END_REF][START_REF] Ginzburg | Steady-state two-relaxation-time lattice Boltzmann formulation for transport and flow, closed with the compact multi-reflection boundary and interface-conjugate schemes[END_REF] , hereafter this paper will be focused on more general and complex situation. Usually, by using terms up to second-order in f eq i , the LB method is sufficient to recover the low-Mach weakly compressible Navier-Stokes equations [START_REF] Lallemand | Theory of the lattice Boltzmann method: Dispersion, dissipation, isotropy, Galilean invariance, and stability[END_REF] . The higher order terms, e.g. O(Ma 3 ) in Eq. ( 2), can be used to account for thermal and compressible effect [START_REF] Feng | Hybrid recursive regularized thermal lattice Boltzmann model for high subsonic compressible flows[END_REF] .

Particularly, the LB equation ( 1) can be directly restored to the macroscopic mass equation by taking its zeroth moment, i.e.,

∂ ρ/∂t + ∇ • (ρu) = 0. (3)
However, it should be pointed out that this restoration is based on a complete stream process, i.e., all the neighbouring nodes are within the fluid domain, thus not taking boundary conditions into consideration. In addition, as discussed in the introduction, the reported boundary treatments within LB framework are not designed to directly satisfy the mass equation but are more likely focused on the boundary velocity condition. Consequently, mass leakage usually occurs along boundaries and cannot be directly described by the macroscopic mass equation.

To unveil the macroscopic mechanisms of mass leakage, it is important to express the nonequilibrium parts of the distribution functions, f ne i = f if eq i , as functions of their macroscopically determined equilibrium counterparts f eq i . By applying Taylor's series expansion, the LB equation

(1) can be rewritten as:

Df Dt + ∆t 2
D 2 f Dt 2 + O(∆t 2 ) = 1 ∆t L (f ne ) (4)
where f ne = ff eq , D/Dt is a diagonal matrix spanned by material derivative operators, e.g. ∂ /∂t + e i • along e i , and L is a matrix that contains all the linear collision operators i .

Applying D/Dt to both sides of Eq. ( 4), one obtains:

D 2 f Dt 2 = 1 ∆t L Df ne Dt - ∆t 2
D 3 f Dt 3 + O(∆t 2 ). (5)
Substituting f = f eq + f ne and Eq. ( 5) into Eq. ( 4), the non-equilibrium part can be expressed as follows

f ne = ∆tL -Df eq Dt + ∆t L -1 + 1 2 Id Df ne Dt + O(∆t 2 ) ( 6
)
where Id is the identity matrix. Applying D/Dt to both sides of this equation, it comes:

Df ne Dt = ∆tL -D 2 f eq D 2 t + ∆t L -1 + 1 2 Id D 2 f ne D 2 t + O(∆t 2 ) (7)
Now substituting Eq. ( 7) into Eq. ( 6), one obtains:

f ne = ∆tL -1 Df eq Dt + O(∆t 2 ) (8)
This relation shows that f ne contains time-dependent informations. Based on Eq. ( 8), f can be expressed as a function of f eq :

f = f eq + f ne = f eq + ∆tL -1 Df eq Dt + O(∆t 2 ) (9)
where f eq is explicitly determined by ρ and u through Eq. ( 2). According to Eq. ( 8), amplitude of the non-equilibrium distribution functions f ne scales as O(∆t Df eq Dt ), which is consistent with the widely used Chapman-Enskog analysis [START_REF] Lallemand | Theory of the lattice Boltzmann method: Dispersion, dissipation, isotropy, Galilean invariance, and stability[END_REF] . Noticeably, external body force is not included from Eq.

(1) to (9) because its effect on mass leakage can be directly expressed by the change of f eq due to the induced change of ρu. This is reasonable considering that the LB simulations have been verified to be insensitive to the adopted external force models [START_REF] Mohamad | A critical evaluation of force term in lattice Boltzmann method, natural convection problem[END_REF] .

At boundary nodes, reconstruction of f using Eq. ( 9) represents a second-order accurate boundary treatment. The mass leakage analysis presented hereafter is related to second-order accurate boundary treatments based on this equation in 2D.

III. DEFINITIONS OF MASS LEAKAGE

The mass leakage at boundaries in LB simulations is defined in this subsection. For the sake of convenience, a discrete velocity together with the associated distribution function is referred to as a LB link hereafter, and a fluid lattice node with wall-cut LB links is referred to as a boundary node (e.g. node B 1 shown in Fig. 1(a)).

Following the widely applied conception [START_REF] Onishi | A lattice Boltzmann model for polymeric liquids[END_REF][START_REF] Osmanlic | Lattice Boltzmann method for Oldroyd-B fluids[END_REF][START_REF] Küng | Comparison of passive scalar transport models coupled with the lattice Boltzmann method[END_REF][START_REF] Krivovichev | Analysis of the parametric models of passive scalar transport used in the lattice Boltzmann method[END_REF][START_REF] Krivovichev | Parametric schemes for the simulation of advection process in finite-differencebased single-relaxation-time lattice Boltzmann methods[END_REF][START_REF] Zhao | Toward fully conservative hybrid lattice Boltzmann methods for compressible flows[END_REF] , the mass leakage is defined as the difference between the sum of the distribution functions streaming out from the fluid domain ( e.g. f 1 and f 3 in Fig. 1(a)) and that of those streaming into the fluid domain (e.g. f 2 and f 4 Fig. 1(a)). In addition, the normal mass flux due to the movement of the boundary should be deducted. Consequently, the local mass leakage per time unit per area unit can be defined as

E(x) = ∆x D ∆S∆t ∑ x+e i ∆t∈S [ f i (x) -f ī (x + e i ∆t)] + ρ w U • n ( 10
)
where x is the coordinate vector of a boundary node, D is the number of spatial dimensions (D = 2 here), S is the solid side of the boundary, f i and f ī corresponds to discrete velocities in opposite directions, i.e. e i and -e i , respectively, ρ w and U are the fluid density and velocity at the projection of x onto the boundary, n is the unit normal vector (pointing to the fluid side) at the projection point, ∆x is the local grid spacing, and ∆S is the projection area of a lattice on the boundary. Clearly, the mass leakage E defined in Eq. ( 10) is a density flux. Actually, it can be considered as an approximation to the conventional density flux loss (-ρu + ρ w U ) • n:

E(x) = ∆x D ∆S∆t ∑ x+e i ∆t∈S [ f i (x) -f ī (x + e i ∆t)] + ρ w U • n ≈ - ∆x D ∆S∆t ∑ i f i Sign(e i • n) + ρ w U • n ≈ -∑ i f i e i • n + ρ w U • n = (-ρu + ρ w U ) • n (11)
It should be noticed that Eq. ( 11) is an effective but not the unique measurement of the potential mass leakage, for example, at a steady flow state the density at t -∆t could be directly used to determine the mass leakage as ρ(t -∆t) -∑ i f i (t). By implicitly included this measurement as a correction source term within a two-relaxation-time LB framework, Ginzburg [START_REF] Ginzburg | Consistent lattice Boltzmann schemes for the Brinkman model of porous flow and infinite Chapman-Enskog expansion[END_REF] proposed a LB scheme successfully reproducing the Poiseuille flow without mass leakage. Apparently, the steady state intended measurement is not suitable for general LB simulations as concerned in this paper. Based on the local mass leakage at individual boundary nodes defined in Eq. ( 10), the averaged mass leakage over a smooth boundary can be now expressed as:

Ē = ∑ E∆S ∑ ∆S (12)
where the summation is conducted over all the involved boundary nodes.

IV. LOCAL MASS LEAKAGE AT INDIVIDUAL BOUNDARY NODES

Based on the fundamentals of the LB method presented in §II, the local mass leakage at individual boundary nodes defined by Eq. ( 10) is now theoretically expressed in terms of macroscopic variables in this section. Firstly, mass leakage at boundary nodes with LB links regularly cut by aligned planar boundaries (referred to as "regular boundary nodes", see Fig. 1(b)) is analysed in §IV A. Based on that, mass leakage at boundary nodes with LB links irregularly cut by general smooth boundaries (referred to as "irregular boundary nodes", see Fig. 2(a)-(c)) is considered in §IV A.

A. Regular boundary nodes

Due to the symmetry of wall-cut LB links at a regular boundary node, the zero and first order moments of the links can be directly expressed by macroscopic variables. Hence, regular boundary nodes are an ideal starting point to express the mass leakage by macroscopic variables.

By applying the symmetry of the wall-cut LB links at regular boundary nodes, the mass leakage can be expressed as (see Appendix A for details):

E(x) ≈ -ρu m + ∆x 6 c∂ m ρ + ∆x 3 ∂ m (ρu m ) + ∆x 6 ∇ • (ρu) + ∆x 4c [∂ m (ρu 2 m ) + ∇ • (ρu m u)] + ρ w U • n (13)
where the sub-index m indicates components along the symmetry axis of the wall-cut LB links (referred to as "main direction" with a unit vector n m pointing to the fluid domain). It is worth noticing that ρu m is the major part of the momentum recovered by the wall-cut LB links, and expression ( 13) is accurate without requiring n m = n and will be used for irregular boundary nodes in the next subsection.

Considering that n m = n is satisfied for aligned boundaries, Eq. ( 13) can be further simplified as:

E(x) = ∆x 3 -∆w ∂ n (ρu n ) + ∆x 6 c∂ n ρ + ∆x 6 ∇ • (ρu) + ∆x 4c [∂ n (ρu 2 n ) + ∇ • (ρu n u)] (14)
where ∆w is the wall distance (∆w < ∆x), and ∂ n and u n are the normal gradient and velocity components, respectively.

The last relation (14) clearly shows that the mass leakage along aligned boundaries mainly arises from the normal gradient of density ρ, normal momentum ρu n , and the corresponding kinetic energy ρu 2 n , as well as the divergence of the momentum ∇ • (ρu) (compressibility effect) and the normal momentum flux ∇ • (ρu n u). Moreover, amplitude of the mass leakage is in the order of O(∆x). As a consequence, in the case of weakly compressible flows with fixed non-slip boundaries, all the terms in Eq. ( 14) are all supposed to be negligibly small, as is the resultant local mass leakage. It should be noticed that this conclusion does not conflict with the reported mass leakage along aligned boundaries [START_REF] Oulaid | On the origin of numerical errors in the bounce-back boundary treatment of the lattice Boltzmann method: A remedy for artificial boundary slip and mass leakage[END_REF][START_REF] Yu | Modified curved boundary scheme for two-phase lattice Boltzmann simulations[END_REF] because they adopted an external force at boundary nodes, thus yielding a significant ∂ n (ρu n ) in Eq. ( 14).

B. Irregular boundary nodes

The mass leakage at irregular boundary nodes is now addressed. The analysis is conducted by approximating irregular boundary nodes as regular ones, and then estimating the mass leakage by applying the analysis proposed for regular nodes in §IV A. 13), the mass leakage can be estimated as (see Appendix B):

Solid

E(x) = ± ρu t tan θ + O(∆xc ∇ρ ) + O (∆x ∇(ρu) ) + O ∆x ∇(ρuu) c s ( 15
)
where u t is the tangential velocity component.

For irregular boundary nodes associated with irregular cell pairs, they can be approximated as regular boundary nodes through two kinds of operations. Firstly, the links between neighbouring boundary nodes could be treated as virtual cut links without causing extra net mass leakage. For example, as shown in Fig. 2(a), by considering the virtually cut links e 5 and e 6 between B 2 and B 3 , the node B 3 can be treated as an approximated regular boundary node with a horizontal main direction. Secondly, wall-cut LB links at one boundary node can be shifted to its neighbouring boundary nodes to complete the remained approximation to regular boundary nodes. For example, as shown in Fig. 2(a), by shifting the cut links e 3 and e 4 at B 1 to B 2 , the total cut links, including the virtual e 5 and e 6 , exactly approximate B 2 as two regular boundary nodes with different main directions (one is horizontal and the other is vertical). Mass leakage at the approximated regular node pairs, e.g. those at B 2 as analysed, can be estimated by Eq. ( 13) as (see Appendix B):

E(x) =ρu t tan(±θ ∓ 45 • ) + O(∆xc ∇ρ ) + O (∆x ∇(ρu) ) + O ∆x ∇(ρuu) c s (16)
Clearly, relations (15) and (16) shows that the mass leakage at both kinds of irregular boundary nodes are dominated by terms proportional to ρu t . Recalling Eq. ( 13), it can be concluded that the dominating terms in Eqs. ( 15) and ( 16 To this point it should be pointed out that, from Eq. ( 13) to ( 16), the potential significant mass leakage ∝ ρu t is inevitable for general non-aligned boundaries because it originates from the intrinsic deficiency of the incomplete distribution functions at boundary nodes in restoring the momentum vector properly. In addition, the dependence of mass leakage on different boundary treatments is supposed to be included in the potentially different coefficients before the qualitative terms (e.g. O(∆xc ∇ρ )), and this can be expected to holds for the corresponding averaged mass leakage discussed in the followed §V.

V. AVERAGED MASS LEAKAGE OVER SMOOTH BOUNDARIES A. Estimates of the averaged mass leakage

Based on the estimates for local mass leakage at individual boundary nodes derived in §IV, the averaged mass leakage over smooth boundaries is analysed in this section.

Firstly, for aligned boundaries, regular boundary nodes are uniformly distributed, and the averaged mass leakage is exactly the same as the local value at an individual node given by Eq. ( 14).

Secondly, for 45 • inclined boundaries, all the boundary nodes are associated with irregular cells, and the mass leakage can be estimated by Eq. ( 16) taking θ = 45

• , i.e. Ē = O(∆xc ∇ρ ) + O (∆x ∇(ρu n ) ) + O ∆x ∇(ρuu) c s (17)
Finally, for general smooth boundaries, the boundary nodes can be decomposed into a series of serrated cells (see Fig. 3(a)-(d)), and the overall averaged mass leakage can be estimated from that of the serrated cells. Specifically, a general serrated cell can be deemed as a linear supperposition of four basic types, and so is the averaged mass leakage. By combining the local mass leakage (Eqs. ( 15) and ( 16)), the definition of averaged mass leakage (Eq. ( 12)) and the geometry relationships shown in Fig. 3

       E sc =O(∆xc ∇ρ ) + O (∆x ∇(ρu) ) + O ∆x ∇(ρuu) c s E ∆s =∆s∂ n (ρu n ), E ∆θ = ρu t tan ∆θ = O(ρu t ∆x L ), E R c = O ρ u ∆x R c (18)
It should be noted that E ∆s , E ∆θ and E R c are not just for a local serrated cell, but have taken the whole smooth boundary into consideration (see Appendix C). Consequently, the resultant averaged mass leakage for smooth boundaries can be estimated as the sum of Eq. ( 18):

Ē = O(ρu t ∆x L ) + O ρ u ∆x R c + O(∆xc ∇ρ ) + O (∆x ∇(ρu) ) + O ∆x ∇(ρuu) c s (19)
As shown by this relation, the averaged mass leakage mainly originates in several sources, including the tangential momentum ρu t , the boundary curvature and the momentum ρu/R c and the gradients of density ρ, momentum ρu and momentum flux tensor ρuu. Moreover, all the terms are proportional to the grid spacing O(∆x), indicating that the averaged mass leakage can be reduced by refining grids. In addition, it should be noticed that all the terms appearing in Eq. ( 19) are significantly dependent on the boundary orientation. For example, the first term caused by ∆θ is zero when the boundary is inclined at

θ = arctan 1 N with N = 1, 2, 3 • • • , e.g. θ = 45 • corresponding to N = 1.
Physically, moving boundaries induce significant velocities at boundary nodes, and high

Reynolds numbers are associated with strong gradients near the boundaries. According to Eq. ( 19), these two situations are expected to induce significant averaged mass leakage, which is consistent with the previous observations reported in [START_REF] Sanjeevi | Choice of no-slip curved boundary condition for lattice Boltzmann simulations of high-reynolds-number flows[END_REF] .

B. Examples on two academic cases: Poiseuille and Taylor-Couette flows

The above expressions for mass leakage can be used to analyse the behaviour of mass leakage in simple flow configurations, for which analytical solutions are available. Hereafter, averaged mass leakage along boundaries for steady laminar Poiseuille flow and Taylor-Couette flow are estimated thanks to Eq. ( 19).

For Poiseuille flow, considering ∇ρ = 0, u n = 0, ∂ t u t = 0 and ∇ • (ρu) ≈ 0, the averaged mass leakage along the straight walls can be expressed as:

Ē =O(ρu t ∆x L ) + O (∆xρ∂ n u t ) + O ∆x ρu t ∂ n u t c s (20)
For the steady laminar weakly compressible Taylor-Couette flow with a static outer cylinder (pressure gradient offseted by external body force), one has ∇ρ = 0, u n = 0 and ∇ • (ρu) ≈ 0, and the averaged mass leakage along the curved inner wall can be expressed as:

Ē = O(∆xρΩ) + O ∆x ρr 1 Ω 2 c s (21)
where the reference length L in Eq. ( 19) is represented by the inner radius r 1 , and Ω is the angular velocity of the inner cylinder.

VI. A FIRST-ORDER ACCURATE MASS CORRECTION SCHEMES

Corresponding to the local and averaged mass leakage estimates, there are two mass correction schemes, i.e. the local mass correction (LMC) scheme and the averaged mass correction (AMC) scheme.

With the LMC, the local mass leakage is used to satisfy mass conservation by applying a correction as:

∆ρ(x) = ∆t∆S ∆x D E(x) = ∑ x+e i ∆t∈S [ f i (x) -f ī (x + e i ∆t)] + ∆t∆S ∆x D ρ w U • n ( 22
)
where E is defined in Eq. (10). For static or tangentially moving boundaries (i.e. U • n), this scheme complies with the way of mass conservation of the well-known bounce-back method, has been used by Sanjeevi et al. [START_REF] Sanjeevi | Choice of no-slip curved boundary condition for lattice Boltzmann simulations of high-reynolds-number flows[END_REF] , and is similar to most of the reported mass correction schemes [START_REF] Bao | A mass conserving boundary condition for the lattice Boltzmann equation method[END_REF][START_REF] Coupanec | A mass conserving boundary condition for the lattice Boltzmann method for tangentially moving walls[END_REF][START_REF] Yu | Modified curved boundary scheme for two-phase lattice Boltzmann simulations[END_REF] .

According to the analysis in §IV, the local mass leakage could be proportional to the tangential velocity amplitude, and thus the LMC scheme adds a zero-order perturbation to the solution of density.

In contrast, the averaged mass leakage is promising to facilitate a first-order accurate correction scheme because it has an amplitude in the order of O(∆x). With the AMC, the mass correction can formulated as:

∆ρ(x) = ∆t∆S ∆x D Ē = ∑ Ω ∑ x+e i ∆t∈S [ f i (x) -f ī (x + e i ∆t)] ∑ Ω ∆S + ∆t∆S ∑ Ω ∆S ∑ Ω (ρ w U • n∆S) ( 23
)
where Ē is defined in Eq. ( 12), and Ω represents a concerned boundary. It should be noticed that ∆ρ(x) is implied locally, but the AMC only provide a global mass conservation, which is consistent with local mass flux correction proposed by Ginzburg [START_REF] Ginzburg | Steady-state two-relaxation-time lattice Boltzmann formulation for transport and flow, closed with the compact multi-reflection boundary and interface-conjugate schemes[END_REF] . The two schemes will be compared in the followed §VII.

VII. NUMERICAL VALIDATION ON THE CHANNEL AND TAYLOR-COUETTE FLOWS

A preliminary validation of the proposed scaling laws for the mass leakage error is achieved by performing 2D simulations of two academic flows with straight and curved boundaries, e.g.

the Poiseuille flow and the Taylor-Couette flow, respectively. Meanwhile, the influence of mass leakage on the numerical solutions as well as the effectiveness of the LMC and AMC schemes are investigated. Without loss of generality, the multiple-relaxation-times (MRT) LB method together with a second-order accurate immersed boundary method will be used for the validation.

An in-house C/C++ solver [START_REF] Xu | A novel geometry-adaptive Cartesian grid based immersed boundary-lattice Boltzmann method for fluid-structure interactions at moderate and high Reynolds numbers[END_REF] applying the LB method based on multi-block Cartesian grids is used to perform the simulations until otherwise specified. In this solver, the multiple-relaxationtimes (MRT) [START_REF] Humières | Multiple-relaxation-time lattice Boltzmann models in three dimensions[END_REF] method together with the D2Q9 lattice is used for two-dimensional simulations.

The no-slip boundary condition is implemented by an interpolation-based scheme similar to the one discussed in Ref. [START_REF] Filippova | Grid refinement for lattice-BGK models[END_REF][START_REF] Feng | Hybrid recursive regularized thermal lattice Boltzmann model for high subsonic compressible flows[END_REF] . In this approach, as shown in Fig. 4, the macroscopic variables (ρ and The averaged mass leakage for the Poiseuille flow is estimated by Eq. (20). Since the serrated cell at the lower and upper walls (denoted by indices l and u, respectively) are rotational symmet-rical by 180 • , the tangential velocity u t in Eq. ( 20) is in opposite signs to estimate their averaged mass leakage. Therefore, Ēl -Ēu is supposed to be contributed by the linear terms O(ρu t ∆x L ) and O(∆xρ∂ n u t ), and the net mass leakage Ēl + Ēu is supposed to be contributed by the non-linear term O(∆xρu t ∂ n u t /c s ) as well as the coupling between the linear terms and the asymmetrical geometry uncertainties (e.g. ∆s in Eq. ( 18) could be different for the two walls).

Effect of geometry orientation

To investigate the mass leakage dependence on geometry orientation, the Poiseuille flow with tangentially moving walls inclined from θ = 0 • to 45 • by an increment of 5 • is simulated. A finite boundary velocity U b = 0.5U re f is enforced to provide significant velocity at the boundary nodes.

To illustrate and validate the numerical solution, grid distribution at θ = 30 mass leakage is consistent with the analysis in Ref. [START_REF] Ginzburg | Consistent lattice Boltzmann schemes for the Brinkman model of porous flow and infinite Chapman-Enskog expansion[END_REF] that the non-linear terms in the equilibrium definition of distribution functions cause unavoidable mass leakage.

Effect of tangential velocity

To emphasise the effect of tangential velocity on the mass leakage, simulations with the walls translating at a finite tangential velocity U b varying from 0 to 0.5U re f are performed in this subsection. The normal gradient of tangential velocity ∂ n u t is set to be constant so that u t is the sole changing variable in Eq. ( 20). In addition, to clarify effect of the angle bias ∆θ , i.e. the term O(ρu t ∆x L ) in Eq. ( 20), configurations at θ = arctan 1 2 (≈ 26.6 • , ∆θ = 0) and θ = 30 • (∆θ = 0) are simulated and compared. 15) and (16). Meanwhile, the net mass leakage Ēl + Ēu in both figures exhibits similar linear dependence versus U b , which can be explained by the nonlinear term O(∆xρu t ∂ n u t ) in Eq. (20). The main difference between the two figures is that Ēl -Ēu remains almost unchanged about 3 × 10 -5 ρU re f for θ = arctan 1 2 (Fig. 8(a)), while that shows a significant linear dependence on U b for θ = 30 • (Fig. 8(a)). This clearly confirms the effect of ∆θ is definitely well estimated by the term O(ρu t ∆x L ).

Effect of velocity gradient

To measure the effect of velocity gradient on the mass leakage, simulations with the maximum relative incoming velocity u max -U b varying from 0 to U re f are now performed. Without loss of generality, the solid walls are set to move tangentially at U b /U re f = 0.5 so that ∂ n u t is the sole changing variable in Eq. ( 20). θ = arctan 1 2 is chosen to vanish the term O(ρu t ∆x) in Eq. ( 20), and thus to avoid the uncertainty induced by ∆θ . velocity (u-U b ) is well reproduced as expected. For the mass leakage shown in Fig. 9(b), the local amplitude |E| max at both walls exhibits linear dependence on ∂ n u t with a small amplitude variation (2%). This is consistent with that ∂ n u t only emerges linearly in minor terms in Eqs. ( 15) and ( 16).

Comparatively, the linear averaged mass leakage Ēl -Ēu is almost proportional to ∂ n u t , which could be attributed to the linear term O(∆xρ∂ n u t ) in Eq. (20). Differently, the net mass leakage Ēl + Ēu is well described by a parabolic curve, which can be explained by the considerable increase of u t due to ∂ n u t in the nonlinear term O(∆xρu t ∂ n u t ). By taking this increment into consideration, the nonlinear term is changed to O[∆xρ(u t + ∆w∂ n u t )∂ n u t ], including eventually a quadratic term for significant ∂ n u t (the relative increase of u t at (u max -U b )/U re f = 1 reaches a considerable amplitude of about 10%U b in Fig. 9(a)).

Effect of grid spacing ∆x

To examine effect of ∆x on the mass leakage, simulations using four different grid spacings near the solid walls, i.e. ∆x/h = 0.025, 0.0125, 0.00625 and 0.0003125, are performed. Without loss of generality, U b is taken equal to 0.5U re f , and the channel is inclined at θ = 30 • to include all the terms in Eq. (20).

Figure 10 displays the mass leakage as functions of ∆x. The local amplitude of the mass leakage |E| max at both walls is slightly affected (1%) in a linear way, which can be attributed to the ∆x∂ n u t -related terms in Eqs. ( 15) and (16). Meanwhile, the averaged mass leakage at both walls exhibits a general linear growth versus ∆x with some observable deviations. This can be deemed as reasonable validation of Ē ∝ ∆x (Eq. ( 20)) considering that refining the grid could generate substantial uncertainties to the boundary node distribution, e.g. ∆s and ∆θ in Eq. ( 18).

Effect of mass leakage on wall-modelled turbulent flow through straight channel

As demonstrated above ( §VII A 1-VII A 4), mass leakage in the laminar Poiseuille flow simulations is well predicted by the proposed theory, but the amplitude is too low (0.01%) to significantly impact the solution. Thus, to emphasise the effect of mass leakage, a wall-modelled RANS simulation of turbulent flow at Re = 5 × 10 4 (Re tau ≈ 2000) is performed. At steady state, this solution is similar with the Poiseuille flow in terms of mass leakage estimate, i.e. still described by Eqs. ( 15), ( 16) and (20). Moreover, the velocity and its gradients at the boundary nodes are supposed to be large enough to involve a significant mass leakage. Finally, the use of wall function increases the sensitivity of the solution to the boundary treatment. Without loss of generality, the channel is inclined at θ = 30 • with L/h = 200, and a uniform grid spacing ∆x/h = 0.05 is used. Instead of a parabolic velocity profile at the inlet, a uniform inlet velocity is imposed.

The simulation is conducted using the ProLB software [START_REF] Feng | Hybrid recursive regularized thermal lattice Boltzmann model for high subsonic compressible flows[END_REF] in 2D by extending one layer of grid along the third axis with periodic setup. In ProLB, the D3Q19 LB model is solved for the fluid dynamics by applying the hybrid recursively regularised strategy [START_REF] Feng | Hybrid recursive regularized thermal lattice Boltzmann model for high subsonic compressible flows[END_REF] As can be observed, without mass correction, the local mass leakage varies significantly with an amplitude about two orders higher than the averaged values. This observation can be well explained by the proposed theory ( Eqs. ( 15), ( 16) and ( 19) ). Meanwhile, the sectional mass flux shown Fig. 11(b) exhibits a loss up to 1.3%, and the AMC scheme apparently outperforms the LMC scheme in controlling the mass flux loss. Fig. 12(a)-(b) display the friction coefficient and sectional velocity profile, respectively. As can be observed, without mass correction, the mass leakage causes significant error to the solution, e.g. the flow symmetry is ruined and both the friction and streamwise velocity significantly deviate from the reference data. Notably, these problems are well cured by the AMC scheme, while the LMC scheme leads to worse result.

B. Taylor-Couette flow

In order to validate the proposed theoretical analysis of mass leakage on curved boundaries, the steady laminar Taylor-Couette flow is considered. The boundary set up is shown in Fig. 13(a). The inner and outer cylinders' radii are denoted as r 1 and r 2 , respectively. The inner cylinder rotates clockwise at an angular velocity Ω. The half channel height, i.e. h = (r 2r 1 )/2, is used as the reference length. A reference value of angular velocity Ω re f is defined as Re = Ω re f h/ν = 50. Accordingly, a reference velocity is defined as U re f = Ω re f h. A volumic forcing term equal to ρu 2 t /r (with r being the local radial distance) is used to ensure negligible density gradients. The minimum grid spacing is fixed at ∆x/h = 0.02 unless otherwise specified. The generated multiblock grid is shown in Fig. 13(b).

The averaged mass leakage is assumed to be described by Eq. ( 21). Due to the symmetry of the geometry setup, the serrated cells are symmetrically distributed around every line inclined at 0 • or 45 • . Therefore, the linear averaged mass leakage (denoted as ĒL , corresponding to O(∆xρΩ) in Eq. ( 21)) is evaluated as the difference between two adjacent 45 • arcs. In contrast, the net mass leakage of a 90 • arc is denoted as ĒNL , corresponding to O(∆xρr 1 Ω 2 ) in ( 21).

Effect of boundary curvature

To validate effect of boundary curvature on the mass leakage, simulations with r 1 ranging from h to 6h are performed in this subsection. The angular velocity is set as Ω = Ω re f . 15) and (16). Meanwhile, the linear averaged mass leakage ĒL exhibits significant oscillations without an obvious increasing or decreasing trend. The oscillations can be attributed to the substantial uncertainties associated with the boundary node distribution caused by the change of geometry configuration. The unclear trend is consistent with the term O(ρ∆xΩ) (see Eq. ( 21))

which is independent of r 1 . In contrast, the net averaged mass leakage ĒNL shown in Fig. 14(c) is almost proportional to r 1 /h, which can be well explained by the term O(∆xr 1 Ω 2 /c s ) in Eq. ( 21).

Effect of angular velocity

To quantify effect of the angular velocity on the mass leakage, simulations with Ω ranging from 0.125Ω re f to Ω re f are performed with r 1 = 2h. 15) and ( 16), and the later is consistent with the term O (∆xρΩ) in Eq. (21). In contrast, the net averaged mass leakage ĒNL shown in Fig. 15(b) is almost proportional to Ω 2 , which is consistent with the term O ∆xρr 1 Ω 2 in Eq. (21).

Effect of grid spacing ∆x

In order to measure the effect of grid spacing ∆x on the mass leakage, simulations are performed considering five different values of ∆x at the wall, i.e. ∆x/h = 0.04, 0.02, 0.01, 0.005 and 0.0025. Considering the substantial uncertainties associated with boundary node distribution caused by the grid refinement, behaviours of both ĒNL and ĒL reasonably validate the theoretical prediction Ē ∝ ∆x (Eq. ( 21)).

Effect of mass leakage on the solution

To emphasise the effect of mass leakage on the simulation of Taylor-Couette flow, a case without body force using a uniform coarse grid spacing ∆x/h = 0.1 is conducted. The body force is cancelled to clarify the influence on the pressure estimate. Other parameters are set as: r 1 = h and Ω = Ω re f . Fig. 17 with an amplitude comparable to ρU re f , but the negative averaged value has an amplitude about three order lower than ρU re f . These observations can be well explained by the proposed theory ( Eqs. ( 15), ( 16) and ( 21) ). Consistent with the negative average mass leakage in Fig. 17(a), without mass correction, the total mass M within the circular channel increases linearly with time, to an increment of 17% at U re f t/h = 600. In contrast, both the local and averaged schemes are effective to remove the total mass leakage. Fig. 18(a) displays the pressure coefficient profile at x/h = 0. As can be observed, without mass correction, a significant shift of pressure is produced due to the significant mass increment shown in Fig. 17 scheme in terms of accuracy and smoothness. The superiority of the averaged scheme is further corroborated by the pressure coefficient contours shown in Fig. 18(b) and(c).

VIII. CONCLUSIONS

In this paper, mass leakage within the general lattice Boltzmann framework is theoretically analysed for boundary treatments of at least second-order accuracy. The local mass leakage, defined from a mesoscopic viewpoint, is theoretically approximated using macroscopic variables at regular (aligned boundaries) and irregular boundary nodes. Based on that, the averaged mass leak-age over smooth boundaries is mathematically estimated from that of a basic serrated structure of boundary nodes.

The present theoretical results show that the mass leakage in LB simulations originates from the deficiency of the wall-cut links in recovering momenta at the boundary nodes. Unexpectedly, from the macroscopic point of view, the local mass leakage at individual boundary nodes is dominated by terms proportional to the tangential momentum (i.e. ρu t ). In contrast, the averaged mass leakage involves several different terms related to the tangential momentum ρu t , product of the boundary curvature and the momentum ρu/R c and the gradients of density ρ, momentum ρu and momentum flux tensor ρuu. Moreover, the present analysis shows its amplitude is proportional to the local grid spacing. Based on that, a mass correction scheme is proposed. Besides, both the local and averaged mass leakage significantly depends on the boundary orientation. Particularly, the mass leakage at fixed no-slip aligned boundaries is negligibly small. The good agreement between the proposed theoretical estimates and numerical results extracted from the simulations of the two-dimensional weakly compressible flows through moving straight and curved channels provides a preliminary validation of the theory. In addition, the numerical results show that LB solutions can be significantly impacted by mass leakage. However, the averaged mass correction scheme designed from the present theoretical analysis performs well in curing the mass leakage problems in the considered cases.

The theoretical analysis proposed in this paper unveils the basic mechanisms of mass leakage in the LB framework. However, its validity is limited to 2D flow academic configurations and its validation in more complex situations as well as its extension to 3D still requires future research efforts as its validation. In addition, the complex dependance of the evitable local mass leakage on local flow characteristics indicates that it might be necessary to includes specific characteristic like "steady flow" into the configuration of boundary treatments as preliminarily demonstrated by Ginzburg & d'Humières [START_REF] Ginzbourg | Local second-order boundary methods for lattice Boltzmann models[END_REF][START_REF] Ginzburg | Consistent lattice Boltzmann schemes for the Brinkman model of porous flow and infinite Chapman-Enskog expansion[END_REF][START_REF] Ginzburg | Study of simple hydrodynamic solutions with the Two-Relaxation-Times lattice Boltzmann scheme[END_REF] .

Appendix A: Local mass leakage at regular boundary nodes

Theoretically, mass leakage defined by Eq. ( 10) can be decomposed into three parts by splitting f into then equilibrium f eq and non-equilibrium f ne parts, i.e.

E(x) =E eq (x) + E ne (x) + ρ w U • n (A1) with            E eq (x) = ∆x D ∆S∆t ∑ x+e i ∆t∈S f eq i (x) -f eq ī (x + e i ∆t) E ne (x) = ∆x D ∆S∆t ∑ x+e i ∆t∈S f ne i (x) -f ne ī (x + e i ∆t) (A2)
It should be noted that f ne i could contribute to mass leakage during the stream process (involving two neighbouring nodes) despite that its zero and first order moments are theoretically zero at a given node.

Applying the definition of f eq i given by Eq. ( 2), the equilibrium error part E eq can be expanded

as:

E eq (x) = ∆x D ∆S∆t ∑ x+e i ∆t∈S 2ω i ρe i • u c 2 s + 2∆tω i e i ∇ • ρe i • u c 2 s -∆te i • f eq i + O(∆t 2 ) = ∆x D ∆S∆t ∑ x+e i ∆t∈S 2ω i ρe i • u c 2 s + ∆tω i e i • ρe i • u c 2 s -∆tω i e i • ρ + O(∆t 2 ) - ∆x D ∆S ∑ x+e i ∆t∈S ω i e i • ∇(ρe i • u) 2 2c 4 s -ω i ρ e i • ∇(ρu • u) 2c 2 s + O ρMa 2 ∇u (A3)
Clearly, Eq. (A3) shows that high-order terms O(ρMa n ) in f eq i could induce a mass leakage that scales as O Ma n-1 ρ ∇u . Without loss of generality, in the rest of this paper, Ma-related high-order terms, e.g. O(Ma 3 ), will be neglected because their effect can be inferred from those of O(Ma) and O(Ma 2 ). Hence, by omitting the high-order terms O Ma 2 ρ ∇u and O(∆t 2 ), Eq.

(A3) simplifies as:

E eq (x) = ∆x D ∆S∆t ∑ x+e i ∆t∈S 2ω i ρe i • u c 2 s -∆tω i e i • ρ + ∆tω i e i • ρe i • u c 2 s - ∆x D ∆S ∑ x+e i ∆t∈S ω i e i • ∇[ρ(e i • u) 2 ] 2c 4 s -ω i e i • ∇(ρu • u) 2c 2 s (A4)
As illustrated in Fig. 1(b), at a boundary node along an aligned boundary, the wall-cut LB links is symmetrical around an axis, e.g. the mth axis, with a unit vector n m pointing to the fluid side (referred to as "main direction"). Accordingly, all wall-cut links share a non-zero mth velocity component e i,k = ±c, and the other components e i, j =m are symmetrically distributed among -c, 0 and c. Considering the symmetry of wall-cut links and ∆S = ∆x D-1 , the linear terms of u on the right-hand side of Eq. (A4) can be evaluated as:

                                         ∆x D ∆S∆t ∑ x+e i ∆t∈S 2ω i ρe i • u c 2 s =c Q ∑ i=1 ω i -cρu m c 2 s = -ρu m ∆x D ∆S∆t ∑ x+e i ∆t∈S (∆tω i e i • ρ) = ∆x 2 Q ∑ i=1 (-cω i ∂ m ρ) = - ∆x 6 c∂ m ρ ∆x D ∆S∆t ∑ x+e i ∆t∈S ∆tω i e i • ρe i • u c 2 s = ∆x 2 Q ∑ i=1 ω i e 2 i,m ∂ m (ρu m ) c 2 s + ∆x ∑ x+e i ∆t∈S, j =m ω i e 2 i, j ∂ j (ρu j ) c 2 s = ∆x 2 ∂ m (ρu m ) + ∆x 6 ∂ j (ρu j ) = ∆x 3 ∂ m (ρu m ) + ∆x 6 ∇ • (ρu) (A5)
where ∂ m = n m • ∇, u m = u • n m , j denotes the axis perpendicular to the main direction (or the mth axis), Q is the number of the discrete velocities (e.g. Q = 9 for the D2Q9 model), and the following relations hold:

                   e 2 i,m = c 2 , c 2 s = c 2 /3, Q ∑ i=1 ω i e 2 i, j c 2 s = 1 e i,m • n m = -c for x + e i ∆t ∈ S 2 ∑ x+e i ∆t∈S, j =m ω i e 2 i, j c 2 s Q ∑ i=1 ω i e 2 i, j c 2 s = c 2 s c 2 (A6)
In the last relation of Eq. (A6), compared to the denominator, the nominator does not include the discrete velocities perpendicular to the main direction, e.g. the links between B 1 and B 2 in Fig.

1(b).

The relations in Eq. (A6) are valid for the D2Q5 and D2Q9 models as well as D3Q19 and D3Q27 models when the solution is invariant along the third axis. Equation (A5) demonstrates that the wall-cut LB links at a regular boundary node are only able to restore velocity and gradient components along the main direction, while the other components are completely omitted. Still, by applying Eq. (A6), the non-linear terms in Eq. (A4) can simplified as:

∆x D ∆S ∑ x+e i ∆t∈S ω i e i • ∇[ρ(e i • u) 2 ] 2c 4 s -ω i e i • ∇(ρu • u) 2c 2 s = ∆x 2 Q ∑ i=1 ω i -ce 2 i,m ∂ m (ρu 2 m ) 2c 4 s + ∆x ∑ x+e i ∆t∈S, j =m ω i e i,m e 2 i, j 2c 4 s [∂ m (ρu 2 j ) + ∂ j (ρu m u j )] - ∆x 2 Q ∑ i=1 ω i -c∂ m (ρu 2 m + ρu 2 j ) 2c 2 s = - ∆x 2 c ∂ m (ρu 2 m ) 2c 2 s + c 2 s c 2 ∂ m (ρu 2 j ) + ∂ j (ρu m u j ) 2c 2 s - ∂ m (ρu 2 m + ρu 2 j ) 2c 2 = - ∆x 2 c ∂ m (ρu 2 m ) 3c 2 s + ∂ j (ρu m u j ) 6c 2 s = - ∆x 4c [∂ m (ρu 2 m ) + ∇ • (ρu m u)] (A7)
Consequently, by substituting Eqs. (A5) and (A7) into Eq. (A3), the mass leakage directly induced by f eq can be expressed as:

E eq (x) = -ρu m + ∆x 6 c∂ m ρ + ∆x 3 ∂ m (ρu m ) + ∆x 6 ∇ • (ρu) + ∆x 4c [∂ m (ρu 2 m ) + ∇ • (ρu m u)] (A8)
The mass leakage due to the non-equilibrium distribution function, i.e. E ne in Eq. (A2), is now analysed. Still considering regular boundary nodes along aligned boundaries (∆S = ∆x D-1 ), by applying the symmetry of discrete velocities, it can be rewritten as:

E ne (x) = - Q ∑ i=1 [e i,m f ne i (x)] - ∆x ∆t ∑ x+e i ∆t∈S [ f ne i (x + e i ∆t) -f ne i (x)] = -∆x ∑ x+e i ∆t∈S D f ne i (x) Dt + ∆t D 2 f ne i (x) Dt 2 + O(∆t 2 ) (A9)
Now using Eq. ( 8), one obtains :

E ne (x) = -c∆x 2 ∑ x+e i ∆t∈S L -1 i D 2 f eq (x) Dt 2 + O(∆x 3 ) (A10)
Substituting Eqs. (A8) and (A10) into Eq. (A1) and omitting high order terms, the resultant mass leakage is dominated by E eq , and can therefore be expressed as:

E(x) ≈ -ρu m + ∆x 6 c∂ m ρ + ∆x 3 ∂ m (ρu m ) + ∆x 6 ∇ • (ρu) + ∆x 4c [∂ m (ρu 2 m ) + ∇ • (ρu m u)] + ρ w U • n (A11) E(x) = ± ρu t tan θ + O(∆xc ∇ρ ) + O (∆x ∇(ρu) ) + O ∆x ∇(ρuu) c s (B5)
For irregular boundary nodes associated with irregular cell pairs, they can be approximated as regular boundary nodes through two kinds of operations. Firstly, the links between neighbouring boundary nodes could be treated as virtual cut links without causing extra net mass leakage. For example, as shown in Fig. 2 The required correction for shifting f j from x + e i ∆t to x can be expressed as:

E s (x) = ∆x D ∆S∆t f j (x + e i ∆t) -f j (x) = ∆x D ∆S e i • ∇ f j + O(∆x 2 ) = ∆x D ∆S e i • ∇ f eq j + O(∆x 2 ) (B7)
where the expression of f by f eq (Eq. ( 9)) has been used. Omitting the high order term O(c∆x 2 ) and applying the analysis of e i • ∇ f eq i given in Eqs. (A3) to (A8), Eq. (B7) can be rewritten as: Interestingly, the averaged mass leakage of two cascaded serrated cells share a same formula with that of an ideal serrated cell (see Eq. (C1)). This similarity is important because it indicates that the mass leakage of two cascaded cells with an angle deviation can be directly described by that of a rotated cell given by Eq. (C4). Based on that, three cascaded serrated cells can be decomposed into one rotated serrated cell and two cascaded serrated cells with an angle deviation, and the total mass leakage is given by Eq. (C5). Recursively, a series of cascaded serrated cells with an angle deviation ∆θ (see Fig. 19(b)) can be estimated thanks to Eq. (C5).

Considering that the distance between the starting and ending "sawtooth" tip nodes (e.g. the "Start" and "End" points in Fig. 19 As a matter of fact, the extra mass leakage can be attributed to the associated velocity disturbances.

Along a curved boundary, the tangential and normal (relative to the virtual planar boundary, e.g.

AC shown in Fig. 3(d)) velocity disturbances at a boundary node due to the boundary curvature can be estimated as:

∆u t = u t (cos ∆ϕ -1)u n sin ∆ϕ, ∆u n = u t sin ∆ϕ + u n (cos ∆ϕ -1) (C10)

where u t and u n are the tangential and normal velocity components relative to the curved boundary, respectively. Applying Eq. C10, the mass leakage can be expressed as: In Eq. (C13), the angle variation ∆ϕ can be approximated as:

∆ϕ ≈ n∆x cos θ R c - ϕ 2 (C14)

FIG. 1 .

 1 FIG. 1. Sketches of the streaming process at boundary nodes considering a D2Q9 lattice: (a) boundary nodes with LB links irregularly cut by a general smooth boundary and (b) boundary nodes with LB links regularly cut by an aligned planar boundary. Nodes in the solid area, e.g. G 1 , are not necessary.

  FIG. 2. Irregular boundary nodes at a smooth boundary. The symbols are the same with those in Fig. 1(a). Lattice cells with boundary going through two perpendicular edges, e.g. B in (b) and C (c), are referred to as "irregular cells". Three typical cases are illustrated: boundary node irregularly cut with the cut-links similar to those cut by aligned boundaries in Fig. 1(b), without irregular cells involved (a); part of the boundary nodes are involved with irregular cells (b); and all boundary nodes are involved with irregular cells when a straight boundary is inclined at 45 • (c),

  ) are due to ρu m = ρ w U •n when n = n m , indicating that the mass leakage is mainly caused by the intrinsic deficiency of the wall-cut LB links in recovering the momenta. Meanwhile, these terms depends on the boundary orientation significantly. Particularly, for θ → 0, only the first kind of irregular boundary nodes emerges, and the dominating term ρu t tan θ → 0; For θ → 45 • , only the second kind of irregular boundary nodes emerges, and the dominating term ρu t tan(±θ ∓ 45 • ) → 0.

Figures 3 (

 3 Figures 3(a)-(d) display the four basic kinds of serrated cells, i.e. an ideal serrated cell with the two sawtooth tip nodes exactly on a planar boundary (see Fig. 3(a)), an ideal serrated cell shifted from the planar solid boundary by a distance ∆s (see Fig. 3(b)), an ideal serrated cell rotated from the planar boundary by an angle of ∆θ (see Fig. 3(c)), and an ideal serrated cell with the sawtooth tip nodes on a boundary with a curvature radius R c (see Fig. 3(d)).

  (a)-(d), the averaged mass leakage of an ideal serrated cell, and those caused by ∆s, ∆θ and R c are estimated as (see Appendix C for details):

  FIG. 3. Sketches of serrated cells including boundary nodes between two successive "sawtooth tips" (e.g. A and C ). Four basic types are shown: an ideal serrated cell with the two sawtooth tips (A and C) on a planar boundary (a), an ideal serrated cell shifted from a planar boundary by ∆s (b), an ideal serrated cell deviating from a planar boundary by an angle of ∆θ (c), and an ideal serrated cell with the sawtooth tips on a curved boundary. AC is a real (solid) or virtual (dashed) planar boundary. θ is the minimal angle between AC and the coordinate axes. u t and u n are the tangential and normal velocity components, respectively. In (d), R c is the local boundary curvature, ϕ is angle corresponding to the arc AC, and ∆ϕ is angle between AC and the local tangential direction.

FIG. 4 .

 4 FIG.4. Sketch of the interpolation based IB method. For interpolation of velocity, the projection point is also used.

FIG. 5 .

 5 FIG. 5. Sketch of the Poiseuille flow through a straight channel inclined at θ . U b is the translating velocity of both walls.

FIG. 6 .

 6 Fig. 7(a) displays the local mass leakage amplitude |E| max at each wall as a function of the inclined angle θ . As can be observed, |E| max at the two walls varies significantly with θ in a

Fig. 7 (

 7 Fig. 7(b) shows the averaged mass leakage Ē as a function of θ . The linear part | Ēl -Ēu | approaches zero at θ = 0 • and 45 • , and displays significant oscillations over 0 • < θ < 45 • . This indicates that | Ēl -Ēu | is dominated by the term O(ρu t ∆x L ) which is zero at tan θ = 1/N, e.g. 0 • , 26.6 • and 45 • (see §V). In contrast, the net mass leakage Ēl + Ēu generally increases with θ , indicating that the term O(∆xρu t ∂ n u t ) has a positive reliance on θ . The observable average

FIG. 7 .

 7 FIG. 7. Mass leakage as a function of the inclined angle θ . The sub-indices l and u indicate the lower and upper walls, respectively. In (a), |E| max is the local mass leakage amplitude. In (b), Ēl -Ēu and Ēl + Ēu correspond to the averaged mass leakage due to the linear terms in Eq. (20) and the net averaged mass leakage, respectively.

Fig. 8 (

 8 Fig. 8(a) and (b) display the mass leakage as functions of U b for θ = arctan 1 2 and 30 • , respectively. As can be observed, for both values of θ , the local amplitude |E| max at both walls is almost proportional to the boundary velocity U b , which can be related to the dominating terms ∝ ρu t in Eqs. (15) and(16). Meanwhile, the net mass leakage Ēl + Ēu in both figures exhibits similar

Figure 9 (FIG. 8 .

 98 Figure 9(a) and (b) display the computed streamwise velocity profiles and the mass leakage over 0 ≤ (u max -U b )/U re f ≤ 1, respectively. As shown in Fig. 9(a), profiles of the relative streamwise

FIG. 9 .FIG. 10 .

 910 FIG.9. Velocity profiles (a) and mass leakage (b) for 0 ≤ (u max -U b )/U re f ≤ 0.5. k is the value of ∂ n u t at the walls, and k re f is associated with u max -U b = U re f . In (b), Ēl + Ēu can be described by a parabolic curve.

  Fig. 11(a) and (b) display the local mass leakage and sectional mass flux variation, respectively.

FIG. 11 .FIG. 12 .

 1112 FIG. 11. Local mass leakage (a) and sectional mass flux (b). "NMC" is short for "no mass correction".

FIG. 13 .

 13 FIG.13. Sketch of Taylor-Couette flow through a circular channel (a) and the adopted grid distribution (b).

Figure 14 (

 14 Figure 14(a), (b) and (c) display the computed velocity profiles, the local and averaged mass leakage, respectively. As shown in Fig. 9(a), profiles of the relative velocity (u -Ωr 1 ) agree well with the analytical solutions. All the Taylor-Couette flow in the rest of this paper is accurately reproduced similarly without showing. For the mass leakage shown in Fig. 14(b)-(c), the local amplitude |E| max is almost proportional to r 1 . Considering r 1 is proportional to the tangential velocity u t = Ωr 1 and Ω is fixed, this observation can be explained by the leading terms (∝ ρu t ) in Eqs. (15) and(16). Meanwhile, the linear averaged mass leakage ĒL exhibits significant oscilla-

(u-Ω r 1 )FIG. 14 .

 114 FIG. 14. Velocity profiles (a), local mass leakage (b) and averaged mass leakage (c) for 1 ≤ r 1 /h ≤ 6.

Figure 15 (

 15 Figure 15(a) and (b) display the local and averaged mass leakage over 0.125 ≤ Ω/Ω re f ≤ 1, respectively. As can be observed, the local amplitude |E| max (Fig. 15(a)) and the linear averaged mass leakage ĒL (Fig. 15(b)) are almost proportional to Ω. The former can be well explained by the leading error terms (∝ ρu t ) in Eqs. (15) and (16), and the later is consistent with the term

4 E L ×10 3 (

 43 FIG. 15. Local (a) and averaged mass leakage (b) for 0.125 ≤ Ω/Ω re f ≤ 1.

r 1 and

 1 Ω are set equal to h and Ω re f , respectively.

Figure 16 (

 16 Figure16(a) and (b) display the local and averaged mass leakage over 0.0025 ≤ ∆x/h ≤ 0.04, respectively. As can be observed, the local amplitude |E| max (Fig.16(a)) is slightly affected within 2%, which is consistent with Eqs. (15) and (16) where the ∆x-related terms are of relatively minor amplitude. In contrast, the net mass leakage ĒNL (Fig.16(b)) is almost proportional to ∆x, and the linear averaged mass leakage ĒL (Fig.16(b)) exhibits a similar trend with some oscillations.

2 E NL ×10 3 (FIG. 17 .

 2317 Fig.17(a) and (b) display the local mass leakage along the inner cylinder and time history of the total mass M, respectively. As shown in Fig. 17(a), the local mass leakage varies significantly

  FIG. 18. Comparison of pressure coefficient profiles (a) and the corresponding contour lines with LMC (b) and AMC (c).

E

  (x) = E(x, ±θ )∆x D-1 cos θ + E(x, ±θ ∓ 90 • )∆x D-1 cos(±θ ∓ 90 • ) ∆x D-1 cos θ + ∆x D-1 cos(±θ ∓ 90 • ) =ρu t tan(±θ ∓ 45 • ) + O(∆xc ∇ρ ) + O (∆x ∇(ρu) ) + O ∆x ∇(ρuu) c s (B6)In addition, using shifted links to derive Eq. (B3) requires extra mass leakage correction terms.

  E s (x) =O(∆xc ∇ρ ) + O (∆x ∇(ρu) ) + O ∆x ∇(ρuu) c s (B8)Since all the terms in Eq. (B8) have been included in Eq. (B6), Eq. (B6) adequately describes the resultant mass leakage at the irregular boundary nodes associated with irregular cells.E sc,s = ∆s∂ n (ρu n ) (C2)Thirdly, extra mass leakage of a rotated serrated cell caused by an angle deviation ∆θ (see Fig.

3 1 1 2 1

 12 (c)) is quantified. The mass leakage can be estimated based on the virtual ideal serrated cell ABC with the normal velocity of the virtual boundary AC being:u n,v = -u t sin ∆θ + u n cos ∆θ (C3)Combined Eq. (C3) with Eq. (C1), the mass leakage of a rotated serrated cell can be expressed as:E sc,r = (E sc,iρu n,v ) AC AC cos ∆θ + ρ w U • n =ρu t tan ∆θ + O(∆xc ∇ρ ) + O (∆x ∇(ρu) ) + O ∆x ∇(ρuu) c s (C4)Since the angle deviation ∆θ may vary significantly for each serrated cell, it is necessary to estimate averaged mass leakage of cascaded rotated serrated cells.Firstly, two cascaded serrated cells are considered. As shown in Fig.19(a), averaged mass leakage of two cascaded serrated cells can be estimated from those of each cells (see Eq. (C4)) as:E sc,2rt = (ρu t tan ∆θ 1 )| x=x C AA tan ∆θ AA (tan ∆θ 1 + tan ∆θ 2 ) -(ρu t tan ∆θ 2 )| x=x A AA tan ∆θ AA (tan ∆θ 1 + tan ∆θ 2 ) + O(∆xc ∇ρ ) + O (∆x ∇(ρu) ) + O ∆x ∇(ρuu) c s =O(∆xc ∇ρ ) + O (∆x ∇(ρu) ) + O ∆x ∇(ρuu) c s (C5)

  FIG. 19. Sketches of cascaded serrated cells at a smooth boundary: two cascaded serrated cells (a) and a serious of cascaded serrated cells (b).

Et (cos ϕ 2 - 1 ) -u n sin ϕ 2 ]

 212 sc,c = E(x C )∆s C + ∑ B P=A E(x P )∆s P R c ϕ = E sc,ct + E sc,cn ct = ρ∆u t cos θ ∆x -∑ tan -1 θ n=1 ρ∆u t sin θ ∆x R c ϕ , E sc,cn = -ρ∆u n sin θ ∆x -∑ tan -1 θ n=1 ρ∆u n cos θ ∆x R c ϕ (C12)where E sc,ct and E sc,cn are the mass leakage caused by the tangential and normal velocity disturbances, respectively, and terms in the order of O(∆x 2 ) are omitted. Substituting Eq. (C10) into Eq. (C12), one obtains :cos θ ∆x -∑ tan -1 θ n=1 ρ[u t (cos ∆ϕ -1)u n sin ∆ϕ] sin θ ∆x R c ϕ , E sc,cn = -ρ[u t sin ϕ 2 + u n (cos ϕ 2 -1)] sin θ ∆x -∑ tan -1 θ n=1 ρ[u t sin ∆ϕ + u n (cos ∆ϕ -1)] cos θ ∆x R c ϕ (C13)

TABLE I .

 I Examples of mass leakage observation using the BGK lattice Boltzmann model. "FH" refers

  (a), by considering the virtually cut links e 5 and e 6 between B 2 and B 3 , the node B 3 can be treated as an approximated regular boundary node with a horizontal main direction. Secondly, wall-cut LB links at one boundary node can be shifted to its neighbouring boundary node to complete the remained approximation to regular boundary node. For example, as shown in Fig.2(a), by shifting the cut links e 3 and e 4 at B 1 to B 2 , the total cut links, including the virtual e 5 and e 6 , exactly approximate B 2 as two regular boundary nodes with different main directions (one is horizontal and the other is vertical). Mass leakage at the approximated regular node pairs, e.g. those at B 2 as analysed, can be estimated by Eq. (B3) with γ being ±θ and ±θ ∓ 90 • , respectively, i.e.

ACKNOWLEDGEMENTS

This work was funded by ANR Industrial Chair ALBUMS (grant ANR-18-CHIN-0003-01).

Appendix B: Local mass leakage at irregular boundary nodes

As depicted in §IV, irregular boundary nodes can be classified into two kinds regarding whether they are associated with irregular cells. Their mass leakage are analysed separately now.

For irregular boundary nodes not associated with irregular cells ( e.g. B 1 shown in Fig. 2(a)), they can be directly approximated as regular boundary nodes because the cut links are exactly the same as those cut by aligned boundaries (see B 1 shown in Fig. 1(b) and Fig. 2(a)). The main difference is that the main direction n m of the approximated regular boundary node is different from the local normal vector n, i.e.

where γ is the angle rotating from n to n m anti-clockwise, and n t is the tangential unit vector.

By substituting Eq. (B1) into Eq. (A11), the mass leakage can be estimated as:

which can be simplified as:

where Firstly, mass leakage of an ideal serrated cell is quantified by considering the one shown in Fig. 3(a) without loss of generality. Apparently, the averaged process is conducted over the boundary nodes along AB (excluding A) and the one at C. According to the analysis proposed in §IV B, the irregular boundary node at B and C can be approximated as two regular ones with a horizontal and vertical main directions, respectively, and those along AB (except B) can be approximated as regular ones with a vertical main direction. Accordingly, their mass leakage can be estimated by Eq. (B3) with γ being equal to 90 •θ and -θ , respectively. Consequently, the averaged mass leakage of an ideal serrated cell can be expressed as:

where the macroscopic variables are evaluated at C, and the terms in the order of O(∆x 2 ) are omitted.

Secondly, extra mass leakage of a serrated cell (see Fig. Consequently, the extra mass leakage caused by ∆s can be expressed as:

Substituting ∆x R c = 2 sin θ sin ϕ 2 (Eq. (C8)) and sin ϕ 2 = ϕ 2 + O( ϕ 3 8 ) into the above expression, it comes :

Now substituting (C15) along with sin

), sin ∆ϕ = ∆ϕ + O(∆ϕ 3 ) and cos ∆ϕ = 1 -∆ϕ 2 /2 + O(∆ϕ 4 ) into Eq. (C12), the mass leakages E sc,ct and E sc,cn can be simplified as:

Consequently, by omitting the terms in the order of O(

), the dominating parts of the total mass leakage caused by boundary curvature can be estimated as:

Notably, the above analysis of averaged mass leakage is derived from the mesoscopic view to the macroscopic view, and thus is generally valid for LB simulation. Whereas, it should be noticed that it is not the only way to quantify the averaged mass flux. For example, Ginzburg [START_REF] Ginzbourg | Local second-order boundary methods for lattice Boltzmann models[END_REF] proposed an exact computation of mass leakage in inclined channel Stokes flow based on distribution functions reconstructed by Chapman-Enskog analysis up to second order, and it has been already exemplified analytically for stair-wise boundaries.