Advances in Water Resources Wetting boundary condition for the color-gradient lattice Boltzmann method: Validation with analytical and experimental data

In the color gradient lattice Boltzmann model (CG-LBM), a ﬁctitious-density wetting boundary condition has been widely used because of its ease of implementation. However, as we show, this may lead to inaccurate results in some cases. In this paper, a new scheme for the wetting boundary condition is proposed which can handle complicated 3D geometries. The validity of our method for static problems is demonstrated by comparing the simulated results to analytical solutions in 2D and 3D geometries with curved boundaries. Then, capillary rise simulations are performed to study dynamic problems where the three-phase contact line moves. The results are compared to experimental results in the literature (Heshmati and Piri, 2014). If a constant contact angle is assumed, the simulations agree with the analytical solution based on the Lucas–Washburn equation. However, to match the experiments, we need to implement a dynamic contact angle that varies with the ﬂow rate.


Introduction
Understanding multiphase flow in porous media is important for several industrial applications such as hydrocarbon recovery and subsurface storage of carbon dioxide. In these applications, at the pore scale, the flow is generally very slow and dominated by capillary forces arising from the interfacial tension between fluid phases. Under these conditions, the displacement of one phase by another is strongly influenced by the wettability of the rock surface.
To compute multiphase flow properties at the pore scale several direct numerical simulation methods have been developed using the volume-of-fluid method ( Raeini et al., 2015;Shams et al., 2018 ), the level set method ( Jettestuen et al., 2013;Prodanovi ća n d Bryant, 2006 ) and the lattice Boltzmann method (LBM) ( Boek et al., 2017;Leclaire et al., 2017;Ramstad et al., 2012;Tölke et al., 2006 ). Several studies have used the LBM to compute flow through porespace images ( Ahrenholz et al., 2008;Boek et al., 2017;Leclaire et al., 2017;Ramstad et al., 2012;. In these cases, a color gradient (CG) LBM was used. The CG-LBM was originally developed by Gunstensen et al. (1991) and it has been improved since then ( Grunau et al., 1993;Halliday et al., 2007;Liu et al., 2012 ). In the CG-LBM, a fictitious density wetting boundary condition proposed by Latva-Kokko and Rothman (2005b) is most widely used to represent the contact angle. In this boundary condition, the solid phase is considered as a mixture of fluid phases for which a fictitious density is assigned. However, a recent study angle implementation is also described and the results are compared with experimental data ( Heshmati and Piri, 2014 ).

Multiphase lattice Boltzmann method
Our 3D immiscible two-phase lattice Boltzmann method is based on the color-gradient approach proposed by Halliday et al. (2007) . In their work, the interfacial tension between two fluids is modeled based on the continuum surface force (CSF) model of Brackbill et al. (1992) for a D2Q9 lattice model. The only difference from their model is that we employ a D3Q19 lattice model. A detailed description of our two-phase lattice Boltzmann method is given in Appendix A .

Wetting boundary condition 2.2.1. The fictitious-density wetting boundary condition
The fictitious-density boundary condition was originally proposed by Latva-Kokko and Rothman (2005b) . In this method, at solid nodes x s , the color function given by Eq. (A.12) is computed using a fictitious density assigned at solid lattice nodes. The color function at solid nodes N ( x s ) takes effect when the color gradient ∇ N of the flow domain next to the solid node is calculated. Based on the force balance at the contact line, the contact angle is expressed as = arccos ( ) . (1) Since this boundary condition is easy to implement, several studies using the CG-LBM adopted this approach ( Gunde et al., 2013;Huang et al., 2014;Liu et al., 2014;Ramstad et al., 2012;. However, Leclaire et al. (2016) showed that this method can be inaccurate because of unphyiscal mass transfer of the wetting phase along the solid. We will also demonstrate unphysical mass transfer using this approach and quantify the resultant error in contact angle and capillary pressure.

Imposing a contact angle directly
The basic idea of our wetting boundary condition is to modify the direction of the color gradient ∇ N at the boundary according to a specified contact angle. Liu et al. (2015) showed an implementation of this wetting boundary condition in the 2D CG-LBM. Their approach is based on the geometrical formulation proposed by Ding and Spelt (2007) . This method accurately simulates wetting phenomena on a flat surface ( Huang et al., 2014;Li et al., 2016;Liu et al., 2015 ) and can be extended to 3D ( Yu et al., 2017 ), but the implementation for arbitrary surfaces is not obvious. Leclaire et al. (2017) proposed a way to find the proper direction of the color gradient ∇ N based on the recurrence relation for the secant method. However, the recurrence relation would in principle require many iterations, which could increase computational costs.
Our wetting boundary condition is similar to the one proposed by Xu et al. (2017) . They proposed a wetting boundary condition for 2D problems. In their work the direction of the color gradient ∇ N is enforced so as to match the prescribed contact angle at the boundary. The direction is obtained by rotating the normal vector to the boundary by the contact angle using the vector transformation equation in 2D. However, this algorithm is only applicable to 2D problems. Our wetting boundary condition is an extension of their method to 3D.
Following the work by Leclaire et al. (2016) and Xu et al. (2017) , lattice sites are divided into four categories, i.e., • C FB : a list of lattice sites that belong to the fluid domain and are in contact with at least one lattice site in the solid domain.
• C Fl : a list of lattice sites that belong to the fluid domain but are not in contact with any lattice sites in the solid domain.
• C SB : a list of lattice sites that belong to the solid domain and are in contact with at least one lattice site in the fluid domain. Fig. 1. Schematic image of the wetting boundary condition. n s is the unit normal vector of the wall while + and − are the possible two unit vectors which have an angle from n s in a counter clockwise and clockwise direction, respectively. These unit vectors, n ± , are obtained by the linear combination of n s and the estimated unit normal vector to the fluids interface, n * , using Eq. (4) , then n * is replaced with either + or − which has the shorter Euclidean distance to n * . Here, ′ is the angle between n s and n * .
• C Sl : a list of lattice sites that belong to the solid domain but are not in contact with any lattice sites in the fluid domain.
A schematic image of the wetting boundary condition is shown in Fig. 1 . First, to define the direction of the wall boundary, the unit normal vector of the boundary n s is calculated for the lattice sites belonging to C FB . For complex geometries, the method presented in Xu et al. (2017) can be used. Because the calculation of the color gradient ∇ N needs to be performed for all the lattice sites belonging to C FB and C Fl by the gradient operator defined by Eq. (A.15) , the color function N at the nodes belonging to C SB is required. This is estimated by the extrapolation of the color function at neighboring lattice nodes which belong to C FB by the following lattice-weighted average scheme: (2) With the values of N at the nodes belonging to C SB , it is possible to estimate the color gradient ∇ N for the lattice nodes in C Fl . This estimated value is denoted as ∇ N * . In the next step, the direction of ∇ N * is modified to match the prescribed contact angle against the wall boundary while keeping the norm of ∇ N * unchanged. The direction of ∇ N * is given by The unit vector n * is understood as the estimated unit normal vector to the red and blue fluid interface at the wall boundary ( x ∈ C FB ). In the modification step, the direction of the color gradient is modified in accordance with the prescribed contact angle . In the 2D problem, Xu et al. (2017) proposed a vector transformation that rotates n s by the angle . However, in the 3D problem, the normal vector that forms an angle to n s makes a circle around n s . We adopted the method employed in the OpenFOAM finite volume library ( OpenFOAM, 2016;Shams et al., 2018 ). In Fig. 1 , two contact lines (points in 2D) are shown. For each contact line, there are two possible unit vectors making an angle with n s (i.e., the unit vector, + , which is rotated by an angle in a counter clockwise direction from n s and the unit vector, − , which is rotated by an angle in a clockwise direction from n s ). These two unit vectors, n ± , are obtained by the linear combination of n s and n * as: The Euclidean distances between n ± and n * are evaluated, then n * is replaced with either + or − which has the shorter Euclidean distance to n * . For instance, at the left contact line in Fig. 1 , n * is replaced with + , whereas at the right contact line, n * is replaced with − . With this, the modified interface normal vector falls into the plane spanned by the unit vectors n * and n s . Note that up to this point, the contact angle is defined as the angle measured through the red fluid to be consistent with the previous literature (for instance Latva-Kokko and Rothman (2005b) ). However from the next section, we refer to the contact angle as the angle measured through the blue fluid since the blue fluid typically represents the denser water phase.

2D static contact angle test
To investigate the difference in accuracy between the fictitiousdensity boundary condition and our method, a 2D static contact angle test is conducted. The equilibrium distributions of wetting and nonwetting fluid placed in a 2D capillary slit are investigated. The simulations are performed in a 2D domain with a mesh containing 101 × 301 lattice nodes. Solid walls are placed on both upper and lower boundaries and periodic boundary conditions are applied for both left and right boundaries.
As shown in Fig. 2 a, blue fluid is initially placed in the center of the domain and red fluid is placed in the rest of the domain. The prescribed contact angle varies from 30 • to 150 • for both the fictitious-density boundary condition and our boundary condition. The other parameters are fixed as = 0 . 1 , 0 = = = 1 . 0 , = 0 . 7 and = = 1 . 0 . The simulations are conducted for 200,000 time steps to achieve an equilibrium state ( Fig. 2 b). Then the simulated contact angle is evaluated using the following geometrical relation:   where is the contact angle, R is the half width of the capillary slit and H is the height of the non-wetting fluid measured from the contact line (point in 2D) as shown in Fig. 2 c. The comparison between the prescribed contact angle and the resultant simulated contact angle for both wetting boundary conditions is shown in Fig. 3 . To quantify the error between an analytical value, X 0 , and a simulated value, X sim , the relative error, E ( X ), is defined by Our method shows an excellent agreement between the prescribed and simulated contact angle for the entire range of values with a maximum relative error of ( ) = 1 . 5% for = 30 • , whereas the fictitious-density boundary condition shows a discrepancy especially for contact angles close to 0 • and 180 • with a maximum relative error of ( ) = 34 . 1% for = 30 • . Fig. 4 shows a zoomed-in view of the color function around the three-phase contact line (point in 2D) for a 30 • prescribed contact angle. For the fictitious-density boundary condition, the interface moves along the wall boundary. This unphysical mass transfer of wetting phase makes the simulated contact angle inaccurate.
We further compare the simulated capillary pressure, P c , for both wetting boundary conditions. The average pressure for each phase is calculated and then the capillary pressure is obtained from the difference. Because of the inaccuracy in the simulated contact angle, the results obtained using the fictitious-density boundary condition show a discrepancy from the analytically evaluated capillary pressure with a maximum  relative error of ( ) = 47 . 5% for = 120 • , whereas our method shows an excellent agreement with a maximum relative error of ( ) = 1 . 1% for = 30 • ( Fig. 5 ). Fig. 6 shows the fluid velocity distribution around the three-phase contact line with a 60 • prescribed contact angle. In these figures, the fluid/fluid interface is indicated by the yellow line and velocity vectors are color coded based on their magnitude (vectors whose magnitude is below 10 −5 lu/ts (lattice units per time step) are shown in black and vectors above 10 −3 lu/ts are shown in white). In both cases, an unphysical, spurious, current persists around the three-phase contact line. With our method, the spurious current around the three-phase contact line is on the order of 10 −5 to 10 −3 lu/ts and the velocity further away is less than 10 −5 lu/ts. In contrast, with the fictitious-density boundary condition, the magnitude of spurious velocity around the three-phase contact line is greater than 10 −3 lu/ts. Table 1 summarizes the comparison of the static contact angle estimation, capillary pressure estimation and the magnitude of the maximum spurious current for both boundary conditions. It can be concluded that our method gives a more accurate estimation for both contact angle and capillary pressure with lower spurious velocity compared to the fictitious-density boundary condition. As we show later, these features are important when we consider slow flows where capillary forces are significant.

3D static droplet test
To demonstrate the applicability of our wetting boundary condition to a 3D problem, a static droplet test in a 3D domain is performed. The simulation domain is composed of a mesh consisting of 101 × 101 × 101 lattice nodes. Solid walls are placed on the upper and lower boundaries and periodic boundary conditions are applied for all other boundaries. Initially, a semi-spherical shape with a radius of 20 lattice units of the red phase is placed on the lower wall. The prescribed contact angle varies from 30 • to 150 • . Simulations are conducted until they reach equilibrium conditions (500,000 ts).
The shapes of the droplet at equilibrium for simulations for different prescribed contact angles are shown in Fig. 7 . Based on the shape of the droplet, the simulated contact angle is obtained and compared with the prescribed contact angle ( Fig. 8 ). As shown in the figure, an excellent agreement with the maximum relative error of ( ) = 2 . 3% for = 60 • is also observed for this 3D test case. Fig. 9. Schematic image of the shape of a red fluid droplet placed on a spherical solid object (green). R 1 and R 2 are the radius of curvature of the red droplet and the radius of the solid object, respectively. The distance between the center of the red droplet and the solid object, R 3 , is determined so as to satisfy = + , where is the contact angle. n p is the unit normal vector to the fluids interface and n s is the unit normal vector of the wall. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

3D static droplet placed on a curved geometry
In this section we validate our method on a 3D surface with a curved geometry. The shape of a red droplet with a radius of curvature R 1 placed on a spherical solid object with radius R 2 can be obtained analytically, as shown in Fig. 9 . The simulation domain is composed of a mesh consisting of 101 × 101 × 101 lattice nodes. A spherical solid object occupies the following region: where R 2 is the radius of spherical solid object that is fixed at 40 lattice units.
The simulations are conducted for contact angles = 60 • and = 120 • . Initially, a known volume of red fluid is placed around the solid object occupying part of a sphere so that the radius of curvature of the red fluid R 1 in equilibrium is 22.5 lattice units. The rest of the domain is filled with blue fluid. Analytical solutions are obtained by calculating R 3 from , R 1 , and R 2 . The simulations are conducted until they reach an equilibrium state (30,000 time steps). Fig. 10 shows the results of the simulations, in which the analytical solution is shown by the white dotted lines. The equilibrium droplet shape agrees with the analytical solution for both = 60 • and = 120 • .

Capillary rise simulation
Our boundary condition allows us to modify the direction of the color gradient at the wall boundary according to the prescribed contact angle at each time step. Therefore, at every time step, it is possible to assign different contact angles. This is particularly useful when we investigate a dynamic contact angle problem.
As an example, capillary rise simulations are conducted. Capillary rise is a phenomenon in which a dense wetting fluid imbibes up a capillary tube until equilibrium between the capillary force and the gravitational force is achieved. Washburn first modeled this behavior and proposed an analytical solution ( Washburn, 1921 ). He assumed a constant contact angle. However, the contact angle during capillary rise does change. There have been several experiments that have measured the dynamic contact angles as a function of capillary number (for instance, Hamraoui and Nylander, 2002;Hoffman, 1975;Rillaerts and Joos, 1980 ). In this section, we refer to Heshmati and Piri (2014) 's experiment. They observed the change of dynamic contact angle during capillary rise with high speed cameras. They performed three experiments using an air/water system in a glass tube with an internal diameter of 0.75 mm, 1.0 mm and 1.3 mm. We define a Bond number by where Δ is the difference in density, g is the gravitational acceleration, is the interfacial tension between two fluids, and r is the characteristic length. This length, r , is the tube radius in 3D or the slit width in 2D. We also define a capillary number by where w is the dynamic viscosity of water and u is the velocity of the fluid/fluid interface. The fluid properties and Bond numbers corresponding to each diameter of capillary tube are shown in Table 2 . Fig. 11 a shows the rise as a function of time and 11 b shows the measured dynamic contact angle as a function of capillary number for the three experimental cases. Sheng and Zhou (1992) showed that the relationship between dynamic contact angle and capillary number is approximated well by where s and d are the static and dynamic contact angle respectively, Ca is the capillary number, K is a constant, L is characteristic length and l s is the slip length. Through least squares fitting of Eq. (10) to the experimental data, = 5 . 59 • and log ( ∕ ) = 39 . 37 are obtained as shown by the black line in Fig. 11 b.
In the literature, capillary rise simulations using the lattice Boltzmann method have been performed ( Latva-Kokko and Rothman, 2005b;Raiskinmäki et al., 2002;Wolf et al., 2010 ). These simulations considered Bond numbers = 10 −1 − 10 0 based on the definition of Eq. (8) , whereas our simulations are performed using ≈10 −2 to compare with the experimental data directly. Hence, compared to previous work, the capillary force in our simulations is stronger compared to gravitational effects. The other key difference is that the experimentally measured dynamic contact angle is used in our simulations.
Our LBM simulations are conducted for capillary rise in a 2D slit to save computation time and to avoid errors associated with the discretization of the circular cross-section of a 3D capillary tube. Assuming laminar flow of an incompressible Newtonian fluids and constant contact angle, the capillary rise is analytically described by the Lucas-Washburn equation ( Hamraoui and Nylander, 2002 ): where r is a characteristic length for the problem and k is a dimensionless flow conductance ( = 8 for the 3D capillary tube and = 12 for the 2D capillary slit). Using the following non-dimensional quantities, * = , ℎ * = 1 ℎ, Based on this equation, we compare our simulation results with Heshmati and Piri (2014) 's experimental data.
First, capillary rise simulations with the constant contact angle implementation are performed. In the simulations, red and blue fluid represent air and water, respectively. A 2D capillary slit with a width of 21 lattice units and 3000 lattice units in length is used. Densities of both red and blue fluids are set at 0 = = = 1 . 0 , but the body force is only applied to the blue fluid. By adjusting the magnitude of the body force, the three Bond numbers corresponding to the experiments are considered. The relaxation times for the red and blue fluids are set at = 0 . 51 and = 1 . 0 , which gives a viscosity ratio similar to the experimental conditions ( = ∕ = 50 ). Initially, the domain is fully filled with air and water is placed at the inlet with a constant pressure boundary condition ( = = 0 ). The outlet boundary is also controlled by a constant pressure boundary condition ( = = 0 ). Since Heshmati and Piri (2014) used a high-speed camera with a large field of view to detect the position of the meniscus and time-stamp the images, there is experimental uncertainty in the determination of the start time corresponding to the frame speed of the camera. Considering this, a time shift is applied to the experimental data points so that the first measured time coincides with the simulated time at the first measured height. This results in 7 ms of time shift for the experiment with a 0.75 mm capillary tube, 5.7 ms for the 1.0 mm capillary tube and 13 ms for the 1.3 mm capillary tube. Fig. 12 shows the simulation results, experimental results and analytical predictions based on the Lucas-Washburn equation with a constant contact angle = 5 . 59 • . The LBM simulations show a good agreement with the analytical prediction by the Lucas-Washburn equation, but both of them show a faster rise compared to the experimental data, especially in the early time region ( t * < 20).
Next, simulations with the dynamic contact angle implementation are conducted. At every time step, the capillary number is computed based on the average velocity of the moving interface. The dynamic contact angle to be imposed as a wetting boundary condition is obtained from Eq. (10) with = 5 . 59 • and log ( ∕ ) = 39 . 37 . The simulation results are shown in Fig. 13 . Compared to the constant contact angle implementation, the simulation results are shifted downward and they show good agreement with the experiments, albeit with a slight difference in the early time region for the case of = 5 . 7 ×10 −3 . Abu-Al-Saud et al. (2017) showed a comparison to Heshmati and Piri (2014) 's capillary rise experiment conducted with an air/Soltrol 170 system in a glass tube with a tube diameter of 1.3 mm using the level-set method. The dynamic contact angle in their model is implemented based on the Cox-Voinov model ( Cox, 1986;Voinov, 1977 ). Their simulations showed good agreement with the trend of measured rise height with a slight difference in the early time region, which is also observed in our simulation. They suggested that this discrepancy was due to inertial effects. Fig. 14 shows a comparison between the prescribed contact angle and the resultant simulated contact angle during the capillary rise simulations. For both the constant imposed angle and the dynamic angle implementation, the simulated contact angles agree with the prescribed values for a wide range of capillary number ( 10 −6 < < 10 −2 ).

Conclusions
We have tested a new implementation of the wetting boundary condition in the CG-LBM framework. Our wetting boundary condition changes the direction of the color gradient according to the prescribed contact angle and it can be applied to 2D and 3D geometries. Therefore, it is suitable for studying multiphase flow in porous media, which will be the subject of future work.
Using the 2D static contact angle test case, the difference in accuracy between the traditional fictitious-density wetting boundary condition and our wetting boundary condition was investigated. When the widely-applied fictitious-density boundary condition was used, because of the unphysical mass transfer along the wall boundary, the simulated contact angle deviated from the prescribed contact angle for contact angles both greater and smaller than 90 • . This deviation resulted in an inaccurate estimation of capillary pressure. In contrast, when our wetting boundary condition was applied, the simulated contact angle and capillary pressure were accurately simulated. In addition, our boundary condition gave smaller spurious currents around the three-phase contact line than the fictitious-density boundary condition. Moreover, we demonstrated that the model worked for a 3D test case of a static droplet on a flat surface and on a curved surface.
To study the applicability of our method to the dynamic problem where the interface and three-phase contact line move, capillary rise simulations were conducted. The simulations were performed with Bond numbers of order 10 −2 to make a direct comparison with previouslypublished experimental data ( Heshmati and Piri, 2014 ). For the constant contact angle implementation, the simulated results agreed with the analytical prediction using the Lucas-Washburn equation. This suggests that the balance between capillary and gravitational forces is properly modeled in our model. Then we demonstrated an implementation of a dynamic contact angle which changes the prescribed contact angle at every time step dependent on the capillary number. The simulated results showed good agreement with the experimental data. For both the constant and dynamic contact angle implementation, the simulated contact angles corresponded to the prescribed values. This suggests that our method precisely controls contact angle for dynamic problems regardless of the fluid velocity.

Acknowledgment
This work is financially supported by Japan Oil, Gas and Metals National Corporation (JOGMEC).

Appendix A. Multiphase lattice Boltzmann method
For the D3Q19 lattice model, the lattice velocity is given as follows: where c = x / t is the lattice speed with x being the lattice length and t the time step size. For simplicity, we set = = 1 . The geometry of the D319Q lattice model can be found in, for instance, Hecht and Harting (2010) . We consider two immiscible fluids, labeled red and blue. The particle distributions of the fluids at position x and time t are denoted as ( , ) and ( , ) . The total particle distribution is given by = + . The density of fluid k is given by the first moment of the particle distribution functions: The fluid velocity u is obtained from the second moment of the particle distribution: where is the total fluid density given by = + . The total particle distribution undergoes collisions as * ( , ) = ( , ) + Ω ( , ) + , (A.4) where * is the post-collision total particle distribution, Ω i is the collision operator which is responsible for the relaxation of the moments towards local equilibrium, and i is the source term and has the effect of introducing a body force in the fluid. In our model, the Bhatnagar-Gross-Krook (BGK) collision operator ( Qian et al., 1992 ) is used: where is the single relaxation time (SRT) and is the equilibrium distribution function which is obtained by a second-order Taylor expansion of the Maxwell-Boltzmann distribution with respect to the local fluid velocity u : where c s is the speed of sound and i is the weight coefficient which is given by Through the Chapman-Enskog expansion analysis of the total particle distribution f i , it is confirmed that the collision operator recovers the continuity equation and the Navier-Stokes equation where the following relations are obtained: where k is the kinetic viscosity of fluid k and p is the hydrodynamic pressure in the fluid. In the interface region, the relaxation time in Eq. (A.5) is determined by the harmonic average of the viscosity of the fluids: The interfacial force is modeled by introducing the body force F at the interface of two fluids. To trace the fluid interface, a color function N is defined as Based on the continuum surface force (CSF), the body force can be expressed as ( Brackbill et al., 1992 ): (A.13) where ∇ N is the gradient of color function N and it is called the color gradient, is the interfacial tension, and is the curvature of the interface which can be calculated by where n is the unit normal vector defined by = − ∇ |∇ | . In Eqs. (A.13) and (A.14) , the partial derivatives are calculated using an isotropic finite difference scheme: where is any function. The calculated body force F which is spatially varying can be imposed through the source term in Eq. (A.4) by the following equation ( Guo et al., 2002 ): [3( − * ) + 9( ⋅ * ) ] ⋅ , (A.16) where u * is the redefined fluid velocity to account for the influence of the spatially varying body force which is given by * = 1 ( ∑ + 1 2 ) . (A.17) Although the body force generates the interfacial force at the interface, it does not guarantee the immiscibility of the fluids. To ensure phase segregation and maintain the sharpness of the interface, the recoloring scheme proposed by Latva-Kokko and Rothman (2005a) is used: where is the segregation parameter which can take a value in (0,1) and is the angle between the color gradient and the lattice vector, which can be obtained by After the recoloring step, the particle distribution functions stream to the neighboring lattice by We briefly describe the boundary conditions used in the paper, except for the wetting boundary condition which is described in detail in the main text. At the solid boundary, the full-way bounce back scheme is applied. In this scheme, the particle distributions at boundary lattice nodes are bounced back into flow domain instead of performing the collision step. This is written as ̃ ( , ) = ( , ) , k = r or b . (A.21) where ̃ is the distribution function in the opposite direction to e i . As for the inlet and outlet boundaries, a periodic boundary condition and a constant pressure boundary condition respectively are used in this paper. For the periodic boundary condition, the outgoing particle distribution from the outlet is re-entered into the flow domain from the inlet. For the constant pressure boundary condition, the method proposed by Zou and He (1997) is used.