# Downstream-Conditioned Maximum Entropy Method for Exit Boundary Conditions in the Lattice Boltzmann Method

Authors: Javier A. Dottori; Gustavo A. Boroni; Alejandro Clausse
Journal: Mathematical Problems in Engineering
 (2015)

---

## 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.


---

## Body

## 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 H 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).Figure 1
Set of grid velocities of the D2Q9 model.According to the kinetic picture, the state of a cell located at positionx at time t is given by a population function, fi(x,t), which is distributed on the 9 directions ei pointing to the neighboring cells. The evolution of fi(x,t) is given by a collision step (see (1)) and a streaming step (see (2)):(1)fi′x,t=fix,t+1τfiex,t-fix,t+Δtgix,t,(2)fix+eiΔt,t+Δt=fi′x,t,where Δt is the time step, τ is a relaxation parameter that controls the viscosity, gi(x,t) is an external source rate, and fie(x,t) is called the equilibrium population or distribution. The equilibrium distribution is constructed in terms of the local momenta of fi(x,t):(3)ρx,t=∑i=08fix,t,(4)ux,t=1ρx,t∑i=08eifix,t,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 offie(x,t) determines the partial differential equations that will be approached by the fields ρ(x,t) and u(x,t) 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 is(5)fie=wiρ1+1cs2ei·u+12cs4ei·u2-12cs2u·u,where wi is 4/9 for i=0, 1/9 for the Cartesian directions (1, 3, 5, and 7), and 1/36 for the diagonal directions (2, 4, 6, and 8). From now on, all spatial coordinates are expressed in units of cell’s length Δx, all times in units of Δt, and all velocities in units of Δx/Δt.The parametercs2 in (5) is the factor relating the pressure with the density; namely, p=cs2ρ, which according to the common use is taken as cs2=1/3 in units of Δx/Δt2. The kinematic viscosity in units of Δx2/Δt is related to the relaxation parameter by [17](6)υ=162τ-1.The external sourcegαx,t is used here to introduce external forces Fx,t acting on each cell. The force model [19] used in this study is(7)gi=3wiF·ei.It is known [14, 22] that the equilibrium populations given by (5) maximize the entropy(8)H=∑i=08filog⁡fiwi,subjected to the constraints given by (3) and (4) up to order u2.


## 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 Figure2, there are three unknown populations. In order to perform the streaming step (see (2)), the values of the “upstream” populations, f4, f5, and f6 (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 f4, f5, and f6 is given by the flow boundary conditions at the channel exit.Figure 2
Diagram showing the unknown populations (dashed arrows) in a border cell.


### 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,(9)∂ux∂x=0,(10)∂uy∂x=0.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:(11)uxNx,y=uxNx-1,y.The problem then is a special case of determining the populations pointing in one direction (in this case upstream),f4, f5, and f6, given all the other populations of the cell and the two components of the macroscopic velocity u. The latter can be written in our case (Figure 2) as(12)f4+f5+f6=f2+f1+f8-uxρ=fx,f4-f6=uyρ+f8-f2=fy,where the unknowns f4, f5, and f6 are isolated in the left hand side.It should be noted that the densityρ in the cell is completely determined by ux and the known populations. The summation over the unknowns in (3) and (4) at x-dimension can be eliminated; namely,(13)ρ=11+uxfo+f3+f7+2f1+2f2+2f8.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, f5; that is,(14)f5-f5e=f1-f1e.However, in steady state, mass conservation requires(15)∂ρux∂x=0,which cannot be satisfied together with (9) since ρ is proportional to the pressure in BGK and the pressure gradient along the channel precludes ∂ρ/∂x=0. In order to account for this mismatch, Tong et al. [25] proposed the following correction to (11):(16)uxNx,y=σuxNx-1,y,where(17)σ=∑yρ1,yux1,y∑yρNx-1,yuxNx-1,yis a renormalization factor of the velocity that imposes the inlet mass flow x=1 at the exit x=Nx. This renormalization factor leads to mass conservation and improves convergence to stationary state. Further, the method proposed to use the boundary values of u and ρ to calculate and force the equilibrium populations on all the particles of the cell, even the known populations. Also uy 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 as(18)fiNx,y=1-αfiNx-2,y+αfiNx-1,y,where i=4,5,6. There are known methods that set α=1 [26] and α=-1 [27], respectively. The problem when α=1 is that pressure drops cannot be accurately represented at the end of channel. Using α=2 is a much better selection but can lead to negative values in fi 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 usef4, f5, and f6 values that satisfy (12) and maximize (8). Accordingly, the problem is equivalent to determining f4,f5,f6 such that the Lagrangian(19)L=∑i=4,5,6filn⁡fiwi+λxf4+f5+f6-fx+λyf4-f6-fyis an extreme. The magnitudes λx and λy are the Lagrange multipliers corresponding to the preservation of ux and uy, respectively.Differentiating with respect to each unknown gives(20)∂L∂f4=ln⁡f4w4+1+λx+λy=0,∂L∂f5=ln⁡f5w5+1+λx=0,∂L∂f6=ln⁡f6w6+1+λx-λy=0,which leads to(21)w52f4f6=w4w6f52.Combining (12) and (21),(22)34f52-2fxf5+fx2-fy2=0,where the standard values w4=w6=1/36 and w5=1/9 were assumed. The solution of (22) gives(23)f5=232fx±fx2+3fy2,and the corresponding values of f4 and f6 are(24)f4=12fx+fy-f5,f6=12fx-fy-f5.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 f4 and f6 with the value of f5;(25)f4+12f5=α,(26)f6+12f5=β,where α=(fx+fy)/2 and β=(fx-fy)/2. It can be remarked that α>0 and β>0 are necessary and sufficient conditions for the existence of solution. The expression of f5 is(27)f5=43β+α±β+α2-3αβ.Using (27) in (25),(28)f4=13α-2β∓2β+α2-3αβ=13α-2β∓α-2β2+3α2.Analogously, using (27) in (26),(29)f6=13β-2a∓β-2a2+3β2.So the minus (-) sign in (27) must be used to ensure positive defined values for both f4 and f6. Finally,(30)f5=43β+a-β+a2-3aβ,f4=13a-2β+a-2β2+3a2,f6=13β-2a+β-2a2+3β2are explicit definition for positive defined values that maximize entropy subject to macroscopic constraints.


## 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,(9)∂ux∂x=0,(10)∂uy∂x=0.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:(11)uxNx,y=uxNx-1,y.The problem then is a special case of determining the populations pointing in one direction (in this case upstream),f4, f5, and f6, given all the other populations of the cell and the two components of the macroscopic velocity u. The latter can be written in our case (Figure 2) as(12)f4+f5+f6=f2+f1+f8-uxρ=fx,f4-f6=uyρ+f8-f2=fy,where the unknowns f4, f5, and f6 are isolated in the left hand side.It should be noted that the densityρ in the cell is completely determined by ux and the known populations. The summation over the unknowns in (3) and (4) at x-dimension can be eliminated; namely,(13)ρ=11+uxfo+f3+f7+2f1+2f2+2f8.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, f5; that is,(14)f5-f5e=f1-f1e.However, in steady state, mass conservation requires(15)∂ρux∂x=0,which cannot be satisfied together with (9) since ρ is proportional to the pressure in BGK and the pressure gradient along the channel precludes ∂ρ/∂x=0. In order to account for this mismatch, Tong et al. [25] proposed the following correction to (11):(16)uxNx,y=σuxNx-1,y,where(17)σ=∑yρ1,yux1,y∑yρNx-1,yuxNx-1,yis a renormalization factor of the velocity that imposes the inlet mass flow x=1 at the exit x=Nx. This renormalization factor leads to mass conservation and improves convergence to stationary state. Further, the method proposed to use the boundary values of u and ρ to calculate and force the equilibrium populations on all the particles of the cell, even the known populations. Also uy 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 as(18)fiNx,y=1-αfiNx-2,y+αfiNx-1,y,where i=4,5,6. There are known methods that set α=1 [26] and α=-1 [27], respectively. The problem when α=1 is that pressure drops cannot be accurately represented at the end of channel. Using α=2 is a much better selection but can lead to negative values in fi 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 usef4, f5, and f6 values that satisfy (12) and maximize (8). Accordingly, the problem is equivalent to determining f4,f5,f6 such that the Lagrangian(19)L=∑i=4,5,6filn⁡fiwi+λxf4+f5+f6-fx+λyf4-f6-fyis an extreme. The magnitudes λx and λy are the Lagrange multipliers corresponding to the preservation of ux and uy, respectively.Differentiating with respect to each unknown gives(20)∂L∂f4=ln⁡f4w4+1+λx+λy=0,∂L∂f5=ln⁡f5w5+1+λx=0,∂L∂f6=ln⁡f6w6+1+λx-λy=0,which leads to(21)w52f4f6=w4w6f52.Combining (12) and (21),(22)34f52-2fxf5+fx2-fy2=0,where the standard values w4=w6=1/36 and w5=1/9 were assumed. The solution of (22) gives(23)f5=232fx±fx2+3fy2,and the corresponding values of f4 and f6 are(24)f4=12fx+fy-f5,f6=12fx-fy-f5.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 f4 and f6 with the value of f5;(25)f4+12f5=α,(26)f6+12f5=β,where α=(fx+fy)/2 and β=(fx-fy)/2. It can be remarked that α>0 and β>0 are necessary and sufficient conditions for the existence of solution. The expression of f5 is(27)f5=43β+α±β+α2-3αβ.Using (27) in (25),(28)f4=13α-2β∓2β+α2-3αβ=13α-2β∓α-2β2+3α2.Analogously, using (27) in (26),(29)f6=13β-2a∓β-2a2+3β2.So the minus (-) sign in (27) must be used to ensure positive defined values for both f4 and f6. Finally,(30)f5=43β+a-β+a2-3aβ,f4=13a-2β+a-2β2+3a2,f6=13β-2a+β-2a2+3β2are 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 inuin=0.01 (in units of grid velocity) in the x-direction and the relaxation parameter τ=0.6. 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 α=1 [26] and α=-1 [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 centerumax=3/2uin. Each method was applied to determine fully developed boundary conditions at the exit. Once the steady state was achieved with a relative tolerance of 10-7, 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.Figure 3
Poiseuille flow. Convergence of the different methods towards the analytical solution as the grid resolution is increased. The relative error of the exit-velocity profile against the analytic solution is plotted against the number of cells covering the channel width.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. Figure4 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 Nx<30, and for longer channels it approaches the accuracy of M1.Figure 4
Poiseuille flow. Contour map of the local deviation of the velocity when compared with the numerical solution in a long channel (Nx=500).Figure 5
Poiseuille flow. Dependence of the maximum error on the channel length.


### 4.2. Channel with a Transversal Volumetric Force

The second case studied is the flow in a rectangular channel with a uniform volumetric forceFy=-ρg 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).Figure 6
Flow in a channel with transversal volumetric force. Contour map of the local deviation of the velocity when compared with the numerical solution in a long channel (Nx=500).Figure 7
Channel with volumetric force. Dependence of the maximum error on the channel length.An analytical reference can be obtained in such a case assuming fully developed steady-state conditions, which leads to(31)-∂p∂x+μd2udy2=0,-∂p∂y+ρg=0.Using the equation of state p=cs2ρ, (31) give(32)∂ρ∂x=μcs2d2udy2,(33)∂ρ∂y=-gcs2ρ.Integrating (33),(34)ρx,y=ρoxexp⁡-gycs2.Differentiating (34) with respect to x,(35)∂ρ∂x=dρodxexp⁡-gycs2.Combining (32) and (35) yields(36)μcs2d2udy2=dρodxexp⁡-gycs2.For small gH/cs2 results(37)u=-12μdpodxyo2-y21-gy3cs2,where the channel center line was set at y=0 and the lateral walls were set at y=±yo.The Poiseuille flow is recovered forg=0, as expected. It will be easier to compare the deviation of the profile with respect to Poiseuille case; that is,(38)Δu=ug-ug=0=uingyo32cs21-y2yo2yyo,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. Figure8 shows the map of Δu/uin(gyo3/2cs2) representing the deviation from the velocity development map in a channel without force. Equation (38) is the analytical reference for the profile of Δu 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.Figure 8
Channel with transversal volumetric force. Map ofΔu/uin(gyo3/2cs2) representing the deviation from the velocity without force.


### 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, f4∗, f5∗, and f6∗, are calculated in each cell using (23) and (24), and then the maximum achievable entropy is evaluated in each cell with (8) but using f4∗, f5∗, and f6∗ 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.Figures9 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.Figure 9
Poiseuille flow. Distance between the actual specific entropy and the maximum achievable specific entropy.Figure 10
Channel with transversal volumetric force. Distance between the actual specific entropy and the maximum achievable specific entropy.Figure 11
Map of the module of the velocity gradient. Poiseuille flow (above), channel with transversal volumetric force (down).


### 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 51×340 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.Figure 12
Map of the velocity module in a rectangular channel with a sudden expansion.


## 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 centerumax=3/2uin. Each method was applied to determine fully developed boundary conditions at the exit. Once the steady state was achieved with a relative tolerance of 10-7, 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.Figure 3
Poiseuille flow. Convergence of the different methods towards the analytical solution as the grid resolution is increased. The relative error of the exit-velocity profile against the analytic solution is plotted against the number of cells covering the channel width.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. Figure4 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 Nx<30, and for longer channels it approaches the accuracy of M1.Figure 4
Poiseuille flow. Contour map of the local deviation of the velocity when compared with the numerical solution in a long channel (Nx=500).Figure 5
Poiseuille flow. Dependence of the maximum error on the channel length.


## 4.2. Channel with a Transversal Volumetric Force

The second case studied is the flow in a rectangular channel with a uniform volumetric forceFy=-ρg 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).Figure 6
Flow in a channel with transversal volumetric force. Contour map of the local deviation of the velocity when compared with the numerical solution in a long channel (Nx=500).Figure 7
Channel with volumetric force. Dependence of the maximum error on the channel length.An analytical reference can be obtained in such a case assuming fully developed steady-state conditions, which leads to(31)-∂p∂x+μd2udy2=0,-∂p∂y+ρg=0.Using the equation of state p=cs2ρ, (31) give(32)∂ρ∂x=μcs2d2udy2,(33)∂ρ∂y=-gcs2ρ.Integrating (33),(34)ρx,y=ρoxexp⁡-gycs2.Differentiating (34) with respect to x,(35)∂ρ∂x=dρodxexp⁡-gycs2.Combining (32) and (35) yields(36)μcs2d2udy2=dρodxexp⁡-gycs2.For small gH/cs2 results(37)u=-12μdpodxyo2-y21-gy3cs2,where the channel center line was set at y=0 and the lateral walls were set at y=±yo.The Poiseuille flow is recovered forg=0, as expected. It will be easier to compare the deviation of the profile with respect to Poiseuille case; that is,(38)Δu=ug-ug=0=uingyo32cs21-y2yo2yyo,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. Figure8 shows the map of Δu/uin(gyo3/2cs2) representing the deviation from the velocity development map in a channel without force. Equation (38) is the analytical reference for the profile of Δu 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.Figure 8
Channel with transversal volumetric force. Map ofΔu/uin(gyo3/2cs2) representing the deviation from the velocity without force.


## 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, f4∗, f5∗, and f6∗, are calculated in each cell using (23) and (24), and then the maximum achievable entropy is evaluated in each cell with (8) but using f4∗, f5∗, and f6∗ 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.Figures9 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.Figure 9
Poiseuille flow. Distance between the actual specific entropy and the maximum achievable specific entropy.Figure 10
Channel with transversal volumetric force. Distance between the actual specific entropy and the maximum achievable specific entropy.Figure 11
Map of the module of the velocity gradient. Poiseuille flow (above), channel with transversal volumetric force (down).


## 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 51×340 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.Figure 12
Map of the velocity module in a rectangular channel with a sudden expansion.


## 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.



---

*Source: 159418-2015-07-30.xml*