{'fraction_non_alphanumeric': 0.03199523052464229, 'fraction_numerical': 0.012370826709062003, 'mean_word_length': 4.624196702989662, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Flow force calculation in Lattice Boltzmann Method

We revisit force evaluation methodologies on rigid solid particles suspended in a viscous fluid and simulated via lattice Boltzmann method (LBM). We point out the non-commutativity of streaming and collision operators in the force evaluation procedure and provide a theoretical explanation for this observation. Based on this analysis, we propose a discrete force calculation scheme with enhanced accuracy. The proposed scheme is essentially a lattice version of the Reynolds Transport Theorem (RTT) in the context of the lattice Boltzmann formulation. Besides maintaining satisfactory levels of reliability and accuracy, the method also handles force evaluation on complex geometries in a simple and transparent way. We run simulations for NACA0012 airfoil for a range of Reynolds numbers ranging from 100 to 0.5x10^6 and show that the current approach significantly reduces the grid size requirement for accurate force evaluation.

simple and transparent way.
We run simulations for NACA0012 airfoil for a range of Reynolds numbers ranging from 100 to 0.5 × 10 6 and show that the current approach significantly reduces the grid size requirement for accurate force evaluation.

I. INTRODUCTION
The lattice Boltzmann method (LBM), a discrete space-time kinetic theory, has made major leaps in solving hydrodynamic problems at low Mach numbers [1][2][3][4][5][6][7][8][9][10][11]. The theoretical foundation of the algorithm is well established and boundary condition,s at the level of discrete populations, are well developed, with a number of variants ranging from the simple and efficient bounce-back (BB) boundary condition for stationary walls to the microscopic diffusive boundary condition [12]. An important asset of the method by now, is its ability to deal with complex shapes and moving boundary problems in a simplified but efficient manner.
For such systems involving fluid-solid interaction, along with appropriate treatment of boundary conditions, an accurate calculation of force (lift and drag) or torque, on the solid body is often crucial.
In order to compute these hydrodynamic forces on objects, two widely used methodologies are: the stress integration (SI) approach [13][14][15] and the momentum exchange algorithm (MEA) [16].
In the stress integration method (as the name suggests), one computes the force by integrating the stress tensor on the surface of the body. In the momentum exchange algorithm (MEA), the force is computed by accounting for microscopic exchange of momentum at the wall between fluid and solid boundary, directly in terms of the discrete probability density function. Though both methods are well established in general, the momentum exchange method is shown to be more accurate than the stress integration method [17] at moderate to high Reynolds numbers.
The momentum exchange method allows computation of forces and torques on solid bodies in a simple and an effective manner by calculating the exchange of momentum between the solid surface and the fluid. Essentially, the momentum exchange method exploit the fact that force exerted on the body should be equal to the difference between momentum carried by populations going into the wall and those coming back from the wall.
In its most basic version (see Fig. 1), the solid boundary is approximated in a discrete sense at the middle of every fluid-solid link, each of which crosses the boundary and connects a fluid node. Even in this case, when the boundary of solid object is laid down approximately, MEA is shown to be quite effective [16,18]. For low resolution simulations, the disadvantage of a cartesian grid in representing a complex shape is often resolved either by interpolations near boundary [19] or by having multi-resolution frameworks.
In this manuscript without loss of any generality, we adhere to the simple and widely used BB treatment of the complex boundary and point out a source of ambiguity in the force computation routine. Essentially, the issue is related to the fact that LBM is a repeated sequence of collision and streaming. In absence of a boundary, the order of these operations does not matter. However, in the presence of a solid boundary, this symmetry breaks and collision and streaming operations, no longer commute with each other. Streaming is the generator of spatial translations while the boundary operator by construction must break translational invariance. Hence, in the presence of solid boundaries, the correct order of operations to calculate the flow forces, is not self evident. We describe two possible sequence of operations indicated by algorithm A and algorithm B in Fig. 1. In Fig. 2,3, we show the effect of these two approaches on two test cases, that is, flow past a 2D circular cylinder and flow past a 3D NACA0012 airfoil, at moderate Reynolds numbers. These are standard non trivial test cases extensively used for validation of numerical schemes and good amount of benchmark experimental and computational data are available for the same [20][21][22][23]. To show that the discrepancy in drag values between the two algorithms, is generic, the force is computed with two discrete velocity models (D2Q9 and RD3Q41). The use of algorithm B is noticed in some literature [16,24], however, to the best of our knowledge, a detailed comparison between the two algorithms and an analysis as to why one is better than the other, is not presented in any work.
In this work, we investigate the difference in the two approaches, theoretically and numerically. In order to understand the differences, we formulate a discrete analogue of the Reynolds Transport Theorem (DRTT) in the discrete space-time setting of the lattice boltzmann equation and devise a simple and effective way for flow force calculation. The main idea in this approach is that the complex-surfaced solid can be enclosed inside a simpler cartesian-friendly bounding box, creating a control volume that facilitates balancing of momentum fluxes and accurate calculation of surface stresses on the solid boundary. The method simplifies to MEA when the outer encapsulating surface falls directly onto the discrete solid surface [25]. Hence, DRTT can be thought of as a modified MEA approach, built to tackle the errors that arise from the inability of cartesian grids to resolve complex surfaces, especially at low resolutions. To demonstrate this approach, we consider a 1D toy model as shown in Fig. 4. The standard lattice Boltzmann equation (in one dimensional space) with a relaxation term is a two step evolutioñ A more detailed description of notations used in lattice Boltzmann literature is given in the next section. It would be convenient for further analysis to define the global momentum calculated over the control volume as, The evolution equation for global momentum derived from Eq.(1), for a semi-infinite control volume with only X 1 as the boundary point, simplifies to : where, F 1 (X 1 , t+1) is the contribution from the boundary condition in a general case. In the simple case of a bounce-back boundary, if the boundary returns an un-collided distribution (F 1 (X 1 , t + 1) = f −1 (X 1 , t)), Eq.(3) simplifies to,  [22].
, where X 1 is the boundary grid point of the control volume and X 2 , X 3 , X 4 are the bulk grid points. In the toy example, we assume that the grid is unbounded on the right end and continues to X 5 , X 6 and so on.
simplified version of RTT, where,f i represents a post-collision population.
To computationally validate our claims, we present simulations and convergence studies for NACA0012 at zero angle of attack for Reynolds numbers ranging from 10 2 to 0.5 × 10 6 .
The outline of this paper is as follows: In Sec. II, we give a brief description of LBM. In Sec. III, we review widely used methods for flow force calculation in LBM literature, namely stress integration and momentum exchange. In Sec. IV, we formulate a discrete analogue of the Reynolds tranport theorem to understand the reason behind the accuracy of Algorithm B. In Sec. V, we demonstrate some numerical simulations to validate our findings.
where, D signifies the dimension of the setup. For these discrete velocity models, typically the kinetic equation is written with a single relaxation (BGK) term as, where, τ is related to mean free time and f eq i (M Slow ) is the discrete analog of Maxwell-Boltzmann distribution with M Slow being, the set of slow hydrodynamic moments of populations. In continuous kinetic theory, the slow moments consists of mass, momentum and energy density (M Slow = {ρ, ρu, ρu 2 + 3ρθ}), while in the discrete case one can often ignore energy conservation to focus on isothermal hydrodynamics only (M Slow = {ρ, ρu}). The lattice boltzmann equation is obtained by further discretizing Eq.(7) in space and time, as a sequence of discrete collisions followed by free flight(streaming) as, where β = δt/(2τ + δt). At this point, it needs to be pointed out that the same discrete evolution can also be written as free flight followed by collision aŝ It is straight-forward to see that the two versions are identical in absence of any solid boundaries. The main theoretical ingredient ensuring accurate hydrodynamics is the construction of the discrete equilibrium. Starting with a zero velocity equilibrium (w i > 0), a second order series approximation to the equilibrium is often written as, Here, it is important that the zero velocity equilibrium at a fixed temperature is positive (w i > 0) and satisfies the condition that lower order moments are same as Maxwell-Boltzmann. In particular, These conditions on zero velocity equilibrium ensures, Later with higher order lattices, the method was extended for finite, but subsonic Mach number case, by ensuring that the contracted sixth order moment [26,27] is, This allows one to write the equilibrium for a small departure from reference temperature, at least at zero velocity (second order in ∆θ) as, which gives a second order (in velocity) approximation to the equilibrium, The definition ofR can be found in the cited reference [26]. The weights and discrete velocity set for two used models (D2Q9 and RD3Q41) are provided in Table I and Table II. The energy shells for both the models is depicted in Fig. 5.
Arguably, the most crucial element of the lattice Boltzmann method is the boundary condition at a solid-fluid interface. Even though, the macroscopic boundary conditions are quite straight-forward, the mesoscopic boundary conditions on the discrete population(f i (x x x, t)), require some explanation. Historically, the two commonly used boundary conditions are bounce-back(BB) and diffusive boundary conditions [6,16,[28][29][30]. In recent years, hybrid formulations like the diffused bounce-back boundary condition [31,32], aimed at expanding the scope of application to a variety of problems, are also being used. In order to keep the message of this paper clear, we stick to the universally used BB boundary condition proposed by Ladd [16]. The basic implementation of BB boundary condition, assumes that a molecule hits the wall and reverses its direction.  Fig. 1). The distribution for the fluid points near the boundary is given by the boundary condition [16], whereĩ is the direction opposite to i and u u u w is the velocity of the solid wall. Here,f represents the post-collision population. There are several higher-order interpolation based curved boundary conditions in literature that address the drawbacks in accuracy of the midway bounce-back boundary condition when dealing with curved surfaces [19]. However, the essential ideas of this paper remain unchanged even with higher-order interpolation based schemes and can be transferred to the same.
where the pressure term p = ρθ (θ = θ 0 for isothermal models) can be computed easily at every grid point. However, the deviatoric part (σ αβ ) of the stress tensor involves the velocity gradient tensor ∂ α u β , that needs to be approximated via finite differences which typically introduces additional error. This inaccuracy in calculation of velocity gradient tensor is circumvented in LBM by observation that the deviatoric stresses is evaluated as the second moment of the non-equilibrium part of the population [33] However, LBM being a cartesian grid based method, another uncertainty in force estimation is added by the calculation of a surface normal (n n n) for complex geometries. Thus, for a typical flow setup the SI calculation results have lower accuracy in contrast to the alternate momentum exchange method [17].
A computationally efficient and simple method is the momentum exchange (ME) algorithm [16] where, the total forces are computed as pairwise sum of momentum difference in discrete directions as populations bounces back from body surface near boundary points.
In its basic version (see Fig 1), boundary is approximated quite effectively in a staircase manner at the middle of every fluid-solid links, each of which crosses the boundary and connects a fluid node. The total force experienced by the solid object is given by, where,ĩ is the direction opposite to i and the summation runs over all fluid points next to discrete boundary. The ambiguity pointed out in the introduction of the paper, regarding the choice of algorithm (see Fig.1), plays an important role here. The momentum exchange algorithm is a limit to the more general Reynolds Transport Theorem [25]. In the next section we derive the discrete version of the Reynolds Transport Theorem and arrive at simplistic force evaluation routine for complex shaped objects.

IV. DISCRETE REYNOLDS TRANSPORT THEOREM
Consider a complex shaped solid object with a continuous boundary marked by ∂Ω as in Fig. 6, placed on a Cartesian grid. The solid is surrounded by fluid nodes (Ω) and the discrete analogue the boundary is represented by ∂Ω D . The total momentum for all the fluid grid nodes in the bulk and boundary, is given by, The momentum balance for the fluid points as the system evolves from t to t+∆t, implies where, J G α (t + ∆t) and J L α (t) represents the gained and lost momentum at the boundary [25], as seen in Fig.6. It must be noted that J G α is calculated over the boundary grid nodes where population is added to the system and J L α is calculated over the boundary grid nodes where population leaves the system, when the system displaces from t to t + ∆t. The gained momentum is simply written as, where,f i denotes a collided population. As per the discussion in the introduction section of this article, there are two ways of calculating J L α (t), based on the choice of algorithm (A or B), The choice of Algorithm B simplifies Eq.(21) to give us the complete momentum balance for fluid nodes as, The boundary surface contribution on RHS of the Eq. (24), is exactly the momentum exchange algorithm proposed by Ladd [16], whereas the choice of Algorithm A would lead to a slightly more complicated expression. Using Newton's second law, one can estimate the forces exerted by the fluid, on submerged solid objects as, The first term in the force formula is dictated by the boundary condition used in the algorithm. For a simple bounce-back boundary condition, the force formula simplifies to, We extend the derivation carried out for a single boundary to a system with two boundaries as given in Fig.7, where the fluid nodes are given by Ω, and are sandwiched between the two boundaries ∂Ω 1 and ∂Ω 2 , such that the momentum balance gives, The construction of the outer boundary ∂Ω 2 is done in such a way that n n n 2 and location of ∂Ω 2 is trivial to resolve on a lattice. A rectangular bounding box aligned with the grid is a good example of the same. This construction helps us bypass calculations at the solid surface (∂Ω 1 ) by rearranging the Eq.(28) to, where, F α (t) is the force exerted by the fluid on the surface ∂Ω 1 . We label this force evaluation method as DRTT and conduct simulations in the next section to validate the method.

V. RESULTS
As a first example, we consider the flow past a circular cylinder in two dimensions, which is a regularly used benchmark case in computational fluid dynamics. The most basic calculation of the pressure coefficient on the surface of the cylinder involves calculation over boundary fluid nodes without any kind of extrapolation. The pressure coefficient at one boundary fluid node is given by,  [20]. Fig. 8a) shows the coefficient of pressure (C p ) vs θ for Re= 100 case. Fig. 8b) compares the drag coefficient (C d ) as a function of number of grid point per diameter for momentum exchange(ME), stress integration(SI) and discrete RTT methods (DRTT). Fig. 8 shows the coefficient of pressure vs θ for Re= 100 case. We see that as expected, the even at a low resolution of D = 15, the C p curve is a close match with literature. Once the correctness of C p is established, the flow force is calculated using all the three methods discussed in the previous section. The drag coefficient over a circular cylinder of diameter D is given by, where, x is taken to be the direction along the flow. The force term (F x ) is calculated using popular methods (MEA and SI), along with with the discrete Reynolds transport theorem (DRTT) algorithm. We do the convergence test for all the methods for Re= 100 and the results are given in Fig. 8. The aspect ratio for the computational domain is taken to be 12 : 1 with the length and width of the domain as 60D and 5D. The cylinder of diameter D is placed at 15D from the inlet. We can clearly see that SI does the worst among all the methods at a given resolution and henceforth we shall not discuss it in the following sections on flow past a NACA airfoil.
The simulations are run for a range of parameters, in order to study the dependence of the drag coefficient on the Reynolds number. As seen in Fig. 9, DRTT does much better than MEA at a given resolution. One can also notice that the difference between DRTT and MEA becomes more prominent at high Re where the flow starts to become unsteady. The airfoil of interest in our work is the NACA0012 airfoil, that has been extensively used for validation cases for turbulence models. NACA0012 airfoil is a skewed geometry with different grid requirements in x and y directions. This makes it an important case study for algorithms like lattice Boltzmann that use uniform meshes (∆x = ∆y). In order to validate our code, we reviewed the coefficient of pressure (C p ) measured over the airfoil surface at a Reynolds number of 1000 and three angles of attack (α) [23]. Fig. 11 shows us the match with literature. We plot for angle of attack= (7,9,15)  The lift coefficient (C l ) is given by, We can see that MEA and DRTT do almost equally well at low angles of attack. But as the angle of attack rises and the flow begins to separate, DRTT is able to more accurately calculate the drag and lift values. We attribute this to the fact that DRTT, being a volumetric method, does not have to resolve the complex boundary of the airfoil at low grid resolution. Reynolds number in order to accommodate a developing flow in the z direction at very high Reynolds. Table III gives details of the computational grid requirements for a converged result. The resolution demands (measured in points per chord) for a converged result vary as a function of Reynolds number. We assign an error percentage of 2 percent in order to call a result converged and tabulate the results in Table III. The converged C d vs Re curve as seen in Fig. 13 shows a good match with literature. Fig. 14 shows a close match with a Re 1/2 scaling which can be explained using boundary layer theory.

VI. CONCLUSIONS
In this work we revisited the force evaluation methodologies in the lattice Boltzmann

VIII. APPENDIX
In this section we present a basic implementation of the discrete Reynolds transport theorem for flow past a two dimensional object. The size of the control volume (Ω as seen in Fig. 7), for the simple case of 2 dimensional cylinder, is taken to be 1.5D × 1.5D and the corners are marked as cv x1 , cv x2 , cv y1 , cv y2 .