Abstract
A method for modeling outflow boundary conditions in the lattice Boltzmann method (LBM) based on the maximization of the local entropy is presented. The maximization procedure is constrained by macroscopic values and downstream components. The method is applied to fully developed boundary conditions of the Navier-Stokes equations in rectangular channels. Comparisons are made with other alternative methods. In addition, the new downstream-conditioned entropy is studied and it was found that there is a correlation with the velocity gradient during the flow development.
1. Introduction
The lattice Boltzmann method (LBM) is a numerical tool that has demonstrated great potential in numerous applications, mainly in fluid dynamics simulations. The most interesting advantages are the capability of direct parallelization and the flexibility to include any type of forces and internal interactions. LBM is particularly advantageous in simulation of flow in porous media, multiphase and multicomponent flows, and hemodynamics, among other domains [1, 2].
Since LBM is based on a kinetic representation of the fluid, the task of implementing boundary conditions consistent with the macroscopic constraints is not trivial and requires special attention. On the one hand, the kinetic picture provides certain degree of flexibility that can be used to design boundary conditions, but unfortunately the precision and the stability of the scheme are influenced by the criteria used to fill the available degrees of freedom. This leads to a large amount of work done focused on the management of boundary conditions. Since the early 1990s, many papers have proposed and investigated the behavior of various boundary conditions [3–5].
The most common boundaries that fluid simulations have to deal with are solid walls [6, 7], constant velocity, and flow input or exits. In general, it is sufficient to specify either velocity or density (not both of them). A ubiquitous approach for solid boundaries management is bounce-back scheme. For density and velocity boundaries, there are several specific methods developed for straight walls [8]. For curved boundaries, unfortunately, most of the existing methods are based on interpolation schemes that violate some local conservation properties, which is one of the most attractive features of LBM [9–12].
On the other hand, several researchers have pointed out that the kinetic nature of LBM calls for the direct application of entropic principles in order to improve the stability of the scheme [13–16]. These studies were mainly focused on the search for more stable algorithms while keeping the attractive numerical properties of the method, namely, explicit scheme and locality. A remarkable advance in this direction is the development of the so-called entropic LBM (ELBM), based on the theorem.
In the present work, a method for modeling outflow boundaries in LBM is presented. The method is based on the maximization of the local entropy to determine the upstream components, provided that the downstream components are known. An explicit algebraic formula is derived, which is very easy to implement without perturbing the performance of the calculation. The method is applied to the simulation of fully developed boundary conditions and compared with other alternatives previously presented in the literature. Moreover, the new downstream-conditioned entropy is shown to have an interesting correlation with the velocity gradient.
2. The Lattice Boltzmann Method
The LBM can be viewed as a class of kinematic transport model with an internal discrete vector variable, whose local averages approach transport equations in the infinitesimal limit [17]. The present study is restricted to the LBM representation of the Navier-Stokes equations in 2D. The method operates on a regular lattice and represents a viscous fluid by means of particle populations which moves between the lattice cells. The most common lattice geometry for simulating Navier-Stokes equations in 2D is the so-called D2Q9 model [18], which corresponds to a square grid with 9 internal velocity vectors (Figure 1).
According to the kinetic picture, the state of a cell located at position at time is given by a population function, , which is distributed on the 9 directions pointing to the neighboring cells. The evolution of is given by a collision step (see (1)) and a streaming step (see (2)):where is the time step, is a relaxation parameter that controls the viscosity, is an external source rate, and is called the equilibrium population or distribution. The equilibrium distribution is constructed in terms of the local momenta of :which represent the fluid density and the velocity, respectively. Sometimes a force correction term is included in (4) for accuracy; however, Mohamad and Kuzmin [19] have shown that the correction term is not necessary and better results are obtained without it.
The functional form of determines the partial differential equations that will be approached by the fields and in the differential limit. In the present work, the BGK scheme will be used, which produces the Navier-Stokes equations [17, 20, 21]. In 2D, the BGK equilibrium function iswhere is for , for the Cartesian directions (1, 3, 5, and 7), and for the diagonal directions (2, 4, 6, and 8). From now on, all spatial coordinates are expressed in units of cell’s length , all times in units of , and all velocities in units of .
The parameter in (5) is the factor relating the pressure with the density; namely, , which according to the common use is taken as in units of . The kinematic viscosity in units of is related to the relaxation parameter by [17]
The external source is used here to introduce external forces acting on each cell. The force model [19] used in this study isIt is known [14, 22] that the equilibrium populations given by (5) maximize the entropysubjected to the constraints given by (3) and (4) up to order .
3. New Method to Model Outlet Boundaries
In the case of an outlet boundary cell at the exit of a channel like the one shown in Figure 2, there are three unknown populations. In order to perform the streaming step (see (2)), the values of the “upstream” populations, , , and (i.e., the stream of particles coming from “outside”), are unknown, due to the fact that they do not belong to the fluid domain. The information available to impose values on , , and is given by the flow boundary conditions at the channel exit.
3.1. Other Methods
The classic methods to model outlet boundaries in LBM can be classified into two groups. On the one hand, there are methods that work at macroscopic level, using a velocity boundary condition and computing the velocity from local neighborhood microscopic values. On the other hand, there are methods that work at microscopic level interpolating population values from neighborhood without taking into account macroscopic variables.
Let us assume that the flow at the exit is fully developed; that is,The first condition is enforced using some numerical estimation of the partial derivative in terms of the upstream values. The simplest way would be a 1-cell interpolation:
The problem then is a special case of determining the populations pointing in one direction (in this case upstream), , , and , given all the other populations of the cell and the two components of the macroscopic velocity . The latter can be written in our case (Figure 2) aswhere the unknowns , , and are isolated in the left hand side.
It should be noted that the density in the cell is completely determined by and the known populations. The summation over the unknowns in (3) and (4) at -dimension can be eliminated; namely,
Since there are three unknowns, an additional relation should be provided in order to fully determine the cell state. This is the main difficulty while implementing boundary conditions in LBM [23]. This choice is not trivial since it can affect the precision and the convergence over the whole domain, as pointed out in [8].
One of the most used assumptions to fill this gap is the proposal by Zou and He [24], which consists in assuming bounce-back conditions in the nonequilibrium population traveling perpendicular to the exit, ; that is,
However, in steady state, mass conservation requireswhich cannot be satisfied together with (9) since is proportional to the pressure in BGK and the pressure gradient along the channel precludes . In order to account for this mismatch, Tong et al. [25] proposed the following correction to (11):whereis a renormalization factor of the velocity that imposes the inlet mass flow at the exit . This renormalization factor leads to mass conservation and improves convergence to stationary state. Further, the method proposed to use the boundary values of and to calculate and force the equilibrium populations on all the particles of the cell, even the known populations. Also is assumed to be zero. This leads to introducing errors in simulation according to our results.
The population level methods can be in general described aswhere . There are known methods that set [26] and [27], respectively. The problem when is that pressure drops cannot be accurately represented at the end of channel. Using is a much better selection but can lead to negative values in when the outflow is not smooth.
3.2. Fully Developed Maximum Entropy Method
The natural criterion that it is proposed here is to impose maximum entropy as the additional condition to determine unknown populations. Thus, we propose to use , , and values that satisfy (12) and maximize (8). Accordingly, the problem is equivalent to determining such that the Lagrangianis an extreme. The magnitudes and are the Lagrange multipliers corresponding to the preservation of and , respectively.
Differentiating with respect to each unknown giveswhich leads toCombining (12) and (21),where the standard values and were assumed. The solution of (22) givesand the corresponding values of and are
The sign in (23) should be chosen to ensure that all the populations are positive.
Equations (12) were solved to get an explicit relation of and with the value of ;where and . It can be remarked that and are necessary and sufficient conditions for the existence of solution. The expression of isUsing (27) in (25),Analogously, using (27) in (26),So the minus sign in (27) must be used to ensure positive defined values for both and . Finally,are explicit definition for positive defined values that maximize entropy subject to macroscopic constraints.
4. Results
The performance of the method presented here was tested in a case of viscous flow development at steady state in rectangular channels with uniform velocity at the inlet. The inlet velocity was set in (in units of grid velocity) in the -direction and the relaxation parameter . Null velocity was imposed at the lateral walls using standard midway bounce-back boundary conditions [3].
In order to compare with other methods for boundary conditions, each test case was solved with different alternatives. The first is the present maximum entropy method. The second (Zou and He) and third (Tong et al.) methods were described in previous section. In the latter and in the maximum entropy method, the compressibility factor given by (17) was used to assess the velocity at the exit. The fourth and fifth methods use (18), setting [26] and [27], respectively. For shortness, the methods are numbered M1 (present method), M2 (Zou and He), M3 (Tong et al.), M4 (1-point interpolation), and M5 (Yu’s 2-cell interpolation).
4.1. Poiseuille Flow
The first test studied is the Poiseuille flow where no volumetric forces are applied in the fluid. In such a case, the analytical solution of the fully developed velocity profile is a parabola with maximum velocity at the channel center . Each method was applied to determine fully developed boundary conditions at the exit. Once the steady state was achieved with a relative tolerance of , the exit profile was compared with the analytical solution. The maximum error at the outlet with respect to the inlet velocity was used as a metric. The channel width was varied maintaining the aspect ratio. Figure 3 shows the convergence of the different methods towards the analytical solution as the grid resolution is increased. The maximum entropy method (M1), Zou and He’s method (M2), and Yu’s 2-cell interpolation (M5) show the best convergence rate (order ~ 2). On the opposite side, the 1-cell interpolation method (M4) degrades the convergence to order ~ 1 and Tong et al.’s proposal of using the equilibrium populations (M3) does not converge.
Since the exit boundary conditions that are imposed at the exit are consistent with fully developed flow, the channel should be longer than the development length of the Poiseuille flow. In order to assess the influence of the exit boundary conditions in the flow development, the flow in a 500-cell long channel was calculated and the solution up to the first 200 cells, where fully developed conditions are already reached, was taken as reference. Then a 200-cells long channel was calculated with each method and the solution was compared cell to cell with the reference. Figure 4 shows the map of Euclidean distances between the local velocity vectors normalized with the inlet velocity, obtained with each method. It can be seen that the lowest errors are achieved with the maximum entropy method (M1), followed by the 2-cell interpolation method (M5). The other methods are more than two orders of magnitude worse than M1, especially at the channel exit. On the other hand, the maximum error in the map for each method is a function of the channel length; that is, the longer the channel, the lower the error. Figure 5 shows the dependence of the maximum error on the channel length. As can be seen, the maximum entropy boundary achieves the lowest errors, ~1 order and ~2 orders of magnitude lower than M2 and M3, respectively. The error of the 1-cell interpolation method (M4) is the highest. The 2-cell interpolation is unstable for , and for longer channels it approaches the accuracy of M1.
4.2. Channel with a Transversal Volumetric Force
The second case studied is the flow in a rectangular channel with a uniform volumetric force applied in each cell. Analogously, as in the Poiseuille flow, the influence of the exit boundary conditions in the flow development was assessed by comparing the solutions in a 200-cell channel with the first 200 cells of a 500-cell channel. Figure 6 shows the map of Euclidean distances between the local velocity vectors normalized with the inlet velocity, obtained with each method (analogous to Figure 4). The results are similar to those obtained in the Poiseuille case, with the lowest errors being achieved with the maximum entropy method (M1). Figure 7 shows the dependence of the maximum error on the length of the channel, and the results are similar to the Poiseuille case (Figure 5).
An analytical reference can be obtained in such a case assuming fully developed steady-state conditions, which leads toUsing the equation of state , (31) giveIntegrating (33),Differentiating (34) with respect to ,Combining (32) and (35) yieldsFor small resultswhere the channel center line was set at and the lateral walls were set at .
The Poiseuille flow is recovered for , as expected. It will be easier to compare the deviation of the profile with respect to Poiseuille case; that is,where the pressure gradient was related with the inlet velocity by means of the Poiseuille profile.
Again, the maximum entropy and the other methods were applied to impose fully developed boundary conditions at the exit. Figure 8 shows the map of representing the deviation from the velocity development map in a channel without force. Equation (38) is the analytical reference for the profile of in the fully developed region. It can be seen that the application of methods M1 (maximum entropy) and M5 (2-cell interpolation) agrees very well with the reference solution in the developed region. In turn, the other methods lead to erroneous solutions.
4.3. Entropies and Velocity Gradients
An interesting result can be obtained by mapping the distance between the actual entropy and the maximum achievable entropy that would be obtained if the method presented here were applied in each cell regardless of it being or not at the exit. That is, once the steady state is reached, the entropy in each cell given by (8) is calculated after a streaming step is completed. Afterwards, the conditioned-maximum entropy upstream populations, , , and , are calculated in each cell using (23) and (24), and then the maximum achievable entropy is evaluated in each cell with (8) but using , , and instead of the actual upstream populations. This magnitude can be viewed as an indicator of the departure from the equilibrium structure of the fully developed flow.
Figures 9 and 10 show the maps of entropy distance to equilibrium for Poiseuille case and the channel with volumetric force. Entropies are divided by the density to discount the pressure gradient; that is, it is the entropy difference per particle. It can be seen that in both cases with the maximum entropy method at the exit boundary the entropy distance decays rapidly to zero after the flow development is completed all the way to the exit. Methods M2, M4, and M5 also get similar trends. However, M3 produces significant disturbances at the exit. Figure 11 shows the maps of the modulus of the velocity gradient in both cases. It is easy to note a correlation between the entropy distance and the velocity gradient, which suggests that the latter is the source of the former.
4.4. Backward-Facing Step Flow
Finally, a case of a rectangular channel with two different cross sections was simulated, a classical geometry often called backward-facing step flow [7]. The inlet cross section is reduced by step sized 2/3 and 3/4 of the channel width and length, respectively. The channel is represented by a grid of cells. Figure 12 shows the maps of velocity module obtained with each method. Method M3 is unstable and does not reach the steady-state solution. Methods M4 and M5 introduce numerical noise at the exit, which prevents them from getting fully developed conditions. The method proposed here, M1, and method M2 perform very well, leading to similar results, with the correct fully developed profile at the outlet.
5. Conclusions
The maximum entropy principle was applied in LBM to complete variables undetermined by streaming in boundary cells with outlet condition. The method was implemented in 2D, but its generalization to 3D is straightforward. The resulting algorithm boils down to explicit algebraic formulas that are easy to implement without compromising the performance of general scheme. The proposed method was tested on simulations of flow development in rectangular channels with and without volumetric forces. The method was compared with other alternative boundary techniques previously presented in the literature, showing good precision and stability in all cases.
Furthermore, the deviation of the actual entropy in each interior cell from the maximum entropy was mapped and shows an interesting correlation with the velocity gradient. Since in many cases numerical instabilities in LBM start in regions where the latter is high, the application of the present method can be further studied for stabilization purposes.
Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.
Acknowledgment
This work was partially supported by PICT 2014-1730 FONCYT-ANPCYT of Argentina.