A study of wall boundary conditions in pseudopotential lattice Boltzmann models

Abstract The effect of fluid–solid interactions on the hydrodynamics of non-ideal fluids and wettability of surfaces is investigated. We integrate the interaction forces, simulated by pseudopotentials, into two on-site boundary conditions: standard bounce-back (SBB) and Zou and He (ZH) [12] to determine the distribution functions of the boundary nodes. Three different interaction forces are tested: pseudopotential-based interaction (ψ), modified pseudopotential-based interaction (mψ), and a ZH-based interaction, which is proposed by this study based on the ZH method. Therefore, the schemes are ψ-SBB, mψ-SBB, mψ-ZH, and ZH-ZH. The first criterion is the achievement of macroscopic Poiseuille flow. The second criterion is the achievement of a wide range of contact angles. The main method of simulation is multipseudopotential interaction [30]. It is found that the scheme of ψ-SBB creates a relatively large fluctuation of density across the channel. Whilst, the schemes of mψ-SBB, mψ-ZH, and ZH-ZH generate much less density variation across the channel. Among them, ZH-ZH treatment is superior based on density fluctuation and the error associated with the resolution, relaxation time, and compressibility. We found that all four boundary conditions can form a wide of range of contact angles. The ψ-SBB scheme creates largest density fluctuation inside a drop on wettable surfaces. The schemes of mψ-SBB and mψ-ZH create almost the same density fluctuation which is larger than ZH-ZH. Moreover, mψ interaction generates spurious velocities as high as six times a free drop with SBB and eight times with ZH while spurious velocities in ψ-SBB and ZH-ZH are very close to the free drop. Therefore, ZH-ZH performs best, also, in wettability tests.


Introduction
Lattice Boltzmann Method (LBM) is a recent and viable solver in a broad range of fluid flow simulations from simple single-phase, to highly complex multiphase flows [1][2][3] . Since LBM originated from kinetic theory, it has several advantages over other solvers. First of all, it has linear convection or streaming process in phase space. Moreover, without solving a Poisson equation for pressure, which is necessary for incompressible Naiver-Stokes solvers, it recovers Navier-Stokes equations. Furthermore, opting for a limited number of velocity directions, instead of covering the whole phase space, it can straightforwardly convert microscopic distribution function to macroscopic quantities [4] .
Setting up proper treatments at the boundaries of the simulation domain is necessary to achieve the desired correct results while satisfying the given flow conditions, such as velocity, pres- [10] found analytical solutions for the lattice Boltzmann fluid in the case of the two-dimensional Poiseuille flow and Couette flow which facilitate the study of the bounce-back scheme.
The extrapolation treatments are able to provide the second order of accuracy. The finite difference-based scheme of Chen et al. [4] requests information from two neighboring node layers. The non-equilibrium extrapolation scheme of Zhao-Li et al. [11] needs information from the nearest fluid nodes to determine a new set of distribution functions for the nodes whose pressure or velocity is defined. To keep the great assets of local collision operations which suits parallel processing, an LBM boundary condition better finds the unknown distribution functions of a node merely by use of the information available at the node itself.
On the other hand, the on-site boundary condition proposed by Zou and He [12] imposes pressure or velocity condition on the focused node regardless of neighbors. It finds the unknowns with the aid of conservation of mass and momentum, and bounce-back rule for the non-ideal part of distributions. Moreover, the method keeps the second order of accuracy for simulations. Inamuro et al. [13] suggested the missed distributions are found from the equilibrium distribution function that is calculated by a fictitious density and velocity. Skordos [14] filled all the distribution functions on the boundary node using the stress tensor which is related to strain rate tensor. The strain rate tensor can be evaluated by second order finite difference approximation of velocity over neighboring nodes. Latt et al. [15] proposed replacing all the distribution functions of the boundary node by the use of a stress tensor evaluated by the bounce-back of the non-equilibrium part of distribution functions. Latt et al. [15] compared boundary treatments proposed by references [12][13][14] and found that despite the fact that they all had the second order of accuracy, the ones which determine the other unknown distributions after streaming step through closure relations have better accuracy at low Reynolds number flows. Additionally, the boundary conditions, which fill all the distributions at the boundary node, demonstrate better numerical stability. Therefore, they are suitable for high Reynolds number flows.
It should be noted that bounce-back applicability is not limited by the geometry of the wall. However, the others need to find the vectors tangent and normal to the wall, for which various probable situations should be assessed beforehand. They also require further treatments for corner nodes. Bounce-back has been going through considerable extensions; it is combined with spatial interpolations to handle moving boundaries [16] . Junk and Yang [17] modified bounce-back to increase the accuracy of velocity and pressure fields whilst keeping the procedure completely local. Nash et al. [18] studied this category of boundary condition treatments regarding accuracy and performance in various types of flow. It should be noted for curved boundaries, the standard bounce-back reduces the accuracy of the system. In such cases, Mei et al. [ 19 , 20 ] proposed a treatment based on the work of reference [21] that keeps the second order of accuracy in LBM.
One of the well-known capabilities of the LBM is the straightforward simulation of multiphase and multicomponent systems without the need for marker-function advection and front tracking methods. The two-phase systems are formed as a result of Van der Waals (VW) loop seen in equations of state (EOSs) of non-ideal fluids. The most popular model is the pseudopotential model proposed by Shan and Chen (SC) [22] at the early stage of LBM development. It is consistently developing and widely used in simulations and analyses of bubble or drop characteristics, multiphase flow in porous media, and nucleate boiling [23][24][25][26][27][28][29] . Multipseudopotential interaction (MPI) [30][31][32] is a new scheme which ameliorates SC model in different aspects; it removes the inconsistency of SC model, can be applied to a liquid-vapor system with large density ratios, gives the flexibility of changing the interface width, and raises the system stability. In spite of the popularity, the boundary treatments for pseudopotential models have not been studied thoroughly. The SC model utilizes a lattice force, based on the pseudopotential of each node, to model interactions of fluids, and also fluids with solids at a wall boundary.
The lattice force at the boundaries can affect fluid flow behavior in channels and change the wettability of solid surfaces. Therefore, it is possible to model wetting and non-wetting multiphase flow in channels and porous media. The aim of this paper is to study and compare several popular pseudopotential boundary conditions and find out how they perform in modeling fluid flows in channels and wettability of surfaces. For the latter, a liquid drop, which is in equilibrium with its vapor, is placed on surfaces with different levels of wettability. For the former, to narrow down the study, we focus only on single-phase pseudopotential fluids in the channels. The first example of this idea is when pseudopotential model is used to decrease the compressibility of a fluid and decrease the compressibility error in the lattice Boltzmann equation [33] . The second example is a half-saturated porous geometry. There exist channels where a single-phase fluid is flowing and its interactions with solid should be correctly modeled. The third example is the simulation of two-phase flows in which some solid walls are merely interacting with one phase such as the flow regimes of slug, annular and dispersed in a pipe. During this, one phase flows in the center of the pipe and the other one is interacting with the pipe wall.
In this paper, we implement different solid-liquid force interactions into the standard bounce-back and on-site Zou and He [12] treatments, to build non-slip wall boundary conditions. Then these methods are assessed for a single-phase fluid, the compressed water, which flows in a two-dimensional channel and recovers macroscopic Poiseuille flow. The accuracy of simulations, density variation, and velocity profiles due to these treatments are studied. We also demonstrate that the choices of solid-liquid interaction can significantly change the flow behavior. Then those boundary conditions are investigated for two-phase systems, a drop on wettable surfaces. In this case, the availability of contact angles, fluctuation of density near the wall, and magnitude of spurious velocities are examined.

Methodology
The BGK lattice Boltzmann equation (LBE), as a particular discretization of Boltzmann equation [34] for simulation of a flow field, with a general forcing term is defined as where f i are particle distribution functions or particle populations at position x and time t , i = 0 ...q − 1 are the indices of neighboring nodes, e i are the discrete velocity vectors to neighboring nodes, τ is the nondimensional relaxation time, f eq i is the local equilibrium distribution function obtained from the Chapman-Enskog expansion of Maxwellian to the second order at constant temperature [35,36] , and F i is the forcing term. For a flow field near its equilibrium state without internal or external forces, the momentum flux tensor can be obtained by employing the equilibrium distribution function and considering its isotropy on the lattice where α and β are Cartesian coordinates, θ is lattice temperature, and u is velocity. The density and velocity are respectively found from zeroth order and first order moments of distribution functions The definition of velocity is modified when internal or external forces are added to the system which will be explained in the following discussions. The lattice temperature depends on the underlying lattice structure. For the case of the D2Q9 lattice, the discrete c = x/ t is micro-velocity and x is lattice spacing. This choice of lattice sets the lattice temperature as θ = c 2 / 3 . It can be found from Eq. (2) that such a system is governed by the ideal gas pressure. To model non-ideal fluids and two-phase systems, Shan and Chen [22] introduced a lattice force into the LBE, to describe the interaction potentials (pseudopotentials) between fluid particles, where G is the amplitude of the force, N the number of neighboring nodes that interact with the node of interest, c i the vectors linking to those neighbors, and w the weight function that makes the interaction force dependent on the distance. The distance to the nearest neighbor is c . The SC force can be incorporated into the LBE through the forcing term and equilibrium distribution functions.
The authors introduced the concept of the multipseudopotential interaction (MPI), over the single-pseudopotential interaction (SPI), to describe the hydrodynamics of practical fluids, under the condition that each potential satisfies a thermodynamic requirement [31] F total = F (1) + F (2) where n is the number of forces (or pseudopotentials), G j the amplitude, and ψ j the consistent pseudopotential of the j -th part of the force where C j and λ j are arbitrary constants, ε j = −2( a 1 + 12 t s j G j ) / a 2 in which a 1 = 0 , and a 2 = 3 for nearest-node interactions [30] , s j is an arbitrary constant, w ( | c i | 2 ) = 0 for | c i | 2 > 2 and from fourthorder isotropy of the force on the lattice w (1) = 1 / 3 and w (2) = 1 / 12 . The physics of MPI is to present the inter-particle potential with the contributions from a set of sub-forces at various potentials, which are functions of particle densities. The forcing term is α ) (8) and the actual fluid velocity, v , applied in the equilibrium distribution function, is defined as, The MPI forces parameters which reproduce the SRK EOS in the LB system.
where Einstein summation convention is adopted. Meanwhile, the MPI EOS is found to be which is suitable to implement cubic EOSs by adopting proper parameters and number of pseudopotentials, the details can be found in Khajepor and Chen [30] . In this study, we chose Soave-Redlich-Kwong (SRK) EOS for the MPI system to present the thermodynamic behavior of water. SRK is the well-known two-parameter cubic equation of state [37] and is defined by rearrangement as, and ω is Pitzer's acentric factor. α can be determined from experimental vapor pressures of non-polar substances. It should be noted that the reduced format of the SRK depends on the acentric factor ω, which helps to treat non-polar substances T R are reduced pressure, reduced density, and reduced temperature respectively. The compressibility of fluid is proportional to the factor, Y = b 2 /a . The relevant parameters of pseudopotentials of the MPI scheme for SRK EOS are listed in Table 1 .

Boundary conditions
We begin with a brief discussion on the Standard bounce-back (SBB) scheme which will be tested. We then demonstrate how the fluid-solid interaction forces can be implemented into Zou and He scheme [12] for the pseudopotential model. Finally, the various schemes of fluid-solid interactions are discussed.

Standard bounce-back
The standard bounce-back (SBB) is the simplest but popular LBM scheme to treat non-slip boundary conditions. It reverses and sends back the distribution functions penetrated to the stationary solid surface. For a fluid node in contact with the wall, the unknown distribution functions, f i , coming from the solid surface can be found from where i is the opposite direction to i and t + is a post-collision time but before streaming. If the nearest solid nodes are the solid surface especially in complex boundaries, the method is of firstorder accuracy [7,8] . However, if the surface line drawn at halfway between the solid and fluid nodes, the second order of accuracy is achieved [9] . It should be noted in the case of a moving boundary, SBB can be applied straightforwardly by taking into account solid to fluid momentum transfer [1] .

Zou and He treatment with interaction force
Zou and He (ZH) [12] proposed a boundary condition scheme, which can be used to set density or velocity at a particular node.
Please cite this article as: S. Khajepor   The idea is to find unknown distribution functions with the aid of macroscopic values and known distribution functions. Here, in a D2Q9 lattice, we revisit ZH treatment to incorporate the pseudopotential interaction force. For a solid or boundary node, A, at the bottom of a domain shown in Fig. 1 , the distribution functions f 2 , f 5 , and f 6 have to be defined after every collision-and-stream step.
For 2D pseudopotential LBM, we consider interaction forces acting on node A, namely F x and F y . We can set three neighboring nodes at the bottom as three ghost nodes. Ghost nodes are the nodes which are not a part of the physical domain but form a virtual layer around the boundary of the domain. The densities of wall surface nodes are copied on to the nearest ghost nodes to calculate interaction force through Eq. (6) . The bounce-back rule is valid for the non-equilibrium part of the distribution function perpendicular to the wall [12] , From Eqs. (3) , (4) and (9) , it can be found that For a stationary solid node, velocity in the above equations are set to v x = v y = 0 . It should be noted that since pseudopotential LBM pressure is directly related to density via EOS, the pressure boundary condition can be applied straightforwardly. However, the velocity component tangent to boundary surface, in this example v x , should be defined along with the density, or a relation between velocity components should be given, for example, the angle of the velocity vector is known. In such a case, Eq. (14) can be solved for velocity and Eqs. (13) , (15) , and (16) are treated the same as velocity boundary condition to find f 2 , f 5 , and f 6 . In addition to the intersection of a solid wall and periodic boundary, corner nodes are almost inevitable situations. For a D2Q9 lattice, a corner node leaves five unknown and three known distribution functions. For example, node B placed at the bottom left of the domain in Fig. 2 has f 1 , f 2 , f 5 , f 6 , and f 8 as unknowns and f 3 , f 4 , and f 7 as known values. To calculate the interaction force on this node, the densities of solid surface nodes are considered for nearby ghost nodes. We set v x = v y = 0 and assume the density of the node is known due to the side pressure (density) boundary or extrapolation over nearby nodes. From the bounceback rule of non-equilibrium part of both perpendicular distribution functions, we find f 2 = f 4 , f 1 = f 3 . Therefore, the rest can be found from Eqs. (3) , (4) and (9) ,

Fluid-Solid interactions
In pseudopotential LBM, the interaction between solid and fluid is considered to control the wettability of a solid surface in contact with a two-phase fluid. This idea, in fact, mimics the interactions observed at the molecular level and scales up to macro-scale. In the cases where pseudopotential fluid is in the form of the single phase, such as a compressed liquid in contact with stationary solid boundaries, the wettability is less focused but obtaining correct density and velocity profiles along and across the channel are of importance.
The most well-known pseudopotential fluid-solid interactions can be formulized as where G fs is the amplitude, φ f and φ s are fluid and solid potentials respectively, s is a switch function gives 0 for fluid-fluid interactions and 1 for fluid-solid interactions. If a non-wetting fluid is simulated G fs should be positive otherwise negative. Martys and Chen [38] defined a fluid-solid interaction force setting φ f (x ) = ρ(x ) and φ s = 1 for a single component system. Raiskinmäki et al. [39] and Sukop et al. [40] proposed to replace the fluid density factor with pseudopotential as φ f (x ) = ψ (x ) and φ s = 1 . Kang et al. [41] consider a constant density for solid nodes φ f (x ) = ρ(x ) and φ s = ρ s where ρ s is an imaginary density set for a solid node. The pseudopotential version of this model is introduced by Benzi et al. [42] where φ f (x ) = ψ (x ) and φ s = ψ ( ρ s ) . Li et al. [43] assumed the solid node has a density equal to the fluid node which it is interacting with, φ f (x ) = ψ 2 (x ) and φ s = 1 . In this way, the fluidsolid force has the same order of magnitude as fluid-fluid interactions if G and G fs are in the same order. They called their interaction model modified pseudopotential-based interaction (m ψbased).
Martys and Chen [38] and Kang et al. [41]   Pressure, density, and temperature of SRK EOS for water at T R = 0 . 7 for at a = 0 . 01 and b = 0 . 2 and ω = 0 . 344 . The numbers are in lattice unit. β R is reduced compressibility. for the solid node is calculated by Eq. (14) . As such, we define the fluid-solid interaction, ZH-based interaction, This force similar to modified ψ-based force keeps the order of magnitude as fluid-fluid interactions which is the most important point of defining such an interaction force [40] .
To apply the fluid-solid interactions in the MPI framework, they are considered in the same shape as fluid-fluid interactions. Therefore, Eq. (20) becomes In this study, the focused boundary interactions are the schemes of ψ-based interaction, modified ψ-based interaction, and ZH-based interaction in combinations with either the bounceback or the ZH boundary treatments.

Results and discussion
The boundary treatments discussed in Section 3 are assessed based on their viability in modeling non-ideal fluid flows in 2D channels ( Section 4.1 ) and drop contact angles for two-phase systems ( Section 4.2 ). The results from simulations are discussed in comparison with the analytical incompressible solution. The focused parameters for the flows are compressibility, resolution, density fluctuation, and viscosity. The important parameters for the drop on surfaces are the availability of contact angles, density fluctuation, and the magnitude of spurious velocities.

Non-ideal Poiseuille flow
All simulations are run on a two-dimensional square lattice including nine velocities (D2Q9). The detail of domain sizes for each case will be given individually. The general setting, unless otherwise stated, is listed as follows. The nondimensional relaxation time and lattice spacing parameter are set to unity τ = 1 , c = 1 .
ZH pressure boundaries opted for driving force of walled Poiseuille flow. The MPI forces are calculated by use of Eqs. (6) and (7) whose parameters, for the SRK EOS, can be found in Table 1 . The internal forces are embedded in the LBE using Eq. (8) , which reproduces Navier-Stokes equations to the second order. The thermodynamic values of SRK EOS at the saturation state T R = 0 . 7 are listed in Table 2 . The amplitudes of fluid-solid interaction are kept at fluidfluid interactions G The simulations are run for at least 10 5 steps and velocity of nodes all over the system is monitored to verify that the equilibrium state is reached by where v max is the velocity magnitude of the node which has the maximum value in the computational domain. The equilibrium time depends on viscosity, domain size, and type of boundary conditions.
We assess the results in comparison with the analytical solution of Poiseuille flow in a no-slip channel. It is assumed that the flow is fully developed and uniform in x direction and pressure drop is linear along the channel d p where L y is the height of the channel, L x is the width of the channel, μ is the dynamic viscosity, p is the pressure drop along the channel. If a body force, F b , is exerted on the fluid with the acceleration of a , the term p / L x is replaced by | F b | = | ρa | .

Periodic Poiseuille flow
The first set of simulations is run for Poiseuille flow where no solid boundary condition is applied [44] as an ideal test case that the fluid itself creates a Poiseuille profile. To do so, we charge MPI fluid in a fully periodic rectangular domain and exert two equal but opposite body forces. The first one only acts on the top half and the other one moves the bottom half. The sketch of the simulation domain is shown in Fig. 3 . The body force is added to Eq. (9) , thus, show that the pseudopotential LBM has the second order of accuracy same as conventional LBM. In fact, after inspection of density throughout the domain, we found that no density change is observed in the whole domain in all simulations. It means the pseudopotential forces are neutralizing each other since pseudopotentials are calculated by node density and, therefore, the model acts similar to conventional LBM.

Pseudopotential-based interaction
To simulate Poiseuille flow in between two stationary plates, we first employ ψ-based interaction and set SBB, Eq.   In the case Y = 40 , which repulsion part of SRK EOS is dominant and the fluid acts more compressible, the density of fluid from wall decreases smoothly in a parabolic distribution to a minimum density slightly lower than the average density. By increasing the attraction parameter or decreasing Y, where the fluid is more approaching to incompressible fluid, the density forms a plateau in the center closer to average density but density fluctuations near the wall increases. The velocity profile of these cases are shown in Fig. 7 (d). The best agreement with analytical results (incompressible fluid flow) comes from Y = 0 . 4 , but with the increase of Y the fluid moves faster than incompressible one.
The relatively same results are observed when ρ s = 0 . 9 ρ right which are demonstrated in Fig. 7 (c). The density of liquid increases from wall density to the average density at the center and decrease of Y moves the profile of density from smooth parabolic shape to flat profile at the center with fluctuations close to the wall. The associated velocity profiles are shown in Fig. 7 (f). The best consistency between simulation and the incompressible theory Eq. (24 ) is for 0.4 < Y < 4. The worse deviation is seen at the highly compressible fluid, Y = 40 . According to Fig. 7 (b), the less density gradient is seen across the channel if wall density is set to average density ρ s = ρ a v e .
However, we show in Section 4.1.3 that this only happens here because the sample is taken at the halfway along the channel. The velocity profiles, seen in Fig. 7 (e), stay relatively the same as previous cases.

ZH-based and modified pseudopotential-based interaction
To simulate walls in Poiseuille flow, here, we investigate three cases of fluid-solid treatment: ZH-based fluid-solid interaction coupled with ZH zero velocity (ZH-ZH), modified ψ-based interaction with bounce-back (m ψ-SBB ), and modified ψ-based interaction with ZH zero velocity (m ψ-ZH ). The MPI fluid and inlet (or left), and outlet (or right) boundaries are set the same as Section 4.1.2 . The tests focus on density fluctuations across the channel and errors due to resolution, relaxation time, and parameter Y. The error is defined as a deviation from macroscopic Poiseuille flow, Eq. (24) . Fig. 8 shows the contours of density field over the computational domain for ψ -SBB, m ψ -SBB, m ψ -ZH, and ZH-ZH treatments. The density variation across the channel for ψ-SBB is much larger than the others. It is seen that only around halfway through the channel an almost vertical contour level, 10 0 0(ρ/ ρ a v e − 1) = 0 , is formed. Since the wall density is set at ρ s = ρ a v e , that is the location where the density of fluid coincides with the density of the wall. Therefore, any density difference between fluid and solid causes the density fluctuation. Due to high fluid fluctuation, this treatment is not comparable with the others, and hereafter we concentrate on the other treatments.  The error is depicted versus vertical resolution of the domain in Fig. 10 . As expected, the error decreases with the increase of domain resolution. ZH-ZH treatment is superior at different resolutions. After that, m ψ-SBB placed which shows lesser error than m ψ-ZH. It should be noted the horizontal resolution is accordingly increased to make sure that the dynamics of the simulations stay the same.
The effect of relaxation time is shown in Fig. 11 . In the range of 0.6 ≤ τ ≤ 1.1, the m ψ-ZH method shows more error than the others but the error almost reaches a plateau after τ = 0 . 8 . Similar behavior is seen for ZH-ZH treatment which has the lowest error.
However, m ψ-SBB error increases proportionally with τ .    Boltzmann. Such behavior can be seen in Fig. 12 which demonstrates the error dependency to parameter Y. ZH-ZH treatment is better than the others and m ψ-SBB is placed after m ψ-ZH.

Wettability assessment
For two-phase systems, the selected boundary conditions are assessed regarding achievable contact angles, density behavior near the wall, and spurious velocities. The domain size is considered as 300 × 100 in lattice unit (lu). Top and bottom of the domain are wall boundary conditions. Left and right boundaries are periodic. The nondimensional relaxation time and lattice spacing parameter are set to unity τ = 1 , c = 1 . A drop with the radius of R = 30 lu is placed on the wall while it is wetting the surface slightly, so the drop center is ( x c , y c ) = (Lx/ 2 , R − 5) . As the initialized system is not in the perfect equilibrium state, there is a slight pressure wave forming at the beginning of the simulation which may cause the drop to separate from the surface. To avoid such a problem, a small gravity force, G = −10 −6 lu, is exerted in the system which is removed after a few thousand steps. For the case of m ψ-SBB, m ψ-ZH, and ZH-ZH the wettability of the surface is controlled via where A G is a multiplier and G j is the amplitude of fluid-fluid interaction. For the case of ψ-SBB, G ( j) f s = G j and the wettability is controlled via density of the solid wall as a function of vapor density, ρ s = A ρ × ρ v , where A ρ is an arbitrary multiplier.
The fluid thermodynamic properties are mentioned in Table 3 .
The results from the simulations of a static drop on surfaces made with different boundary conditions can be found in Figs. 13-16 . It is identified that the drop forms a wide range of contact angles from ∼15 °to ∼180 °. When the contact angle approaches 180 °, it becomes difficult to keep the drop near the wall because of hydrophobicity. The density of the drop and the vapor around it are deviated from the free drop state due to the change in solid-liquid and solid-vapor surface tension forces, γ sl and γ sv respectively.
The effect of fluid-solid interaction strength, which is presented by A G for the schemes of m ψ -SBB, m ψ -ZH, and ZH-ZH and A ρ for the scheme of ψ-SBB, on contact angles are summarized and plotted on Figs. 17 and 18 . It shows that both A G and A ρ play the same role because the stronger the fluid-solid interactions, the smaller the contact angles. The sound correlations between contact angles and the fluid-solid interaction strengths can be formulated as the  Pressure, density, and temperature of SRK EOS for water at T R = 0 . 7 , a = 0 . 04 , b = 0 . 2 and ω = 0 . 344 . The numbers are in lattice unit. β R is reduced compressibility.  Table 4 The parameters of the correlations of the contact angle with interaction strengths.
Boundary condition scheme Except for the scheme of ZH-ZH, a parabolic correlation works well for the curve-fitting. The contact angle created by ZH-ZH scheme is larger than those by schemes of m ψ-SBB and m ψ-ZH when A G > 1.1 and it is smaller for A G < 1.1. The contact angle is more sensitive to the interaction strength as the slope of contact angle-A G for ZH-ZH scheme is much steeper than those of the other schemes. Fig. 19 depicts the density change about the drop-wall contact line. Higher than y/ y = 7 , all schemes showing the density of their drops. However, getting closer to the wall at y/ y = 6 , hydrophobic ψ-SBB deviates quickly from the drop density. The more hydrophilic, the less deviation is seen from ψ-SBB. The width of the solid-liquid interface is in the range of 1-5 lattices. The boundary schemes of m ψ-SBB and m ψ-ZH affect density profiles almost identically. Their first deviation is seen at y/ y = 4 in their hydrophobic version. They get closer to drop density by decreasing hydrophobicity. However, the very hydrophilic surface of them can create densities greater than inside of the drop. In these cases, the width of the solid-liquid interface is 1-3 lattices. This scenario is generally true for the ZH-ZH method but with considerable less deviation. The hydrophobic version of ZH-ZH separates at y/ y = 3 . The width of the solid-liquid interface, in a wide range of contact angle, is only 1-2 lattices. It should be noted that a smaller density deviation is generally more desirable as it describes a scale much bigger than molecular levels.     | . For the cases ψ-SBB and ZH-ZH and θ > 50 °, the maximum spurious velocity is almost equal to that of a free drop. For smaller angles θ ≤ 50 °, spurious velocities of ZH-ZH grow linearly in contrast to ψ-SBB which decreases linearly. Both m ψ-SBB and m ψ-ZH show a considerable increase of spurious velocities at lower angles about 6 times and 8 times, respectively, of what is seen in the free drop. Their magnitude decreases with the decrease of hydrophilicity, but they are still much bigger than those from the ψ-SBB and ZH-ZH.
We test the applicability of ZH-ZH in dynamic wetting processes. A drop is placed on a hydrophobic surface to form a contact angle about θ ≈ 180 °. When the system reaches equilibrium, it is considered as the initial state for three different surfaces: In each case, we monitor dynamic spreading of the drop and record contact radius, r , which is the radius of the wetted area. Based on Tanner's law [45][46][47] at the final stage, a drop spreads on a very wettable surface θ = 0 • with the relationship r ∼ t 1/10 in 3D cases and r ∼ t 1/7 in 2D cases. However, during the initial stage of spreading, the capillary force is at its maximum to reshape the drop. Therefore, the drop spreads faster Please cite this article as: S. Khajepor  at the early stage. Experiments and simulation results of references [4 8,4 9] showed that, for a fully wetting surface, relationship r ∼ t 1/2 is dominant during the initial stage and, for partially wetting surfaces, the exponent, 1/2, decreases with the increase of the static contact angle. The results from our simulations on this dynamic wetting are given in Fig. 21 . The data are analyzed with the aid of r/R = C (t/ t ρ ) n where t ρ = ρR 3 /γ is the inertial time scale, R is the initial drop radius, and γ is surface tension, and C and n are free parameters, which can be determined from curve-fitting. In all three cases, we see the linear behavior at the beginning; it has the exponent of n = 0 . 5 for the very wettable surface, θ = 18 . 7 • . The exponent decreases to n = 0 . 46 for θ = 64 . 2 • and further down to n = 0 . 38 θ = 114 . 9 • which is generally in good agreement with experiment observations.

Conclusion
We have chosen the two most popular boundary treatments, SBB and ZH. They are combined with fluid-solid interaction forces originating from SC model. These boundary treatments are compared and assessed in two scenarios: Poiseuille flow for compressed liquids and contact angles achieved by putting a drop on wettable surfaces. The MPI scheme with SRK EOS is utilized to model water. Three force interactions between fluid and solid are considered, the ψ-based, the modified ψ-based, which are found in the literature, and the ZH-based which is proposed and defined by this study. For a hydrodynamic study, the criterion is closeness to single-phase macroscopic Poiseuille flow and, for a hydrostatic study, the availability of drop contact angles is considered.
The periodic Poiseuille flow is studied as an ideal Poiseuille flow which is purely made of LB equation without boundary conditions. It is shown that the pseudopotential LB method is secondorder accurate as density is constant all over the channel and the method coincides with conventional LB without pseudopotential interactions.
ψ-based interaction along with SBB gives high-density gradient and fluctuations near the wall. The choice of wall density plays the key role. The closer it is to neighboring fluid nodes, the less fluctuation is observed which makes it difficult to use in complex geometries. The velocity profile is very close to analytical Poiseuille flow unless the parameter Y is greatly increased which makes the fluid more compressible and consequently strengthens the compressibility error of LB method.
The density fluctuations of ZH-ZH, m ψ-SBB, and m ψ-ZH, in the channel, are far less than and not comparable to those from ψ-SBB.
Therefore, those three treatments are studied together. The performance of them assessed based on density variation across the channel, and error because of domain resolution, relaxation time, and Y parameter. ZH-ZH in all cases is superior. It has one order of magnitude less variation of density across the channel in comparison with others. It shows less error in different resolutions, a low relaxation time dependency and less compressibility error. In all cases, m ψ-ZH takes the third place after m ψ-SBB.
For a drop on wettable surfaces, all boundary treatments successfully form a wide range of contact angles. The profiles of the contact angle versus the interaction strength can be formulated with cubic and parabolic polynomial. ψ-SBB shows largest density deviation near the wall. ZH-ZH keeps the density close to drop density. The schemes of m ψ-SBB and m ψ-ZH show the same density behavior and are placed between ψ-SBB and ZH-ZH. Spurious velocities of m ψ -SBB and m ψ -ZH are considerably higher (about 6 and 8 times at lower contact angles) than a free drop while ψ-SBB and ZH-ZH effects on spurious velocities are negligible. ZH-ZH is tested for a dynamic drop-spreading process. It successfully follows the power law which is seen in the literature.
To summarize, ZH-ZH shows best results in the Poiseuille flows and wettability tests. The schemes of m ψ-SBB and m ψ-ZH act well in single-phase flow, however, if there is two-phase interface near the wall, they tend to intensify spurious velocities. The scheme of ψ-SBB demonstrates a large density fluctuation/deviation near the wall when there are density gradients in a simulation domain.