Choice of boundary condition for lattice-Boltzmann simulation of moderate-Reynolds-number flow in complex domains

Modeling blood flow in larger vessels using lattice-Boltzmann methods comes with a challenging set of constraints: a complex geometry with walls and inlets and outlets at arbitrary orientations with respect to the lattice, intermediate Reynolds (Re) number, and unsteady flow. Simple bounce-back is one of the most commonly used, simplest, and most computationally efficient boundary conditions, but many others have been proposed. We implement three other methods applicable to complex geometries [Guo, Zheng, and Shi, Phys. Fluids 14, 2007 (2002); Bouzidi, Firdaouss, and Lallemand, Phys. Fluids 13, 3452 (2001); Junk and Yang, Phys. Rev. E 72, 066701 (2005)] in our open-source application hemelb. We use these to simulate Poiseuille and Womersley flows in a cylindrical pipe with an arbitrary orientation at physiologically relevant Re number (1–300) and Womersley (4–12) numbers and steady flow in a curved pipe at relevant Dean number (100–200) and compare the accuracy to analytical solutions. We find that both the Bouzidi-Firdaouss-Lallemand (BFL) and Guo-Zheng-Shi (GZS) methods give second-order convergence in space while simple bounce-back degrades to first order. The BFL method appears to perform better than GZS in unsteady flows and is significantly less computationally expensive. The Junk-Yang method shows poor stability at larger Re number and so cannot be recommended here. The choice of collision operator (lattice Bhatnagar-Gross-Krook vs multiple relaxation time) and velocity set (D3Q15 vs D3Q19 vs D3Q27) does not significantly affect the accuracy in the problems studied.

Figure 1

Convergence of velocity error residuals (${\epsilon }_u^2$) as defined in Eq. (18) vs diameter in lattice units. (a) Re number $\mathrm{Re}=1$; (b) $\mathrm{Re}=30$; (c) $\mathrm{Re}=100$; (d) $\mathrm{Re}=300$. Colors (black, red, blue) indicate boundary condition (SBB, BFL, GZS) and shapes (crosses, triangles, squares) indicate the velocity set (D3Q15, D3Q19, D3Q27). The solid lines are guides to the eye showing first-order (dashed) and second-order (solid) convergence. The simulations with $\mathrm{Re}=300$ and $D=12\Delta x$ became unstable due to the high Ma number and are not shown.

Figure 2

Simulated (D3Q15, LBGK, BFL) axial velocity (normalized by $U_{\max }$) vs the radial coordinate. The Womersley number $\alpha =12$ and the Re number $\mathrm{Re}=300$. The data is shown at four moments during the pulsatile cycle: $t=T/4$ (magenta, maximum positive pressure gradient), $T/2$ (blue, zero pressure gradient), $3T/4$ (red, maximum negative pressure gradient), $T$ (black, zero pressure gradient). The analytical solutions are shown with solid lines.

Figure 3

Error residuals ${\epsilon }_u^2$ as defined in Eq. (18) for simulations of pulsatile flow vs Womersley number $\alpha$. Colors (black, red, blue) indicate boundary condition (SBB, BFL, GZS) and shapes the collision operator and velocity set. LBGK is indicated by polygons (square, D3Q15; triangle, D3Q19; diamond, D3Q27) and MRT by crosses ($+$, D3Q15; $\times$, D3Q19).

Figure 4

Parameters and coordinate systems for the torus simulations. The radius of the tube is $a$ and the distance from the center of the tube to the center of the torus is $c$. We define a toroidal coordinate system $(r,\theta ,\varphi )$ where $(r,\theta )$ are polar coordinates in the cross-section plane and $\varphi$ is the angle that the cross-section plane makes with the $x$-$z$ plane.

Figure 5

Streamlines for the secondary flow and color maps of the modulus of the error field in a torus, with Dean number $\kappa \approx 112$ and Re number $\mathrm{Re}\approx 200$. The inside of the pipe is on the left of the figure. Across columns, we vary the boundary condition (SBB, BFL, and GZS from left to right). Down the rows, we vary the combination of velocity set at collision operator (D3Q15+LBGK, D3Q19+LBGK, D3Q15+MRT, D3Q19+MRT, D3Q27+LBGK from top to bottom). The figures in color show the simulated results. The color field indicates the absolute error Eq. (31) where light gray (yellow online) is zero and dark (red online) is $0.15U_{\varphi }(0,0)$. Since the combination of MRT and GZS is unavailable, we use these spaces to show the predicted secondary flow calculated using Eqs. (23, 24, 25, 26, 27) in grayscale.

Figure 6

Streamlines for the secondary flow and color maps of the modulus of the error field in a torus, with Dean number $\kappa \approx 224$ and Re number $\mathrm{Re}\approx 400$. The inside of the pipe is on the left of the figure. Across columns, we vary the boundary condition (SBB, BFL, and GZS from left to right). Down the rows, we vary the combination of velocity set at collision operator (D3Q15+LBGK, D3Q19+LBGK, D3Q15+MRT, D3Q19+MRT, D3Q27+LBGK from top to bottom). The figures in color show the simulated results. The color field indicates the absolute error Eq. (31) where light gray (yellow online) is zero and dark (red online) is $0.3U_{\varphi }(0,0)$. Since the GZS simulations were unstable, we use these spaces to show the predicted secondary flow calculated using Eqs. (23, 24, 25, 26, 27) in grayscale.