Moment‐based boundary conditions for straight on‐grid boundaries in three‐dimensional lattice Boltzmann simulations

In this article, moment‐based boundary conditions for the lattice Boltzmann method are extended to three dimensions. Boundary conditions for velocity and pressure are explicitly derived for straight on‐grid boundaries for the D3Q19 lattice. The method is compared against the bounce‐back scheme using both single and two relaxation time collision schemes. The method is verified using classical benchmark test cases. The results show very good agreement with the data found in the literature. It is confirmed from the results that the derived moment‐based boundary scheme is of second‐order accuracy in grid spacing and does not produce numerical slip, and therefore offers a transparent way of accurately prescribing velocity and pressure boundaries that are aligned with grid points in three‐dimensional.

boundary method that does not have the numerical slip. This drawback was first discovered by Ginzbourg and Adler 20 who quantified the effective location of the straight/diagonal boundaries in Poiseuille flow. He et al 21 derived explicit expressions for the slip velocity for different schemes in Poiseuille and Couette flows. Later, Bennett 22 gave an interpretation of some existing lattice Boltzmann boundary conditions for the distribution functions in terms of their moments. It was shown that on-node bounce-back has an entirely equivalent formulation of setting u n = Q nnt = Q ntt = 0 at a flat boundary aligned with gridpoints. Here, u is the fluid velocity (the first moment), Q is the third moment (and is a nonhydrodynamic moment), and subscripts n and t refer to normal and tangential components, respectively. Thus, it is seen that this method does not explicitly impose a constraint on the tangential component of velocity, leading to the potential for the artificial slip error. The TRT method may be used to enforce compatibility between the conditions on the third moment with the expected solution.
Extrapolation/interpolation schemes can be viewed as separate types of boundary methods. They require additional points to set constraints for velocity and pressure or their gradients. First proposed by Chen et al 15 the extrapolation scheme for flat walls was later modified and its numerical stability improved by Zhao-Li et al 23 using the idea of extrapolating the nonequilibrium parts of the distribution function at a wall node. Like the interpolation schemes, 24, 25 Guo's nonequilibrium extrapolation scheme can also be applied to curved walls. 26 It maintains simplicity, second-order accuracy and good numerical stability.
Another approach uses the idea of the bounce-back of the nonequilibrium parts of the distribution function, first proposed by Zou and He. 16 This approach, which was originally formulated in terms of distribution functions, also has for the D2Q9 lattice an equivalent interpretation in terms of moments: it can be viewed as setting the constraints u n = u t = Q ntt = 0 at solid boundaries. We see an explicit condition on the tangential velocity but it is more difficult to give a mathematical or physical justification for the third constraint. The extension of nonequilibrium bounce-back to the 3D D3Q15 lattice was first proposed by Zou and He 16 to handle pressure and velocity boundaries. It was then generalized by Hecht and Harting 27 allowing an arbitrary angled flow at the open boundary. Just like Zou and He, they used tangential momentum corrections, introduced by Maier et al, 14 at the inlet of the D3Q19 model. Their method is second-order accurate, explicit, spatially local, relaxation time-independent, and it allows to implement exact boundary conditions on the wall which lies exactly on the lattice nodes. However, one might find this method to be overcomplicated-applying the nonequilibrium bounce-back rule and then modifying the tangential distribution functions by introducing the transverse momentum corrections, recalculating the resting particle distribution function at the pressure boundary edges. Furthermore, because it uses a variation of the bounce-back rule, the imposed conditions for the velocity moments still lack justification from a physical sense. Performing an analysis of the moments at the face boundary reveals that the closure is somewhat arbitrary. The conditions are set for all three momentums, u n = U n , u t1 = U t1 , u t2 = U t2 , and two third-order moments, Q nt1t1 = U n ∕3 and Q nt2t2 = U n ∕3, where U n ∕3 is the equilibrium approximation, see (7). While the selection of the former three moments makes sense, the last two conditions seem odd because they both state the same and they do not have a clear physical interpretation like the second-order moments, momentum flux and stress, for example.
High-order accurate boundary conditions have been proposed by Ginzbourg et al. [28][29][30][31] Local second-order boundary (LSOB) method 28 is derived from the locally known distribution functions via second-order Chapman-Enskog expansion and Dirichlet boundary conditions with a given momentum. To the best of the authors' knowledge, this remains the only local boundary condition implementation method for the LBM that computes the exact solution of Poiseuille flow for arbitrary channel wall inclination with respect to the lattice. Despite it offering a highly accurate and fairly general treatment of boundaries, the LSOB method is often overlooked due to its relatively complex implementation. Multireflection (MR) boundary method, 30,31 however, is implementation-wise much simpler than the LSOB method and is formally third-order accurate for general flows. Derived for arbitrary shaped boundaries, the MR method has been specified for corner treatment, 31 complex porous flow, 32 slip conditions, 33,34 and extended for boundary/interface in advection-diffusion equations. Downsides of the MR method include nonconservation of the local mass and nonlocality. Local on-grid boundary conditions are important for efficient parallel computations. To confront the aforementioned drawbacks, we have taken the first steps of generalizing the moment-based approach by extending it to the 3D lattices with boundaries aligned with grid points, where the boundary conditions are imposed directly onto the hydrodynamic moments of the LBM.
The first hydrodynamic moment-based scheme for velocity boundaries was proposed by Noble et al 13 who used hydrodynamic moments, more precisely the velocity and energy 35 to solve for the unknown distribution functions at the boundaries. Their motivation was simple and valid-the bounce-back boundary condition has a relaxation time-dependent slip, and it cannot be easily generalized for mass inflows or moving walls. They were employing the F I G U R E 1 D3Q19 lattice. Rest particle velocity is not shown [Colour figure can be viewed at wileyonlinelibrary.com] hexagonal D2Q7 lattice where only two distributions functions are unknown on the boundaries aligned with grid points, and two conditions are required to solve for them. Other commonly used lattices, such as D2Q9 and D3Q19 (see Figure 1) have more unknown functions, and therefore they require more conditions that come from linearly independent moments and a more general approach.
The more general moment-based method for imposing hydrodynamic boundary conditions was proposed recently by Bennett. 22 He used the fact that since there is a one-to-one linear mapping from the distribution functions to its moments (m = Mf), this mapping can be inverted (f = M −1 m). One can switch between moments m and particle distribution functions f very easily and therefore can impose a condition on all m to find all f. At a boundary, not all the moments are independent, but the idea is to impose conditions on linearly independent moments only and convert these into conditions on the distribution functions. The intention is to use as far as possible the hydrodynamic moments only because we are simulating hydrodynamics. In some cases, choosing the higher order velocity moments over the hydrodynamic ones might offer better stability, however, exploring the stability of the selected moments requires further research. Some other insights into stability could also come from References 31,36, which discuss numerical behavior of lattice Boltzmann boundary conditions at nodal points, or References 37,38 which discuss higher order contributions to the lattice Boltzmann stress at boundaries.
The LBM with moment-based boundary conditions has been successfully applied to various physical systems, where exact hydrodynamic boundary conditions must be employed. [39][40][41][42] Furthermore, its stability and accuracy has also been commented on briefly, 43 and some theoretical analysis has been performed, 37 but there is still room for an in-depth stability analysis. Studying stability of the present 3D method is left for future research. The results presented in these works show that the method is second-order accurate for velocity and pressure, which matches the accuracy of the LBM, and it recovers the solution to the NSEs in unidirectional channel flow where boundaries are aligned with grid points exactly on all grid sizes. Another important finding is that the moment-based method in combination with the SRT collision operator works very well in the region of low to moderate Reynolds numbers, but more sophisticated collision operators (eg, TRT or MRT) are preferable for flows at high Reynolds numbers.
At present, the moment-based approach assumes walls are straight and aligned with grid point. Under such conditions the method is: • exact-pressure, velocity or shear stress is specified directly on nodes, and thus satisfies the boundary conditions by construction. Bounce-back and diffuse reflection do not fully possess this property. 44 However, it is unknown at present if the exactness property of the moment-based method can be extended to irregular geometries.
• local-only the information from the boundary cell is used in the calculations (like Zou and He, 16 LSOB 28 ). The method can be parallelized easily for efficient computing. This is in contrast to the interpolation/extrapolation schemes. It should be noted that, at present, the application of the moment-based method is restricted to straight boundaries and thus these comments should not be taken out of context. The extension to more complicated geometries is a subject for further research.
• straightforward-the idea and the implementation of the method is relatively simple. The unknown distribution functions are calculated using the hydrodynamic moments that the boundary conditions are imposed on. The trickiest part might be finding a physical interpretation for the higher order moments that may need to be at edge or corner boundaries, see Section 3. Complexity wise it is not as simple as the bounce-back rule and lacks geometric flexibility as opposed to LSOB and MR meaning that it is currently limited to boundaries aligned with the grid points. However, it is conceptually simpler and more straightforward than other methods that involve a mixture of the bounce-back rule, hydrodynamic moments, momentum corrections and other modifications to the distribution functions.
Although we have not considered diagonal boundaries here we expect the methodology to be directly transferrable. That is, we can still find incoming distributions by imposing constraints on an equal number of linearly independent moments precisely and locally at grid points. However, a detailed consideration is reserved for future study.
The remainder of this article is organized as follows. In Section 2, the LBM is described for the D3Q19 lattice. In Section 3, the moment-based method for imposing boundary conditions is extended to three dimensions and conditions from this method are explicitly derived. The results from the validation of the method are presented in Section 4, and conclusions are drawn in Section 5. Derivation details of the method are given in Appendix A1. The moment-based formulation of boundary conditions for the D3Q15 model is given in the Appendix B1.

LATTICE BOLTZMANN METHOD
The LBM is a discretized form of the Boltzmann equation that governs the evolution of the particle distribution function via streaming and collision processes. 45 It can simply be written as where f i is the discrete particle distribution function and F i is the source term. The BGK approximation 46 is used to describe the collision operator. It assumes that the distribution functions relax to their local equilibria f eq i over a short collision time . The equilibria f eq i are prescribed to be where is the fluid density, w i and c i are, respectively, the weights and discrete particle velocities given in Table 1, u is the fluid velocity c = x∕ t is the lattice speed, with x and t being the space and time steps. The sound speed c s is a constant. The hydrodynamic variables, such as density, momentum and momentum flux can be calculated from the distribution functions by taking velocity moments as redefining the local momentum with the half forcing term.

From lattice Boltzmann to Navier-Stokes
The fact the NSEs are embedded within the LBM can be shown by performing a multiscale analysis. One such method that has been an integral part of all the major comprehensive literature on the topic of the LBM 47-50 is the Chapman-Enskog expansion, where f i is expanded formally about a small parameter as The second-order Taylor series expansion in time and space is performed on (1) without the source term and neglecting the higher order terms. Separating the different scales and taking moments of the LBM allows us to find order-by-order contributions to the nonconserved moments, and this eventually yields hydrodynamic partial differential equations. Only the equilibrium and first-order correction parts of are needed to furnish the NSEs and higher order contributions are considered to be negligible here. It is found that the equilibrium and nonequilibrium part of the momentum flux tensor can be expressed as and the equilibrium part of the third-order moment, which appears in the recovery of Π 1 and hence the NSEs, is Combining these expressions into the moment equations yields to order 2 the weakly compressible NSEs, where the pressure p = c 2 s and the kinematic viscosity = (2 − t) ∕6. The term weakly compressible means that the incompressible NSEs are recovered in the low Mach number limit.

Two relaxation time LBM
As the name suggests, the TRT collision operator uses two relaxation rates, one for the even order moments and one for the odd order moments. The relaxation time for the even moments is directly linked to the viscosity of the physical system, while the other is a free parameter that can be fine-tuned for optimal accuracy and stability. Since the truncation errors depend on the product of these two parameters, 51 one can effectively tune the LBM for optimal performance with a judiciously chosen value of the so-called "magic parameter": where + is the symmetric relaxation time that is related to viscosity and − is the antisymmetric time. While different values of Λ offer improved stability for different problems, Λ = 1∕4 appears to be the most stable choice according to Fourier advection-diffusion equation analysis. 51,52 Still, any fixed value of Λ assures the correct dimensionless scaling and independence of the relative errors on the viscosity. 51 Because all the velocity sets are symmetric, distribution functions can be paired up in terms of their velocities as c i = −c ı forming the so-called link. 36 Any link can be decomposed into its symmetric and antisymmetric parts as The rest population is its own opposite so it only has a symmetric part, f + 0 = f 0 and f eq+ 0 = f eq 0 , and a zero antisymmetric component, f − 0 = 0 and f eq− 0 = 0. Using the introduced components, the TRT postcollision distribution function can be written as The relaxation time related to the viscosity is the symmetric one: Due to symmetry the preparatory calculations in the collision step are performed only for one half of the populations. That makes TRT almost as computationally efficient as SRT, but with an improved control over stability and accuracy. 50 And although TRT can be used to correct the numerical slip in bounce-back schemes, in the moment-based method it is chosen purely for numerical stability.

MOMENT-BASED BOUNDARY CONDITIONS IN THREE DIMENSIONS
The D3Q19 model has exactly 19 independent moments, which are all listed in Equation (14), starting from the zeroth velocity moment, which is otherwise known as density, and going all the way to the third and fourth-order moments, whose physical interpretation are not as clear. These moments are used in calculating the incoming particle distribution functions at the local domain boundaries.
The information in the LBM can travel diagonally from one site to another. Therefore, apart from face boundaries we also need to separately consider edges, where two faces meet, and corners, where three faces meet.
There are five unknown distribution functions at every face boundary, nine unknowns at every edge boundary and 12 unknowns at every corner boundary. It means that five, nine, and 12 linearly independent moments are required at every face, edge, and corner, respectively, to solve for the unknown distribution functions. However, not all of the moments in Equation (14) are linearly independent at a boundary. In fact, they can be placed into groups of unique combinations of distribution functions for any given boundary whether it is at the face, edge, or corner boundary.
Next, the derivation process for each of these different cases will be described, distinguishing between velocity and pressure type boundaries, and the combination of them. Note that the pressure boundary considered here only allows the inflow or outflow normal to the face boundary. This is not to be confused with an open type boundary where the velocity can have arbitrary directions and the pressure distribution may not be uniform.

Face velocity
From Equation (14), five hydrodynamic moments are chosen to impose a conditions on a boundary face. If, for example, the west boundary is chosen, see Figure Table 2 do not contain the informations of the unknown functions. They are Π yz , Q yyz , Q yzz , and S yyzz . By looking at their respective expressions in (14) one can confirm that they do not contain the unknown functions of interest. The moments in a row are not linearly independent so only one moment can be picked from each row to impose a constraint on it and to solve the system for the unknown f i at the boundary. The aim is to pick hydrodynamic moments only and avoid selecting the higher order moments as much as possible because they do not have a clear physical meaning.
For the velocity boundary it is logical to select the three momenta, u x , u y , and u z . The remaining two moments are chosen to be Π yy and Π zz due to a natural physical interpretation (normal stresses) compared with the higher order moments. Now that there are five linearly independent equations for the five unknowns, the system can finally be solved. Before solving it, the momentum fluxes need to be defined. Using the first two terms from the Chapman-Enskog multiscale expansion, see Section 2.1, the momentum flux is approximated using the Chapman-Enskog expansion used to derive the NSEs from the LBM: Since the higher order contributions to the moments, which include temporal derivatives of conserved moments, are assumed negligible, replacing the terms in the above expressions with (5) and (6) gives For simple boundaries, such as velocity inlet, slip walls and no-slip walls moving with a constant velocity, the tangential velocity derivatives in (16), which come from the O( ) contributions to the moment, can be discarded giving the following expressions for the momentum fluxes at the velocity face boundary: Setting the velocities at the face boundary to U x , U y , and U z and using the selected moment expressions from Table 2, the solution for the unknown distribution functions can be obtained as, where the density is expressed through the consistency condition, which relates the density and the momentum normal to the west face boundary, as More detailed step-by-step derivation process is given in Appendix A.1.

Face pressure
Pressure boundaries require the density to be specified, leaving out the normal momentum as it is now an unknown moment, see Table 2. So, the only change from the velocity type boundary is the selection of the density, , in the first group. Other momenta, u y and u z , and momentum fluxes, Π yy and Π zz , remain unchanged.
Restricting the pressure inlet boundary to normal flow, the tangential velocities are set to zero. Due to the velocities being zero and their derivatives being zero, only the first terms in the momentum flux expressions remain from (16) giving Setting the pressure value at the boundary to p = 0 c 2 s = 0 3 in lattice unites ( t = x = 1), where 0 is being imposed, and solving the new system, gives the following unknown functions for the west face boundary: The normal velocity u x can be calculated from the set density 0 and the known distribution functions through the consistency condition as See Appendix A0.2 for more detailed solution.

Edge velocity
For the edge boundary, exactly nine linearly independent combinations are required to solve for the unknowns.   Table 3. Ideally one would like to pick the first nine appropriate moments and be done with it, however, the combinations appearing to be different are not all linearly independent. One can easily check that by looking at the rows 4, 9, and 10 in Table 3, for example. The matrix composed from these expressions is a rank-two matrix meaning that only two of the involved different row moments can be selected to specify a boundary condition.
The same simply noticeable restriction applies to the rows 7, 13, and 14. rank So, in a situation where there are more unknown combinations than unknowns, the linearly independent rows have to be selected prioritizing the physically interpretable ones.
By looking at the rank of the matrix consisting of the unknown combinations for the south-west edge boundary, it turns out that the first seven rows from Table 3 are all linearly independent. Other dependencies are less obvious. The row 8 is a linear combination of the rows 1, 2, and 3. rank The rows 9 and 10 have already been covered earlier. The rows 11 and 12 are a linear combination of the rows 1, 3, 5 and 1, 2, 6, respectively.
And finally, the row 15 is a linear combination of the rows 1, 5, and 6. rank Considering the available options from the analysis above, for the velocity boundary at the south-west edge, the first nine appropriate moments are the three momenta, u x , u y , u z , the momentum fluxes and shear stresses, Π xx , Π yy , Π zz , Π xy , Π xz or Π yz , and one higher order moment, Q xzz or Q yzz . There is still some freedom in selecting the moments to complete the system, however, no matter how the nine moments are chosen, having the higher order moment in the selection is inevitable. By basing the choice of the two moments on the symmetry of the components, meaning that either Π xz and Q yzz or Π yz and Q xzz are selected, the final system is written as F I G U R E 4 An example of two perpendicular pressure inlets [Colour figure can be viewed at wileyonlinelibrary.com] Similar to using the truncated approximation (16) for the momentum fluxes, Π xx , Π yy , and Π zz , and shear stresses, Π xy and Π xz , the higher order moment Q yzz is approximated using its equilibrium value (7), where the terms of order O(u 3 ) are neglected. This can be justified by the fact that only the equilibrium value of Q is used in the recovery of the NSE up to the second-order through the Chapman-Enskog analysis. The equilibrium approximation can be written as The system of the unknowns and its solution is included in Appendix A0.3.

Edge pressure
Specifying the pressure inlet at the edge is not straight forward. The conditions here have to agree with the ones at both adjacent faces, but that cannot be achieved due to the uncertainty of the velocity values. At the south-west edge, the normal pointing into the domain has components on x and y axis. The conditions for velocities on the west face read u x = unknown, u y = 0 and u z = 0, and on the south face they are u x = 0, u y = unknown and u z = 0. Problems arise when trying to merge these conditions. Tangential velocity is easy, u z = 0, but how can we know what value to set for the other velocities? Are they both unknown or both zero, or maybe one is unknown while the other is zero? This is not a common physical setup, in fact it is far from it. Rarely, if at all two pressure inlets are encountered being perpendicular to each other. One possible setup is shown in Figure 4 where two channels form the perpendicular pressure inlets. In this situation, there is a solid wall separating the two openings meaning that all the velocities are zero at that point. This is nothing else but a velocity edge discussed above in Section 3.3.

Edge pressure-velocity
The combination of the pressure and velocity type boundaries can be encountered when considering a simple channel flow with a pressure inlet or outlet, or both, for example. So, in situations where the two adjacent faces have different boundary conditions imposed on them, namely, velocity and pressure, the density together with the three momenta have to be specified simultaneously at the edge boundary. It means that , u x , u y , and u z are definitely selected from Table 3. Momentum fluxes Π xx , Π yy , and Π zz also get included. Only one of the shear stresses, Π xz or Π yz , can be selected because of (23). And only one of the third-order moments, Q xzz or Q yzz , is a viable option due to (24). Following the choice made earlier when talking about the velocity boundaries, the moments Π xz and Q yzz are selected to complete the system. The only difference from the velocity edge boundary is that the pressure is known. Because of this and (25), the shear stress component Π xy is left out of the selection. The chosen moments are given in (31), and the solution is given in Appendix A0.4.

Corner velocity
For the corner boundary, the number of the unknown distribution functions and therefore the number of the required linearly independent equations is 12. The low-south-west corner node is considered here. It means that the 12 unknown functions are f 1 , f 3 , f 5 , f 7 , f 9 , f 10 , f 11 , f 13 , f 14 , f 15 , f 17 , and f 18 , see Figure 5. Listing all the unknown combinations for the low-south-west corner boundary gives a total of 19 different combinations. They are shown in Table 4. Every moment has a unique combination of the unknown distribution functions so one is left with no choice but selecting the first 12 appropriate moments to impose the boundary conditions on them. See Appendix A0.5 for the description of the linear dependencies.
For the velocity boundary, the main thing is to pick the three momenta, u x , u y , and u z , followed by the momentum fluxes and shear stresses, Π xx , Π yy , Π zz , Π xy , Π xz , and Π yz . Then ideally choosing the third-order moments before considering anything else. Therefore, the last three fourth-order moments from Table 4, S xxyy , S xxzz , and S yyzz , are overlooked for now. There are enough moments to impose a boundary condition at the corner without them. So, not counting the density, the next ten moment combinations (rows 2-11) are linearly independent. Together with the rows 13 and 14 they make the basic complete system of equations ready to be solved, see Appendix A0.5 for details. Again, there is still some freedom left in choosing which moments will be included in the final system, but there is no clear reason why one would be chosen over the other. For instance, which is better out of the two in each case? Is it Q xxy or Q xxz , Q xyy or Q yyz , Q xzz or Q yzz ? One could probably argue that there are two mathematically explainable options. From a symmetry point of view, either Q xxy , Q yyz , and Q xzz are selected or Q xxz , Q xyy , and Q yzz make the cut. This is as far as the mathematical reasoning can take. Any further choices are left to be made subjectively.
The final system of the moment combinations includes the momenta u x , u y , and u z , the momentum fluxes and shear stresses Π xx , Π yy , Π zz , Π xy , Π xz , and Π yz , and three higher order moments Q xxy , Q yyz , and Q xzz . Alternatively, one can choose the other trio of the third-order moments. The details of the solution at the low-south-east corner boundary are given in Appendix A0.5.1.

Corner pressure-velocity
As discussed earlier, for a system with mixed type boundaries, both the density and velocity conditions have to specified. If density is to be included into the linearly independent moment selection then one of the moments must be left out. Velocities are set, which means that one of the shear stresses must be discarded. There is no clear sign for which is the odd one out as they form a closed symmetry group. But the choice has to be made so Π yz is discarded from the selection of moments, see Appendix A0.5.2 for more details.

RESULTS
In this section, the newly derived 3D moment-based boundary conditions are validated on various 2D and 3D benchmark cases testing the set conditions at the faces, edges and corners of the domain. The 2D simulations are used to verify our method, while the 3D simulations are used for the assessment of the method's accuracy. Velocity boundaries are tested for stationary and moving walls with the no-slip condition applied. The relaxation time + is defined via the Reynolds number, Re = UL∕ . The LB velocity is set to U = 1∕c = t∕ x = 0.1, L = N − 1 so that x = 1∕L and relates to + through (13). From (10), − is always chosen so that Λ = 1∕4 for the moment-based method in all examples. One of the tests in 2D is the relaxation time independence study performed using the 2D Poiseuille flow case. Moreover, it is linked to the exact recovery of the no-slip condition for the velocity on the wall. Because the solution of the developed 2D Poiseuille flow is essentially one-dimensional, a relatively small numerical grid of 33 × 3 × 3 can be used in the calculations to test for the dependence and the no-slip recovery. A grid size of 129 2 × 3 is selected for the 2D lid-driven cavity flow to match the meshes employed by other authors. 53,54 Velocity profiles along the centerlines as well as the extreme values of the stream functions are compared in order to validate the method.
For the 3D convergence studies, the grid size is being varied proportionally with the relaxation time fixing the Reynolds number. A changing grid of N x × N y × N z up to 129 2 × 3 is used for the 3D developed duct flow in 2D while fixing the

2D validation of the method
The relaxation time dependence study is performed to see how the moment-based boundary conditions compare to the bounce-back scheme using SRT. Although the -dependence of the bounce-back relative error is removed with any fixed Λ when using TRT, only Λ = 3∕16 locates the solid walls midway in straight Poiseuille flows. The results are shown in Figure 6. For a fair comparison, the duct flow is considered with the on-grid wall location for the moment-based method and midway location for the bounce-back scheme. The moment-based boundary conditions with either of the collision schemes show a whole range of the relaxation time values for which the solution only has accumulated machine precision type error and does not change in general, apart from the small region where approaches its asymptotic lower limit, min = 1∕2. The lower limit, min , simply restricts the fluid viscosity from becoming negative. Similar behavior to that of the moment-based boundary conditions can be seen for the bounce-back scheme with TRT, but not with SRT. This is the famous numerical slip error 20,21 coming from the uncertainty of the wall position that depends on the relaxation time.
Only when the wall is exactly in the middle between the fluid and boundary nodes, = √ 3∕16 + 1∕2, the numerical error vanishes resulting in a close-to-machine precision accumulated error.
A relaxation time dependence study confirms that the moment-based boundary conditions generate no slip at the walls, and they are -independent, see Figure 6. The channel width and Reynolds number are both fixed while adjusting the magnitude of the flow driving force when sweeping through the LB viscosity values. The relative slip velocity for the bounce-back scheme using SRT increases monotonically with crossing the horizontal axis at = √ 3∕16 + 1∕2 which corresponds to the middle position of the wall. The slip velocity is scaled by the maximum flow velocity in the duct giving the relative slip velocity.
The 2D lid-driven cavity flow is selected as one of the benchmark cases to test the ability of the moment-based boundary conditions to describe more complicated flows. The performed grid convergence study reveals that the combination of the LBM and moment-based boundary conditions is of second-order accuracy, as shown in Figure 7. At Re = 100, the TRT scheme yields slightly better accuracy than SRT (not shown) due to its improved stability. At Re = 1000, the TRT results still show second-order convergence, but relatively larger errors. In addition, no stable SRT results were obtained at Re = 1000. For comparison, the grid convergence results of the 2D moment-based boundary conditions obtained by Mohammed and Reis 43 using MRT and Reis 38 using TRT have been included. The results of Mohammed and Reis show better accuracy at Re = 100, but interestingly the results of the present 3D TRT-LBM moment-based boundary conditions at Re = 1000 show slightly improved overall accuracy in comparison to the 2D lattice models. Velocity profiles along the centerlines as well as the extreme values of the stream functions have been compared with the data from the literature 53,54 at Reynolds numbers Re = 100 and Re = 1000. The results shown in Figure 8 and Table 5  with the results obtained by Ghia et al, 53 who were the first ones to do a comprehensive study on the lid-driven cavity flow, but more so with the results from Botella and Peyret, 54 who used a spectral method in their work.

3D validation of the moment-based boundary conditions
3D cases are the most appropriate representatives of the method's accuracy due to all boundaries being included. Because the LBM together with the moment-based boundary conditions can recover the simplest 2D Poiseuille flow velocity profile exactly, unlike Guo's nonequilibrium method, 23 a 3D duct flow is chosen to test the viscosity (or relaxation time) independence of the method's accuracy and the grid convergence of the algorithm. The former study is performed at a fixed grid verifying whether the solution is set only by the fixed Reynolds number in combination with fixed magic parameter, in this case Re = 100 and Λ = 1∕4. The latter study considers two flow driving mechanisms-external force and pressure difference. The analytical formula for the velocity in a developed 3D square duct flow has been adopted from the theory of elasticity when talking about the deflection surface of a membrane. 55 It has the following form: The expression for the velocity includes the information of the geometry of the square duct or the width L, fluid properties or the dynamic viscosity , and applied conditions or the driving force F z . For pressure-driven duct flow the driving force is the pressure gradient −dp∕dz. The infinite sum is used to account for the rectangular shape of the duct. The coordinates in the velocity solution (32) range from −L∕2 to L∕2 so that the origin (0, 0) is placed in the middle of the duct. If the zeros are substituted into (32) then the formula for the maximum velocity can be obtained as Analytical velocity values of the 3D duct flow are recovered beyond the machine precision and are not affecting the grid convergence study. To evaluate the deviation from the exact solution, the L 2 relative error norm is calculated for the velocity field, where u * i is the exact solution at the calculation domain node i. The results shown in Figure 9 are from 3D pressure-and force-driven flow simulations in Stokes regime meaning that the second-order terms in the equilibrium function (2) are omitted. The viscosity independence study shows that at a fixed grid, Re and Λ the accuracy of the solution remains the same suggesting correct parameterization of the method proposed. In theory, for a force-driven periodic duct flow the density should be constant, however, density variations in the cross-section are observed in the present model. This is a consequence of nonlinear truncation errors that are anisotropic for reduced lattices such as D3Q15 and D3Q19. 56 We remind the reader that this spurious behavior in the numerical solution is observed with bounce-back conditions with the same LBM parameterization and collision model 57,58 -it is a result of the lattice, not the proposed boundary implementation. We note that with the D2Q9 lattice the LBM with moment-based conditions solves the 2D NSEs exactly for Poiseuille flow with no spurious behavior on the minimal grid resolution. 37,38 Exploring potential density variations of the LBM with moment-based conditions on a full D3Q27 lattice is a topic for future studies. A brief parametric study of different Λ values, Λ = {1∕4; 3∕16; 1∕6; 1∕12}, showed that Λ = 1∕12 offers best results in terms of suppressing the density variations and hence the spurious currents in the cross-section of a 3D duct flow on a D3Q19 lattice. This is in accordance with the findings of Bauer and Rüde. 57 Despite the presence of the spurious currents, the solution in the flow direction is unaffected. The grid convergence plot shows that the moment-based boundary conditions are at least second-order accurate in the region where the grid size error is dominant. This is true for both the force-driven and pressure-driven flow case. The bounce-back scheme in combination with TRT and Λ = 3∕16 used in modeling the force-driven duct flow also shows the same order accuracy. The convergence results for the 3D duct flow have been compared with those reported by Mei et al 59 and Hecht and Harting. 27 Although both methods show second-order convergence, the present method overall shows better accuracy. In addition to the 3D duct flow case, the grid convergence study is also performed using the 3D lid-driven cavity flow. Reynolds number is fixed at Re = 1000 and the lid velocity is kept constant at u lid = 0.1 in lattice units so that only the viscosity changes proportionally with the grid size. Figure 10 shows the results of the convergence study. The data are gathered from comparing the velocity field values at different grid sizes. The LBM results of the present work are compared with those obtained by Albensoeder and Kuhlmann 60 (AK) who used a spectral method. The results of the present LBM at Re = 100 are also added for comparison. Because the lid driven cavity flow does not have an analytical solution, the variable field values obtained on the finest grid of 257 3 are used as the reference. For this purpose, (34) is expanded to a 3D vector field calculation as where u * i is the velocity field value on the finest grid, here 257 3 . The results show that the present method remains second-order accurate for the 3D lid-driven cavity flow. Although the AK results yield better accuracy than the present model at Re = 1000, their spectral method shows only O(N −1.5 ) grid convergence. The velocity profiles along the center lines are plotted and compared with the results obtained by Albensoeder and Kuhlmann 60 and Mei et al, 59 see Figure 11. In the present case, the velocity field at Re = 100 and Re = 1000 is calculated using a 33 3

CONCLUSIONS
For the first time, the moment-based boundary conditions have been extended to three dimensions and explicitly derived for the D3Q19 model. Velocity, pressure, and pressure-velocity conditions have been specified for the face, edge, and corner boundaries, where appropriate. Due to the large 3D stencil, the use of the higher order moments on the boundaries is inevitable. For some cases, even when using the more sophisticated TRT collision scheme, it has been observed that the more 'physical' or the asymmetric choice of moments on the corner boundary can lead to instabilities and further to the solution divergence at higher Reynolds numbers. However, the problem can be resolved by choosing a symmetric set of velocity moments, as predicted by Mohammed and Reis. 43 Unfortunately, sometimes it means dropping the condition for the hydrodynamic moment in favor of the third-order moment, such as at edges and corners. It is noted that the LSOB method 28 solves the rectangular matrix formulation and seems to avoid any symmetry breaking. Some preliminary stability analysis for the moment-based boundary method in 2D has been conducted, 37,43 however, an in-depth stability analysis of the current 3D method requires further studies and is a topic for future. At present, a shortcoming of the moment-based method for boundary constraints is its limitation to boundaries aligned with grid points. Overcoming this is a subject for ongoing research. The newly derived 3D moment-based boundary conditions have been validated on benchmark cases testing for the relaxation time dependence, grid convergence and solution accuracy at various flow regimes. The results of the proposed 3D method show excellent agreement with the data obtained using highly accurate spectral methods. Moreover, the method shows the expected second-order grid convergence for all the benchmark cases described. Unlike Guo's nonequilibrium method, 23 the moment-based method can recover the simplest 2D Poiseuille flow velocity profile exactly. In addition, contrary to the bounce-back rule, the present 3D hydrodynamic scheme recovers the no-slip boundary condition for velocity exactly on the wall, the position of which is independent of the relaxation time. This is all achieved locally and thus inherits the main computational advantages of the LBM.
This work comprises the derivation of the 3D moment-based method for the D3Q19 model, but it can be applied to other 3D lattices, D3Q15 and D3Q27 that is. A brief insight into the moment selection for the D3Q15 model is given in

SUPPORTING INFORMATION
Additional supporting information may be found online in the Supporting Information section at the end of this article.

APPENDIX A. DERIVATIONS FOR THE BOUNDARY CONDITIONS
Here the missing steps from the derivation process of the conditions for the velocity and pressure face, edge, and corner boundaries are shown.

A.1 Face velocity
For the west face boundary, five selected moments from Table 2 are rearranged and written out to express the unknown distribution functions. By setting the velocities at the west face boundary to U x , U y , and U z and putting all the unknowns on the left-hand side, the system of equations takes the following form: (A1) Solving the system (A1) yields the unknown distribution functions, which can be further simplified to obtain the solution at the west face velocity boundary given in (18).

A.2 Face pressure
For the west face pressure boundary, the system of equations is composed of the five moments from Table 2. Separating the knowns and unknowns yields the following, where the density 0 is being imposed on the west face boundary. Solving (A3) leads to (21).

A.3 Edge velocity
For the south-west edge boundary, nine selected linearly independent moments from Table 3 are rearranged and written out to express the unknown distribution functions. Setting the south-west edge boundary velocities to U x , U y , and U z , the unknown distribution functions can be calculated as + f 2 + f 5 + f 6 + 2(f 4 + f 16 + f 17 ) + 4(f 8 + f 12 + f 13 ), The density in the above equations is given by the formula which is expressed in terms of the known distribution functions and the relevant moments at the south-west edge boundary.

A.4 Edge pressure-velocity
For the south-west pressure-velocity edge boundary, setting the density to 0 and velocities to U x , U y , U y , and solving for the unknown distribution functions gives the following expressions: − (A6)

A.5 Corner boundaries
The first nine rows of the moment combinations in Table 4 are all linearly independent. The row 10, not obvious at all, is a linear combination of the rows 1, 2, 3, 4, 8, and 9. Setting the density to 0 and velocities to U x , U y , and U z , the unknown distribution functions for pressure-velocity at the low-south-west corner boundary are given as To obtain the unknown distribution functions for the multiple pressure inlet at the low-south-west corner all the velocities in the above system of Equation (A16) have to be set to zero to realize the no-slip boundary condition for velocity, which greatly simplifies the equations.

APPENDIX B. MOMENT COMBINATIONS FOR D3Q15
Here, a possible selection of the velocity moments on the boundaries is shown for the D3Q15 model. There are five unknowns at the face boundary, eight at the edge and 10 at the corner boundary. Having fewer unknowns at the edge and corner boundaries theoretically allows for the selection of the hydrodynamic moments only. That is the case for the face (Table B1) and edge (Table B2) boundaries, and the corner pressure-velocity boundaries. However, the corner velocity boundary still requires a third-order moment to be included into the selection because the density is unknown, see Table B3.
In Table B1, only the hydrodynamic moments are used for both pressure and velocity face boundaries. Although the selection is physically sound, it will produce an antisymmetric solution that might lead to instabilities at higher Reynolds numbers. A symmetric solution can be achieved by replacing the Π yy with Q xzz . By including the third-order moment, the system of equations becomes symmetric and is less prone to developing instabilities.