Influence of numerical resolution on the dynamics of finite-size particles with the lattice Boltzmann method

We investigate and compare the accuracy and efficiency of different numerical approaches to model the dynamics of finite-size particles using the lattice Boltzmann method (LBM). This includes the standard bounce-back (BB) and the equilibrium interpolation (EI) schemes. To accurately compare the different implementations, we first introduce a boundary condition to approximate the flow properties of an unbounded fluid in a finite simulation domain, taking into account the perturbation induced by a moving particle. We show that this boundary treatment is efficient in suppressing detrimental effects on the dynamics of spherical and ellipsoidal particles arising from the finite size of the simulation domain. We then investigate the performances of the BB and EI schemes in modeling the dynamics of a spherical particle settling under Stokes conditions, which can now be reproduced with great accuracy thanks to the treatment of the exterior boundary. We find that the EI scheme outperforms the BB scheme in providing a better accuracy scaling with respect to the resolution of the settling particle, while suppressing finite-size effects due to the particle discretization on the lattice grid. Additionally, in order to further increase the capability of the algorithm in modeling particles of sizes comparable to the lattice spacing, we propose an improvement to the EI scheme, the complete equilibrium interpolation (CEI). This approach allows us to accurately capture the boundaries of the particle also when located between two fluid nodes. We evaluate the CEI performance in solving the dynamics of an under-resolved particle under analogous Stokes conditions and also for the case of a rotating ellipsoid in a shear flow. Finally, we show that EI and CEI are able to recover the correct flow solutions also at small, but finite, Reynolds number. Adopting the CEI scheme it is not only possible to detect particles with zero lattice occupation, but also to increase up to one order of magnitude the accuracy of the dynamics of particles with a size comparable to the lattice spacing with respect to the BB and the EI schemes.

Figure 1

Sketches of two particles at different resolution on the lattice grid. On the left (a) a particle with radius 2.5 lattice units is represented: momentum exchange is calculated for every fluid node (${\mathbf{x}}_i$) and its adjacent solid nodes (${\mathbf{x}}_b$) along the fluid-solid links (black dashed lines). The inward component along the direction $a$ of the velocity distribution function $f$ interacts with the surface of the particle during the streaming from ${\mathbf{x}}_i$ to ${\mathbf{x}}_b$ and the outward component in the direction $\overline{a}$ is generated. Using standard BB, the interaction point is approximated in the middle position of the fluid-solid link (blue squares), while using an interpolation method the surface of the particle is evaluated in its exact location (red squares). If the particle resolution drops, as shown on the right sketch (b), there might be no lattice nodes enclosed by the solid volume and both BB and EI schemes fail to detect the particle. Using complete equilibrium interpolation it is now possible to detect the particle position and surface (green squares) also between adjacent fluid node links.

Figure 2

Comparison of the slip velocity at the walls $U_s$ normalized with respect to the center channel velocity $U_0$ of a rectangular 2D Poiseuille channel flow for standard LBM (blue diamonds) and for projection-based regularized LBM (PR-LBM, red triangles), as a function of $Kn$ (bottom axis) and $\tau$ (top axis), for a fixed channel width $L=64$. When a strict no-slip boundary condition on the walls is imposed, the slip velocity should be identically zero at any $Kn$ number, but due to the discrete nature of the LBM a numerical slip velocity can appear. A PR-LBM strongly suppresses the numerical slip produced at the channel walls, especially for large values of $Kn$. Here the walls are located at the middle position between the fluid and the solid nodes, so that BB and EI schemes are equivalent.

Figure 3

Sketch of the application of the adaptive boundary condition for the problem of a settling sphere. The analytical value of the velocity field (and of the pressure field) produced by a spherical particle settling at velocity ${\mathit{U}}_P$ is enforced on the boundary nodes in order to emulate the properties of an unbounded fluid. At infinite distance from the particle the flow velocity is the unperturbed one, defined with $u_{\infty }$.

Figure 4

Particle velocity ${\mathit{U}}_P$ normalized with the Stokes velocity $U_0$, defined in Eq. (24), is plotted for no-slip (blue), periodic (green), and adaptive (red) boundary conditions applied on the simulation box of different sizes $L$ with respect to a fixed particle radius $R$. The ideal Stokes settling velocity is recovered only using an adaptive domain boundary that mimics the behavior of an unbounded fluid. Results here represented correspond to $R=10$, $\tau =0.85$, ${\rho }_p=2{\rho }_f$, and $Re={10}^{-3}$.

Figure 5

Terminal velocity $U_P$, averaged over the last 100 time steps, of two vertically aligned spheres settling under gravity as a function of the distance between the centers of the spheres $d$ over their radius $R$. The simulation domain size is fixed to $L=180$, while $\tau =0.6$ and ${\rho }_p=2{\rho }_f$. Results using the adaptive treatment on the simulation domain boundaries (red squares) are compared with the ones obtained using periodic (green squares) and no-slip (blue squares) conditions. For the adaptive case, the general solution from Lamb, truncated to the first term, is applied on the domain boundaries to simulate Stokes settling of two vertically aligned spheres. Results are in fairly good agreement with the predictions from the method of reflections, computed at first reflection, from Eq. (25) (purple solid line), except at the smallest value of $d/R$ where the truncation to the first order shows its limitations. The error bars represent a confidence interval of our measurements and are obtained from the velocity difference, converted in relative units, between the adaptive case and in the no-slip case at $d/R=10$, for which we expect maximum accuracy.

Figure 6

Accuracy of the particle terminal velocity using EI (red triangles) and BB (blue diamonds) boundary conditions for $R=10$, $\tau =0.6$, and $Re={10}^{-3}$ obtained varying the simulation domain size $L$ from 48 to 408. Fit functions obtained using Eq. (26) (colored dashed lines) show a convergence of the error to a limit value (black dashed lines) represented by the fit parameter $\varepsilon$. This limit value represents the maximum accuracy achievable for a given particle resolution.

Figure 7

Relative error in the flow as defined in Eq. (27), $e_u=\sqrt{|u_{\text{an}}\left(U_P\right)-u_{\text{sim}}{|}^2/u_{\text{an}}^2\left(U_P\right)}$, for a less resolved (top) and better resolved (bottom) particle with respectively $R=5$ and $R=20$ at the same confinement ratio $2R/L$. The analytical velocity field $u_{\text{an}}\left(U_P\right)$ is computed using Eqs. (15) and (19). Discretization effects induce higher deviation from the ideal Stokes flow field around the particle that strongly depend on particle resolution.

Figure 8

Velocity $U_P$ of a settling particle with $R=10$, normalized with respect to the Stokes settling velocity $U_0$ using the EI (red solid line) and the BB (blue dashed line) schemes, as a function of time. Other relevant parameters are $L=200,\tau =0.6,{\rho }_p=2{\rho }_f$. In the inset a zoom on the $y$ axis of the selected region is performed to show the fluctuations in the velocity profile due to the volume variation of the particle during its motion.

Figure 9

Asymptotic error (top) of the average velocity $\langle U_P\rangle /U_0$ and average velocity fluctuations $\langle \mathrm{\Delta }{U}_P\rangle /\langle U_P\rangle$ as computed from Eq. (28) (bottom), with respect to particle radius $R$. Convergence rates are obtained through a power-law fit (colored dashed lines) of the measured asymptotic error and reported as the slope of the best fit in the semilog plot. EI strongly reduces errors, keeping an accuracy scaling comparable with the one observed in [22]. Using BB, the accuracy scaling is instead only first order. Second-order convergence is maintained in both cases for the velocity fluctuations.

Figure 10

Number of link nodes for momentum exchange (top) and time evolution of particle velocity (bottom); EI (red triangles), CEI using D3Q19 lattice (blue diamonds), and CEI using D3Q27 lattice (green squares) for a particle with radius $R=0.6$. Standard EI is not able to capture a subgrid particle when initialized in the domain center (zero occupation on top plot), leading to a free fall velocity profile (and thus omitted from the bottom plot for the sake of readability). Using CEI, the particle is correctly detected, but some peaks appear in the velocity profile due to the drop in the number of link nodes during the motion of the particle when a D3Q19 lattice is used. This effect can be mitigated by using the higher order lattice such as D3Q27.

Figure 11

Sketch of particles with $R=1$ initialized with random offsets of $\pm 0.5$ with respect to a particle (black solid line) located at the center of the simulation domain (black square). The colored squares represent the center of each different particle (dashed colored lines). Different initializations lead to different discretizations of the particles on the lattice grid, represented by the number of lattice nodes connected by the black dashed lines that are enclosed by the surface of the particles.

Figure 12

Averages (solid lines) of the accuracy on particle velocity (top) and of the velocity fluctuations (bottom) for $0.6\le R\le 2$. Results are averaged over an ensemble of ten simulations having different random initialization of the particle center within $\pm 0.5$ from the domain center. Minimum and maximum values within the ensemble are represented by the error bars. The CEI scheme (green squares) allows us to describe the dynamics of particles with small radii at increased accuracy with respect to EI (red triangles) and BB (blue diamonds) schemes, while strongly reducing velocity fluctuations. With the CEI scheme it is now possible to detect solid features on the order of one lattice grid spacing.

Figure 13

Time required to perform one call of the momentum exchange routine, averaged over 1000 time steps, for BB (blue diamonds), EI (red triangles), and CEI (green squares) as a function of the particle radius $R$. For small values of $R$ the performances of the three investigated algorithm are similar, while for larger values of $R$ CEI introduces some overhead that can be reduced with a higher level of code optimization.

Figure 14

Analytical prediction given by Eq. (31) of angular frequency as a function of time of an ellipsoid with $a/b=a/c=2$ originally oriented with the major axis $a$ parallel to the flow and a sphere with $a/b=a/c=1$. The $x$ and $y$ axis are normalized with respect to the shear rate $C=2U_{\text{shear}}/L$.

Figure 15

Angular velocity for a sphere of radius $R=10$, rotating in a shear flow with shear rate $C=0.002$. The simulation domain size is set to $L=24$, so that the domain boundaries of the simulation box are almost in contact with the particle. Using the adaptive boundary conditions (red solid line) we are able to suppress detrimental boundary effects and the analytical angular frequency (black dashed line) is exactly reproduced. Using a standard approach (blue point-solid line) consisting of no-slip moving walls in the orthogonal direction to the flow and periodic boundaries along the flow the detrimental effects due to the vicinity of walls on Jeffery's ideal dynamics is evident.

Figure 16

Relative error on the analytical period, given by Eq. (33) and measured by fitting our data with Eqs. (36) and (37), of the oscillations of the angular frequency for an ellipsoid with aspect ratio $a/b=a/c=2$, at fixed major radius $a=10$ for different values of the simulation domain size $L$. We compare the results using the BB (blue diamonds) and EI (red triangles) schemes for two different treatments of the simulation domain boundaries. The influence on the dynamics of the particle of the simulation boundaries is strongly reduced using our adaptive treatment. Results for $\tau =0.6$, $U_{\text{shear}}=0.01$.

Figure 17

Angular frequency evaluated from Eq. (31) using BB (blue diamonds), EI (red triangles), and CEI (green squares) schemes for an ellipsoidal particle with major radius $a=1.5$ and minor radii $b=c=0.75$ at $\tau =0.6$ and $L=30$. The particle center is initialized at the point (15,15,15). For the BB case the algorithm is not able to distinguish the ellipsoid from a spherical particle, while the EI and CEI schemes are increasingly better in capturing shape features. The best fits, obtained using Eqs. (36) and (37), are used to evaluated the accuracy of every method in reproducing the ideal rotational dynamics.

Figure 18

Error on the period (top plot) and error on the trajectory of the angular frequency (bottom plot) of an ellipsoid rotating in a shear flow. The ellipsoid has internal aspect ratio $a/b=a/c=2$, where $a$ is the major axis and $b=c$ are the minor axes. The resolution of the ellipsoid ranges between $0.5\le b\le 2.5$. Results are averaged over an ensemble of ten simulations, each with a different offset of the position of the particle from the domain center. The confinement ratio between particle major axis and simulation domain is kept constant, with $2a/L=0.1$. While EI (red triangles) provides better accuracy with respect to BB (blue diamonds), both methods fail in detecting a particle with $a=1$ and $b=0.5$, where CEI (green squares) succeeds. The latter leads to a general improved accuracy with respect to EI. The error bars represent the minimum and maximum value in the ensemble.

Figure 19

Schematic representation of (a) a cylindrical particle of radius $R$ released near one wall of a two-dimensional channel set into motion by a gravity force $\mathbf{g}$, and (b) a cylindrical particle free to rotate in a shear flow.

Figure 20

Settling trajectory of a cylinder released at $x=0.19L$ of a vertical channel at $Re=1.03$ and $Re=8.33$. The particle radius is $R=12$, while the channel width $L=8R$. The horizontal migration towards the center of the channel is an effect related to the inertia of the particle. The dynamics observed by Li et al. [48] (black dashed lines) is correctly reproduced with all the investigated boundary schemes (colored point lines). Other relevant simulation parameters are $\tau =0.6$ and ${\rho }_p=2{\rho }_f$.

Figure 21

Average velocity fluctuations $\langle \mathrm{\Delta }{U}_P\rangle$ [as defined in Eq. (28)] of the terminal velocity $U_P$ of a cylinder settling in a vertical channel as a function of its radius $R$ for BB (blue diamonds), EI (red triangles), and CEI (green squares). The particle is released at $x=0.19L$ of a channel with width $L=8R$ at $Re=1.03$. EI and CEI slightly reduce the velocity fluctuations related to the volume fluctuations of the particle, but the effect is noticeable only for small values of $R$.

Figure 22

The $\widehat{x}$ (dashed lines) and $\widehat{y}$ (solid lines) velocity components $U_x$ and $U_y$ of the flow versus $y/R$ at the domain section $x/R=-1.8467742$ from a rotating cylinder with radius $R=64$ centered at (0,0). The Reynolds number is $Re=79.6$. The same simulation is performed using BB (blue), EI (red), and CEI (green) and numerical results are in good agreement with the experimental data from Zettner and Yoda [52] (black circles and diamonds) for all the investigated schemes.