A comparative study of fluid-particle coupling methods for fully resolved lattice Boltzmann simulations

TL;DR

This paper systematically compares two fluid-particle coupling methods in lattice Boltzmann simulations, highlighting their accuracy differences and emphasizing the importance of resolution choice for reliable fluid-solid suspension modeling.

Contribution

It provides a detailed comparison of the momentum exchange and partially saturated cells methods, including their accuracy and resolution requirements for fully resolved particulate simulations.

Findings

Momentum exchange method is more accurate for streamwise velocity.

Partially saturated cells method better predicts spanwise velocity.

Resolution of at least 24 cells per sphere diameter is needed for errors below 5%.

Abstract

The direct numerical simulation of particulate systems offers a unique approach to study the dynamics of fluid-solid suspensions by fully resolving the submerged particles and without introducing empirical models. For the lattice Boltzmann method, different variants exist to incorporate the fluid-particle interaction into the simulation. This paper provides a detailed and systematic comparison of two different methods, namely the momentum exchange method and the partially saturated cells method by Noble and Torczynski. Three subvariants of each method are used in the benchmark scenario of a single heavy sphere settling in ambient fluid to study their characteristics and accuracy for particle Reynolds numbers from 185 up to 365. The sphere must be resolved with at least 24 computational cells per diameter to achieve velocity errors below 5%. The momentum exchange method is found to be more accurate in predicting the streamwise velocity component whereas the partially saturated cells method is more accurate in the spanwise components. The study reveals that the resolution should be chosen with respect to the coupling dynamics, and not only based on the flow properties, to avoid large errors in the fluid-particle interaction.

Tables

Table 1. Table 1 : Results obtained for the LBM with different coupling methods for case A, with nominal Galileo number G ​ a = 144 𝐺 𝑎 144 Ga=144 . The error is computed with respect to the results of the reference case AL in [ 48 ] .

	$u_{pV}$	$ℰ\left(u_{pV}\right)$	$L_r$	$ℰ\left(L_r\right)$
Reference	$-1.285$		$1.383$
A-BB-18	$-1.2161$	$0.0536$	$1.4126$	$0.0214$
A-BB-24	$-1.2386$	$0.0361$	$1.3884$	$0.0039$
A-BB-36	$-1.2487$	$0.0283$	$1.3701$	$0.0093$
A-BB-48	$-1.2538$	$0.0243$	$1.3634$	$0.0142$
A-CLI-18	$-1.2150$	$0.0545$	$1.3462$	$0.0266$
A-CLI-24	$-1.2383$	$0.0363$	$1.3616$	$0.0155$
A-CLI-36	$-1.2538$	$0.0243$	$1.3540$	$0.0210$
A-CLI-48	$-1.2586$	$0.0205$	$1.3511$	$0.0231$
A-MR-18	$-1.2172$	$0.0528$	$1.4919$	$0.0787$
A-MR-24	$-1.2510$	$0.0265$	$1.4248$	$0.0303$
A-MR-36	$-1.2643$	$0.0161$	$1.3775$	$0.0040$
A-MR-48	$-1.2646$	$0.0159$	$1.3610$	$0.0159$
A-M1B1-18	$-1.1646$	$0.0937$	$1.4685$	$0.0618$
A-M1B1-24	$-1.1924$	$0.0721$	$1.4330$	$0.0361$
A-M1B1-36	$-1.2182$	$0.0520$	$1.4015$	$0.0134$
A-M1B1-48	$-1.2299$	$0.0429$	$1.3873$	$0.0031$
A-M2B2-18	$-1.2363$	$0.0379$	$1.3820$	$0.0008$
A-M2B2-24	$-1.2490$	$0.0280$	$1.3665$	$0.0119$
A-M2B2-36	$-1.2588$	$0.0204$	$1.3521$	$0.0224$
A-M2B2-48	$-1.2616$	$0.0182$	$1.3482$	$0.0251$
A-M3B2-18	$-1.2128$	$0.0562$	$1.4128$	$0.0216$
A-M3B2-24	$-1.2309$	$0.0421$	$1.3876$	$0.0033$
A-M3B2-36	$-1.2453$	$0.0309$	$1.3692$	$0.0099$
A-M3B2-48	$-1.2513$	$0.0262$	$1.3612$	$0.0158$

Table 2. Table 2 : Results obtained for the LBM with different coupling methods for case B, with nominal Galileo number G ​ a = 178.46 𝐺 𝑎 178.46 Ga=178.46 . The error is computed with respect to the results of the reference case BL in [ 48 ] . No steady solution could be obtained for simulations marked with (*).

	$u_{pV}$	$ℰ\left(u_{pV}\right)$	$u_{pH}$	$ℰ\left(u_{pH}\right)$	${\omega }_{pH}$	$ℰ\left({\omega }_{pH}\right)$	$L_r$	$ℰ\left(L_r\right)$
Reference	$-1.356$		$0.1245$		$0.0137$		$1.629$
B-BB-24	$-1.3097$	$0.0341$	$0.0494$	$0.0554$	$0.0605$	$0.0345$	$1.6593$	$0.0186$
B-BB-36	$-1.3265$	$0.0218$	$0.0652$	$0.0437$	$0.0490$	$0.0260$	$1.6324$	$0.0021$
B-BB-48	$-1.3325$	$0.0173$	$0.0722$	$0.0386$	$0.0367$	$0.0170$	$1.6181$	$0.0067$
B-CLI-24	$-1.3056$	$0.0371$	$0.0883$	$0.0267$	$0.0230$	$0.0069$	$1.6023$	$0.0164$
B-CLI-36	$-1.3267$	$0.0216$	$0.0934$	$0.0229$	$0.0175$	$0.0028$	$1.6058$	$0.0143$
B-CLI-48	$-1.3328$	$0.0171$	$0.0991$	$0.0187$	$0.0104$	$0.0025$	$1.6025$	$0.0163$
B-MR-24(*)	$-1.3037$	$0.0385$	$0.1178$	$0.0049$	$0.0092$	$0.0033$	$1.6625$	$0.0206$
B-MR-36(*)	$-1.3358$	$0.0149$	$0.1097$	$0.0109$	$0.0038$	$0.0073$	$1.6486$	$0.0120$
B-MR-48(*)	$-1.3384$	$0.0130$	$0.1108$	$0.0101$	$0.0013$	$0.0091$	$1.6216$	$0.0046$
B-M1B1-24	$-1.2455$	$0.0815$	$0.1171$	$0.0055$	$0.0332$	$0.0144$	$1.7104$	$0.0500$
B-M1B1-36	$-1.2795$	$0.0564$	$0.1153$	$0.0068$	$0.0220$	$0.0061$	$1.6634$	$0.0211$
B-M1B1-48	$-1.2944$	$0.0454$	$0.1159$	$0.0063$	$0.0186$	$0.0036$	$1.6432$	$0.0087$
B-M2B2-24	$-1.3040$	$0.0384$	$0.1329$	$0.0062$	$0.0437$	$0.0221$	$1.6225$	$0.0040$
B-M2B2-36	$-1.3218$	$0.0252$	$0.1238$	$0.0005$	$0.0255$	$0.0087$	$1.6047$	$0.0149$
B-M2B2-48	$-1.3278$	$0.0208$	$0.1209$	$0.0027$	$0.0198$	$0.0045$	$1.5972$	$0.0195$
B-M3B2-24	$-1.2904$	$0.0484$	$0.1150$	$0.0070$	$0.0223$	$0.0064$	$1.6497$	$0.0127$
B-M3B2-36	$-1.3112$	$0.0331$	$0.1128$	$0.0086$	$0.0135$	$0.0001$	$1.6233$	$0.0035$
B-M3B2-48	$-1.3192$	$0.0271$	$0.1124$	$0.0089$	$0.0108$	$0.0021$	$1.6120$	$0.0104$

Table 3. Table 3 : Results obtained for the LBM with different coupling methods for case C, with nominal Galileo number G ​ a = 190 𝐺 𝑎 190 Ga=190 . The error is computed with respect to the results of the reference case CL in [ 48 ] .

	${\overline{u}}_{pV}$	${\overline{u}}_{pH}$	${\overline{\omega }}_{pH}$	$u_{pV}^{\prime }$	$u_{pH}^{\prime }$	${\omega }_{pH}^{\prime }$	$f$
Reference	$-1.376$	$0.136$	$0.012$	$0.008$	$0.033$	$0.008$	$0.071$
C-BB-36	$-1.3479$	$0.0624$	$0.0452$	$0.0108$	$0.0624$	$0.0439$	$0.0415$
${ℰ}^{\left(\text{C-BB-36}\right)}$	$0.0206$	$0.0536$	$0.0241$	$0.0018$	$0.0210$	$0.0258$	$0.4152$
C-BB-48	$-1.3540$	$0.0678$	$0.0351$	$0.0126$	$0.0678$	$0.0350$	$0.0356$
${ℰ}^{\left(\text{C-BB-48}\right)}$	$0.0162$	$0.0496$	$0.0167$	$0.0030$	$0.0250$	$0.0193$	$0.4984$
C-CLI-36	$-1.3461$	$0.1001$	$0.0284$	$0.0053$	$0.0333$	$0.0120$	$0.0637$
${ℰ}^{\left(\text{C-CLI-36}\right)}$	$0.0220$	$0.0262$	$0.0118$	$0.0022$	$0.0001$	$0.0026$	$0.1030$
C-CLI-48	$-1.3527$	$0.1076$	$0.0195$	$0.0058$	$0.0307$	$0.0087$	$0.0668$
${ℰ}^{\left(\text{C-CLI-48}\right)}$	$0.0172$	$0.0207$	$0.0054$	$0.0019$	$0.0020$	$0.0003$	$0.0599$
C-MR-36	$-1.3547$	$0.1091$	$0.0143$	$0.0115$	$0.0538$	$0.0143$	$0.0636$
${ℰ}^{\left(\text{C-MR-36}\right)}$	$0.0158$	$0.0197$	$0.0016$	$0.0023$	$0.0148$	$0.0043$	$0.1042$
C-MR-48	$-1.3590$	$0.1178$	$0.0080$	$0.0078$	$0.0362$	$0.0080$	$0.0667$
${ℰ}^{\left(\text{C-MR-48}\right)}$	$0.0126$	$0.0134$	$0.0030$	$0.0004$	$0.0021$	$0.0003$	$0.0607$
C-M1B1-36	$-1.2978$	$0.1295$	$0.0223$	$0.0015$	$0.0065$	$0.0019$	$0.0680$
${ℰ}^{\left(\text{C-M1B1-36}\right)}$	$0.0571$	$0.0049$	$0.0074$	$0.0050$	$0.0195$	$0.0047$	$0.0421$
C-M1B1-48	$-1.3131$	$0.1306$	$0.0173$	$0.0013$	$0.0049$	$0.0013$	$0.0669$
${ℰ}^{\left(\text{C-M1B1-48}\right)}$	$0.0459$	$0.0041$	$0.0038$	$0.0052$	$0.0207$	$0.0051$	$0.0576$
C-M2B2-36	$-1.3389$	$0.1363$	$0.0250$	$0.0015$	$0.0072$	$0.0018$	$0.0690$
${ℰ}^{\left(\text{C-M2B2-36}\right)}$	$0.0272$	$0.0001$	$0.0094$	$0.0050$	$0.0191$	$0.0047$	$0.0287$
C-M2B2-48	$-1.3455$	$0.1355$	$0.0185$	$0.0013$	$0.0052$	$0.0013$	$0.0712$
${ℰ}^{\left(\text{C-M2B2-48}\right)}$	$0.0224$	$0.0005$	$0.0046$	$0.0051$	$0.0205$	$0.0052$	$0.0024$
C-M3B2-36	$-1.3283$	$0.1245$	$0.0117$	$0.0016$	$0.0086$	$0.0021$	$0.0733$
${ℰ}^{\left(\text{C-M3B2-36}\right)}$	$0.0349$	$0.0084$	$0.0003$	$0.0049$	$0.0181$	$0.0046$	$0.0328$
C-M3B2-48	$-1.3371$	$0.1240$	$0.0084$	$0.0017$	$0.0095$	$0.0021$	$0.0729$
${ℰ}^{\left(\text{C-M3B2-48}\right)}$	$0.0285$	$0.0089$	$0.0027$	$0.0048$	$0.0174$	$0.0045$	$0.0273$

Table 4. Table 4 : Results obtained for the LBM with different coupling methods for case D, with nominal Galileo number G ​ a = 250 𝐺 𝑎 250 Ga=250 . The error is computed with respect to the results of the reference case DL in [ 48 ] .

	$⟨u_{pV}⟩$	${⟨u_{pV}^{\prime \prime }{u}_{pV}^{\prime \prime }⟩}^{1/2}$	${⟨u_{pr}^{\prime \prime }{u}_{pr}^{\prime \prime }⟩}^{1/2}$	${⟨{\omega }_{pV}^{\prime \prime }{\omega }_{pV}^{\prime \prime }⟩}^{1/2}$	${⟨{\omega }_{px}^{\prime \prime }{\omega }_{px}^{\prime \prime }⟩}^{1/2}$
Reference	$-1.4604$	$0.0087$	$0.0854$	$0.0013$	$0.0067$
D-BB-36	$-1.4111$	$0.0049$	$0.0342$	$0.0005$	$0.0200$
${ℰ}^{\left(\text{D-BB-36}\right)}$	$0.0338$	$0.0026$	$0.0351$	$0.0005$	$0.0091$
D-CLI-36	$-1.4114$	$0.0075$	$0.0701$	$0.0016$	$0.0358$
${ℰ}^{\left(\text{D-CLI-36}\right)}$	$0.0336$	$0.0008$	$0.0105$	$0.0002$	$0.0199$
D-MR-36	$-1.4187$	$0.0037$	$0.0348$	$0.0008$	$0.0116$
${ℰ}^{\left(\text{D-MR-36}\right)}$	$0.0286$	$0.0034$	$0.0347$	$0.0004$	$0.0033$
D-M1B1-36	$-1.3701$	$0.0109$	$0.0913$	$0.0008$	$0.0160$
${ℰ}^{\left(\text{D-M1B1-36}\right)}$	$0.0618$	$0.0015$	$0.0041$	$0.0003$	$0.0064$
D-M2B2-36	$-1.4093$	$0.0102$	$0.1045$	$0.0009$	$0.0209$
${ℰ}^{\left(\text{D-M2B2-36}\right)}$	$0.0350$	$0.0011$	$0.0131$	$0.0003$	$0.0097$
D-M3B2-36	$-1.4010$	$0.0066$	$0.0881$	$0.0158$	$0.0082$
${ℰ}^{\left(\text{D-M3B2-36}\right)}$	$0.0407$	$0.0014$	$0.0019$	$0.0099$	$0.0010$

Equations

$${\tilde{f}}_q(\mathbf{x},t)=f_q(\mathbf{x},t)+{\mathcal{C}}_q(\mathbf{x},t)+{\mathcal{F}}_q(\mathbf{x},t),$$

$$f_q(\mathbf{x}+{\mathbf{c}}_q\Delta t,t+\Delta t)={\tilde{f}}_q(\mathbf{x},t).$$

$$f_q^{\text{eq}}({\rho }_f,\mathbf{U})=w_q\left({\rho }_f+{\rho }_0\left(\frac{{\mathbf{c}}_q\cdot \mathbf{U}}{c_s^2}+\frac{({\mathbf{c}}_q\cdot \mathbf{U}{)}^2}{2c_s^4}-\frac{\mathbf{U}\cdot \mathbf{U}}{2c_s^2}\right)\right)$$

$${\rho }_f(\mathbf{x},t)=\sum _q{f}_q(\mathbf{x},t),\mathbf{U}(\mathbf{x},t)=\frac{1}{{\rho }_0}\sum _q{f}_q(\mathbf{x},t){\mathbf{c}}_q.$$

$${\mathcal{C}}_q^{\text{BGK}}(\mathbf{x},t)=\frac{\Delta t}{\tau }\left(f_q^{\text{eq}}({\rho }_f,\mathbf{U})-f_q(\mathbf{x},t)\right).$$

$$\nu =(\tau -\frac{\Delta t}{2})c_s^2.$$

$${\mathcal{F}}_q(\mathbf{x},t)=\Delta t{w}_q\left[\frac{{\mathbf{c}}_q-\mathbf{U}}{c_s^2}+\frac{{\mathbf{c}}_q\cdot \mathbf{U}}{c_s^4}{\mathbf{c}}_q\right]\cdot {\mathbf{f}}^{\text{ext}},$$

$$\mathbf{u}(\mathbf{x},t)=\mathbf{U}(\mathbf{x},t)+\frac{\Delta t}{2{\rho }_0}{\mathbf{f}}^{\text{ext}},$$

$${\mathcal{C}}_q^{\text{TRT}}(\mathbf{x},t)={\lambda }_{+}(f_q^{+}(\mathbf{x},t)-f_q^{\text{eq},+}({\rho }_f,\mathbf{U}))+{\lambda }_{-}(f_q^{-}(\mathbf{x},t)-f_q^{\text{eq},-}({\rho }_f,\mathbf{U})).$$

$$f_q^{+}=\frac{1}{2}(f_q+f_{\bar{q}}),f_q^{-}=\frac{1}{2}(f_q-f_{\bar{q}}),$$

$${\Lambda }_{\pm }:=\left(\frac{1}{2}+\frac{1}{{\lambda }_{+}}\right)\left(\frac{1}{2}+\frac{1}{{\lambda }_{-}}\right),$$

$$f_{\bar{q}}(\mathbf{x},t+\Delta t)={\tilde{f}}_q(\mathbf{x},t)-2\frac{w_q}{c_s^2}{\rho }_0\mathbf{v}({\mathbf{x}}_b,t)\cdot {\mathbf{c}}_q,$$

$$\mathbf{v}(\mathbf{x},t)={\mathbf{U}}_p(t)+{\mathbf{\Omega }}_p(t)\times (\mathbf{x}-{\mathbf{X}}_p(t)),$$

$$f_{\bar{q}}(\mathbf{x},t+\Delta t)={\tilde{f}}_q(\mathbf{x},t)+{\kappa }_0{\tilde{f}}_q(\mathbf{x}-{\mathbf{c}}_q\Delta t,t)-{\kappa }_0{\tilde{f}}_{\bar{q}}(\mathbf{x},t)-\alpha \frac{w_q}{c_s^2}{\rho }_0\mathbf{v}({\mathbf{x}}_b,t)\cdot {\mathbf{c}}_q,$$

$$\displaystyle f_{\bar{q}}(\boldsymbol{x},t+\Delta t)=\$$

$${\bar{\kappa }}_{-1}{\tilde{f}}_{\bar{q}}(\mathbf{x},t)+{\bar{\kappa }}_{-2}{\tilde{f}}_{\bar{q}}(\mathbf{x}-{\mathbf{c}}_q\Delta t,t)-\alpha \frac{w_q}{c_s^2}{\rho }_0\mathbf{v}({\mathbf{x}}_b,t)\cdot {\mathbf{c}}_q+f_q^{\text{pc}},$$

$${\mathbf{F}}^{q_{\text{f-s}}}({\mathbf{x}}_b,t)=\frac{(\Delta x{)}^3}{\Delta t}\left[({\mathbf{c}}_{q_{\text{f-s}}}-\mathbf{v}({\mathbf{x}}_b,t)){\tilde{f}}_{q_{\text{f-s}}}(\mathbf{x},t)-({\mathbf{c}}_{{\bar{q}}_{\text{f-s}}}-\mathbf{v}({\mathbf{x}}_b,t))f_{{\bar{q}}_{\text{f-s}}}(\mathbf{x},t+\Delta t)\right].$$

$$\mathbf{F}(t)$$

$$\mathbf{T}(t)$$

$$f_q({\mathbf{x}}_{\text{new}},t)=3f_q({\mathbf{x}}_{\text{new}}+{\mathbf{c}}_{q_n}\Delta t,t)-3f_q({\mathbf{x}}_{\text{new}}+2{\mathbf{c}}_{q_n}\Delta t,t)+f_q({\mathbf{x}}_{\text{new}}+3{\mathbf{c}}_{q_n}\Delta t,t).$$

$${\mathcal{C}}_q^{\text{PSM}}(\mathbf{x},t)=(1-B(\mathbf{x},t)){\mathcal{C}}_q^{\text{BGK}}(\mathbf{x},t)+\sum _s{B}_s(\mathbf{x},t){\mathcal{C}}_{q,s}^{\text{solid}}(\mathbf{x},t).$$

$$\text{B1: }{B}_s(\mathbf{x},t)$$

$$\text{B2: }{B}_s(\mathbf{x},t)$$

$$\text{M1: }{\mathcal{C}}_{q,s}^{\text{solid}}(\mathbf{x},t)$$

$$\text{M2: }{\mathcal{C}}_{q,s}^{\text{solid}}(\mathbf{x},t)$$

$$\text{M3: }{\mathcal{C}}_{q,s}^{\text{solid}}(\mathbf{x},t)$$

$${\tilde{f}}_q(\mathbf{x},t)=f_q(\mathbf{x},t)+{\mathcal{C}}_q^{\text{PSM}}(\mathbf{x},t)+(1-B(\mathbf{x},t)){\mathcal{F}}_q(\mathbf{x},t),$$

$$\mathbf{u}(\mathbf{x},t)=(1-B(\mathbf{x},t))\left(\mathbf{U}(\mathbf{x},t)+\frac{\Delta t}{2{\rho }_0}{\mathbf{f}}^{\text{ext}}\right)+\sum _s{B}_s(\mathbf{x},t){\mathbf{v}}_s(\mathbf{x},t).$$

$$\mathbf{F}(t)$$

$$\mathbf{T}(t)$$

$$C=\frac{\frac{1}{|\mathbf{a}|}{\mathbf{a}}^{\top }({\mathbf{F}}_d+{\mathbf{F}}_b)}{3\pi {\rho }_0\nu D\bar{u}},$$

$$\bar{u}=\frac{{\mathbf{a}}^{\top }}{|\mathbf{a}|L^3}\sum _{\mathbf{x}\in \Omega }(\Delta x{)}^3\mathbf{u}(\mathbf{x}).$$

$$Ga=\frac{\sqrt{|\frac{{\rho }_p}{{\rho }_f}-1|g{D}^3}}{\nu },$$

Full text

Keywords: lattice Boltzmann method; particulate flow; direct numerical simulation; fluid-structure interaction; momentum exchange method; partially saturated cells method

1 Introduction

To study particulate flows, direct numerical simulations (DNS) have become a viable and important tool. Resolving the flow structures and the submerged particles offers the unique possibility to trace the motion of single particles, to evaluate the hydrodynamic forces acting on each individual particle, and to investigate the flow field in detail. As an alternative to conventional computational fluid dynamics (CFD) methods that solve the Navier-Stokes equations, the lattice Boltzmann method (LBM) has been applied successfully for DNS studies of such particulate systems. Examples include coupled simulations of several thousand [34, 46] up to millions of spherical particles [54, 15] as they typically appear in fluidized beds at a laboratory scale. Also non-spherical objects can be handled by this method which allows to study flows with blood cells [20], the motion of elongated particles [28, 2] and ice floes [29]. Extensions exist that incorporate electrostatic forces between the particles [3, 23], enable deformable objects [53, 22] or model self-propelled swimmers [38].

A crucial component of such simulations is the accurate and efficient description of the interaction between the solid and the fluid phase which requires appropriate coupling methods. Avoiding the costly remeshing due to changes in the geometry and topology that are caused by the moving particles is a key aspect here. Therefore, classical CFD methods often make use of the immersed boundary method to satisfy the no-slip boundary condition on the particle surfaces and to compute the hydrodynamic forces acting on them [37, 47, 44]. The standard LBM employs uniform Cartesian meshes with cubic cells and thus it is natural here that coupling methods are utilized which do not alter the mesh. This results in an excellent parallel performance and scalability of the LBM method [13, 16] and with suitable data structures also the fluid-solid coupling can be implemented with high parallel efficiency [15].

Several different fluid-solid coupling approaches have been proposed for the LBM: the momentum exchange method [25, 1], the partially saturated cells method [31, 19], and methods that rely on Lagrangian marker points like the immersed boundary method [9, 30, 51] or the external boundary force method [53]. Beyond this general classification of the methods, a large number of subvariants can be found in literature.

It is apparent that a rigorous systematic comparison of these methods for systems with moving solid particles is needed to assist the choice for a suitable method. Challenges here are to define appropriate benchmark scenarios for which reliable reference data exist and to cover the flow regimes at larger Reynolds numbers as they are most relevant for typical engineering applications. One standard benchmark is to evaluate the drag force acting on an infinite periodic array of spheres in Stokes flow, as e.g. in [36, 5, 21, 8]. This benchmark can make use of an analytic solution in the form of a series expansion [17, 41] and is thus particularly suitable to validate implementations and investigate certain characteristics of the fluid-solid coupling methods as will be demonstrated later on in this article. However, the significance of this benchmark is limited when the goal is to simulate scenarios with moving particles at larger Reynolds numbers since the motion of the particles will affect the accuracy and stability of the coupling mechanism. A suitable benchmark has to take into account all possible sources of numerical errors that might appear in a coupled simulation. Besides the fluid simulation itself, errors can be introduced by the particle simulation and the coupling from the fluid to the solid phase and vice versa.

Some studies exist that aim at directly comparing different fluid-solid coupling methods for LBM but they are mostly restricted to two-dimensional setups: In [36], the momentum exchange method with an interpolated bounce back scheme and the immersed boundary method are compared with respect to their accuracy and efficiency for a laminar flow over a stationary cylinder [42]. The settling of circular and elliptical particles is simulated in [49] to establish a comparison between three different LBM collision operators in combination with the momentum exchange method. Various coupling methods are compared in [6] for two-dimensional objects with a prescribed motion and with the focus on the accuracy of acoustic properties. In a recent study, sedimentation of a single circular particle is simulated with different coupling methods [52].

The benchmark proposed by Uhlmann and Dušek [48] features a single sphere settling in an ambient fluid in a three-dimensional domain. High accuracy results obtained with spectral elements are available for flow regimes with particle Reynolds numbers ranging from $185$ up to $365$. This test case can be set up independent of the used CFD method and it is well suited to compare different methods. Its applicability to LBM has already been demonstrated in a preparatory study by [40]. In particular, it can be used to evaluate which resolution is required to achieve sufficiently accurate results in these flow regimes. In this article, we therefore use this benchmark to establish a systematic comparison between two commonly used coupling approaches with the LBM, namely the momentum exchange method and the partially saturated cells method. For each method, three distinct subvariants are applied and validated to illustrate their respective features and strengths. For the momentum exchange method, the simple bounce back and two interpolated versions are compared, whereas different solid collision operators are used for the partially saturated cells method. Methods that employ Lagrangian markers are not included in this study as they are typically more suitable for deformable objects in contrast to the rigid particles here. However, future work should also investigate those methods to establish a complete overview over all available coupling methods.

The remainder of this paper is structured as follows: First the numerical methods are introduced in Sec. 2. This includes a brief summary of the lattice Boltzmann method in Sec. 2.1 as basis for the momentum exchange method and its variants presented in Sec. 2.2. The partially saturated cells method with its variants is introduced in Sec. 2.3. A short study in Sec. 3 compares the obtained force on a fixed sphere in Stokes flow for the resulting six approaches. Section 4.1 elaborates the setup of the benchmark test case from Uhlmann and Dušek. The six coupling approaches are then evaluated for four different flow regimes in Secs. 4.2–4.5. The results are discussed in Sec. 4.6. Finally, the most important findings are summarized in Sec. 5.

2 Numerical method

2.1 Lattice Boltzmann method

For the simulation of hydrodynamics, the lattice Boltzmann approach [7] with the D3Q19 lattice model [39] is utilized. Having its origin in statistical mechanics, the evolution of particle distribution functions (PDFs) on a Cartesian lattice is computed. Each of these PDFs $f_{q}$, with $q\in\{0,\dots,18\}$, is associated with a lattice velocity $\boldsymbol{c}_{q}$. In its most general form, the lattice Boltzmann equation is then given by the collision step

[TABLE]

specified by the collision operator $\mathcal{C}_{q}$ and the external forcing operator $\mathcal{F}_{q}$. In the succeeding stream step, the post collision values $\tilde{f}_{q}$ are distributed to the corresponding neighbor lattice cells via

[TABLE]

The most commonly applied collision model is the BGK model [4] that uses a single relaxation parameter to linearly relax the PDFs towards their equilibrium values $f_{q}^{\text{eq}}$. Those can be computed as

[TABLE]

for incompressible flows [18]. The fluid density $\rho_{f}=\rho_{0}+\delta\rho_{f}$, with the mean density $\rho_{0}$ and the fluctuation $\delta\rho_{f}$, and the velocity $\boldsymbol{U}$ are cell local quantities and calculated via moments of the PDFs:

[TABLE]

The lattice weights $w_{q}$ are as given e.g. in [39] and $c_{s}$ is the lattice speed of sound. The collision operator for the BGK model is then

[TABLE]

It features the relaxation time $\tau\in(\frac{1}{2},\infty)$ which determines the kinematic fluid viscosity $\nu$ via

[TABLE]

The forcing operator in Eq. (1) is used to incorporate external forces and can be written as

[TABLE]

with a constant force density $\boldsymbol{f}^{\text{ext}}$ [24]. The cell local macroscopic velocity $\boldsymbol{u}$ is then obtained via

[TABLE]

and thus differs from the velocity $\boldsymbol{U}$, which is used to calculate $f_{q}^{\text{eq}}$ and $\mathcal{F}_{q}$ in Eqs. (3) and (7), respectively, by a shift depending on $\boldsymbol{f}^{\text{ext}}$ [14].

However, it is well-known that the BGK collision operator has shortcomings with respect to stability and accuracy [27]. In particular, the location of solid boundaries exhibits an undesired $\tau$ dependency. To overcome this shortcoming, the two-relaxation-time collision operator from [12] can be applied, which is given by

[TABLE]

Here, the PDFs and their equilibrium values are split into their symmetric and anti-symmetric parts,

[TABLE]

where $\bar{q}$ denotes the opposite direction of $q$, such that $\boldsymbol{c}_{\bar{q}}=-\boldsymbol{c}_{q}$. This model features two independent relaxation parameters, $\lambda_{+},\lambda_{-}\in(-2,0)$, where the first one is again related to the kinematic viscosity and $\lambda_{+}=-\tfrac{\Delta t}{\tau}$ holds. The boundary location can be rendered effectively independent of the viscosity by choosing $\lambda_{-}$ to satisfy

[TABLE]

where $\Lambda_{\pm}$ is a constant value [12, 21].

In the LBM context, it is common to use lattice units which results in the lattice spacing $\Delta x=1$, the time step size $\Delta t=1$, $\rho_{0}=1$ and $c_{s}=\tfrac{1}{\sqrt{3}}$. To allow for simulations of fluid-solid flows with fully resolved rigid objects, the LBM has to be extended by a suitable coupling method. Two possibilities are briefly explained in the following sections, namely the momentum exchange method in Sec. 2.2 and the partially saturated cells method in Sec. 2.3. All features are implemented in the LBM framework waLBerla222http://walberla.net/ [13] that is used for the simulations in Secs. 3 and 4.

2.2 Momentum exchange method

The momentum exchange method (MEM) for fluid-solid interactions was first presented by Ladd [25] and then extended by Aidun et al. [1]. The main idea is to explicitly map the particle into the fluid domain by marking containing cells as solid and to apply suitable no-slip boundary conditions along its surface. Using the link-based momentum exchange idea, the hydrodynamic force acting on the particle can be obtained. Next, the different parts of this method are presented.

2.2.1 No-slip boundary condition

On the particle’s surface, the macroscopic flow velocity must match the particle’s velocity. In the momentum exchange method, this is realized by imposing no-slip boundary conditions between the fluid and solid cells. The simplest variant of this boundary condition is the bounce back (BB) scheme [25]:

[TABLE]

where $\boldsymbol{v}$ is the particle’s velocity calculated at the boundary location $\boldsymbol{x}_{b}$ midway between fluid and solid cell center. This velocity is obtained by

[TABLE]

with the particle’s translational velocity $\boldsymbol{U}_{p}$, rotational velocity $\boldsymbol{\Omega}_{p}$, and center of gravity $\boldsymbol{X}_{p}$. The BB scheme results in a staircase approximation of the particle shape. However, as the actual shape and thus the exact surface position of the particle is known, more advanced boundary conditions can be applied. They typically use the information about the actual boundary location by introducing a variable $\delta_{q}$ which is the ratio between the distance from the fluid cell center to the boundary and the distance from the fluid cell center to the solid cell center. The boundary location is thus given as $\boldsymbol{x}_{b}=\boldsymbol{x}+\boldsymbol{c}_{q}\delta_{q}$. For example, the central linear interpolation (CLI) scheme from [12] interpolates the PDFs along the direction of $\boldsymbol{c}_{q}$ to realize the no-slip boundary condition:

[TABLE]

with the coefficients $\kappa_{0}=(1-2\delta_{q})/(1+2\delta_{q})$ and $\alpha=4/(1+2\delta_{q})$. It can be observed that this scheme is a generalization of the BB scheme which is recovered by setting $\delta_{q}=\tfrac{1}{2}$.

Another variant of the boundary condition is the multi-reflection (MR) scheme which applies a quadratic interpolation strategy for the PDFs [12]:

[TABLE]

where the expressions for $\kappa$, $\alpha$, and the correction term $f_{q}^{\text{pc}}$ can be found in Tabs. 3 and 4 from [12].

In contrast to other existing interpolation schemes, CLI and MR yield viscosity independent results in combination with the TRT collision operator [33, 8], as will be demonstrated in Sec. 3.

2.2.2 Fluid-solid interaction force

Following [25] and the proposed improvement by [50], the force $\boldsymbol{F}^{q_{\text{f-s}}}$ acting on a solid object at position $\boldsymbol{x}_{b}$ over a fluid-solid link $q_{\text{f-s}}$ can be calculated with the help of the momentum exchange idea:

[TABLE]

This force is thus given as the difference between the momentum towards the boundary, i.e. multiplying the lattice velocity with the post-collision value of the PDF, and the momentum away from the boundary, calculated by multiplying the inverse lattice velocity with the PDF after the boundary handling from Eqs. (12), (14) or (15). To ensure Galilean invariance for moving solid boundaries, it has been shown that the boundary velocity $\boldsymbol{v}(\boldsymbol{x}_{b},t)$ has to be subtracted from the lattice velocities [50]. By summing up all the local contributions, i.e. all corresponding fluid-solid links in the boundary cells of a specific object, the total hydrodynamic force $\boldsymbol{F}$ and torque $\boldsymbol{T}$ acting on this object can be computed via

[TABLE]

2.2.3 PDF reconstruction

As the momentum exchange method requires an explicit mapping of the object into the fluid domain, the cells regarded as solid cells are not updated in the LBM time step, i.e. when executing Eqs. (1) and (2). However, as the particle moves across the computational domain, cells will eventually be converted from solid to fluid and vice versa. This requires a consistent reconstruction of the missing PDF information for the cells that returned to the fluid state. It has been examined that this refilling can lead to oscillations in the interaction forces and the pressure, e.g. in [35] where an overview of several refilling schemes is given. Following those findings, the normal extrapolation refilling with an explicit enforcing of the velocity constraint is applied here. In a first step, the missing PDFs are obtained by extrapolating the PDFs of the neighboring fluid cells in an approximate normal direction $\boldsymbol{c}_{q_{n}}$ of the particle’s surface. This direction corresponds to the lattice direction $q_{n}$ that maximizes $\boldsymbol{n}\cdot\boldsymbol{c}_{q}$ where $\boldsymbol{n}$ is the local surface normal of the particle. Using quadratic extrapolation, the PDFs in a new fluid cell $\boldsymbol{x}_{\text{new}}$ are reconstructed as

[TABLE]

Secondly, the newly initialized PDFs are transformed into moment space where the local fluid velocity can be set explicitly to match the boundary velocity without affecting the other moments [35]. This ensures that the no-slip boundary condition is fulfilled and is a key aspect of the refilling. If such an extrapolation is not possible, as due to other objects or a boundary nearby, less fluid neighbor cells are included into the extrapolation procedure. If no additional fluid cell at all is available in this direction, the refilling is done by setting the PDFs to their equilibrium values $f_{q}^{\text{eq}}(\bar{\rho}_{f},\boldsymbol{v}(\boldsymbol{x}_{\text{new}},t))$ from Eq. (3) based on a local average fluid density $\bar{\rho}_{f}$ and the particle velocity $\boldsymbol{v}$.

2.2.4 MEM algorithm summary

An overview of the complete algorithm is displayed in Alg. 1. Besides the already mentioned parts of the algorithm, it features subcycling loops for the LBM and the rigid body solver part, i.e. to execute this part repeatedly within a single global time step. For the LBM, it was proposed in [26] to use two LBM subcycles to avoid oscillations in the force resulting from the bounce back nature of the boundary treatment. Since thus two different values for the force are calculated, an averaging step after the LBM subcycling is needed. If the bodies are in close contact, as e.g. in sediment beds, it might be beneficial to also use subcycling for the rigid body solver to avoid large collision forces because of large overlaps. During this type of subcycling, the hydrodynamic forces acting on the particles are kept constant.

Summarizing, the momentum exchange coupling algorithm is a convenient and well-established approach for the simulation of particulate flows. All parts of the algorithm can be carried out in a strictly local fashion allowing for an efficient implementation for parallel supercomputers. The different parts profit directly from new developments in the respective areas like more accurate boundary conditions or more stable collision operators. As they are applied in a modular fashion, the user can select different combinations of the components to find the most suitable configuration. However, in all cases we employ an explicit, binary mapping of the particle into the fluid domain as a key concept. This has obvious drawbacks like the necessary PDF reconstruction procedure for newly uncovered cells. Additionally, the number of solid cells that belong to one particle can vary over time. This is one of the reasons why the force acting on a particle as it traverses the computational domain exhibits oscillations that may affect the numerical stability and accuracy of the simulation.

2.3 Partially saturated cells method

The partially saturated cells method (PSM), proposed in [31], is also referred to as Noble-Torczynski’s method or the immersed moving boundary method. It differs from the momentum exchange method significantly, as summarized below. Instead of an explicit mapping of the object, local solid volume fractions are computed for each computational cell. This fill level is then used to modify the LBM collision operator and to evaluate the hydrodynamic forces.

2.3.1 Modified collision operator

The PSM collision operator is constructed as a weighted average of the BGK collision operator $\mathcal{C}_{q}^{\text{BGK}}$, Eq. (5), and a special solid collision operator $\mathcal{C}_{q}^{\text{solid}}$, i.e.

[TABLE]

The local weighting factor $B$ is directly related to the cell local solid volume fraction $\varepsilon$. Both quantities are summations over the individual contributions from all particles $s$ intersecting with the cell, such that $B(\boldsymbol{x},t)=\sum_{s}B_{s}(\boldsymbol{x},t)$ and $\varepsilon(\boldsymbol{x},t)=\sum_{s}\varepsilon_{s}(\boldsymbol{x},t)$. At least two versions for the calculation of $B_{s}$ exist [31], given by

[TABLE]

respectively. Also for $\mathcal{C}_{q,s}^{\text{solid}}$, which is the collision operator only acting when an object $s$ intersects with the cell, different variants have been proposed [31, 19]:

[TABLE]

Here, the velocity of particle $s$, $\boldsymbol{v}_{s}=\boldsymbol{v}_{s}(\boldsymbol{x},t)$ from Eq. (13), is evaluated at the cell center. The modifications applied in the PSM to the original LBM approach also require to incorporate a weighting in the collision step, now given as

[TABLE]

and in the calculation of the macroscopic velocity

[TABLE]

Noble and Torczynski suggested to use either (M1) with (B1) or (M2) with (B2) in [31]. In [19], (M1) has been modified to become (M3) and [43] suggested the usage of (M3) with (B2) for accuracy reasons.

2.3.2 Solid volume fraction computation

A crucial and expensive part of the PSM algorithm is the evaluation of the solid volume fraction $\varepsilon_{s}$ in each cell. When the particles move, $\varepsilon_{s}$ must be recomputed in each time step. As summarized in [32], different variants can be used. Ideally, an analytical formula exists for the intersection volume of a cubic cell and the particle. Unfortunately, this is not the case for almost all particle shapes in three dimensions. Alternatively, the volume can be approximated by subdividing the cell into smaller sampling points and evaluating whether these points are inside or outside of the object. This so called supersampling method, see e.g. [3], works for any shape for which one can determine whether a point is contained inside. Since this method is quite efficient and sufficiently accurate when implemented in a recursive fashion, it is used here in the following. As another approximation, we refer to [10], where the cell’s edges are intersected with the particle and the ratio of covered to total edge length is taken as $\varepsilon_{s}$. As mentioned in [55], the partially saturated cells method can also be used if the minimum normalized distance to the particle surface is employed instead of the actual solid volume fraction.

2.3.3 Fluid-solid interaction force

The hydrodynamic force acting on the particle is calculated with the help of $B_{s}$ and $\mathcal{C}_{q,s}^{\text{solid}}$. It is applied at the cell center $\boldsymbol{x}_{s}$ of all cells that intersect with the particle. By summing up over all lattice directions of all intersecting cells, the total force and torque,

[TABLE]

respectively, can be evaluated.

2.3.4 PSM algorithm summary

The partially saturated cells method is summarized in Alg. 2. The LBM subcycling is omitted here since in our experiments it was not necessary. Even though the description of the algorithm might appear shorter than Alg. 1, the LBM step here is much more complex than in the momentum exchange case, since it involves the computation of $\mathcal{C}_{q,s}^{\text{solid}}$ and the local force contributions. Additional complexity is created when several objects intersect with one single cell since this requires to treat the contributions to $B_{s}$ and $\mathcal{C}_{q,s}^{\text{solid}}$ for each object separately [32].

The advantage of this coupling method is the smooth variation of the forces acting on the particle when it traverses the domain. This originates from the use of the solid volume fraction that varies continuously between cells inside and outside of the particle and thus smoothens the staircase approximation of the momentum exchange method. However, as it essentially relies on the BGK collision operator, the simulations suffer from the same inaccuracies and, in particular, exhibit a $\tau$ dependency as shown in Sec. 3. A thorough analysis of the method is still not available since it is characterized by the special solid collision operator and to our knowledge no improvements of its stability and accuracy are yet available in the literature. Additionally, the calculation of the solid volume fraction is complex and expensive and may thus consume a major portion of the computing time.

3 Force on fixed sphere in Stokes flow

In a short study, we investigate the dependency of the fluid-solid interaction force on a particular choice of the relaxation time $\tau$ and thus the fluid viscosity $\nu$. For that reason, we evaluate the force on a periodic array of spheres in Stokes flow for which an analytical solution in the form of a series expansion exists [17, 41] . Such a study has been carried out earlier for different variants of the MEM approach, e.g. in [36, 5, 21, 8], but to the best of our knowledge no results are available for the PSM yet. The simulation setup features a periodic and cubic domain $\Omega$ of size $L\times L\times L$ with a fixed sphere of diameter $D=L/2$ positioned in its center. This corresponds to the so called SC setup. The flow is driven by an external force density $\boldsymbol{a}=(a,0,0)^{\top}$ with $a=10^{-5}$ which enters Eq. (7) as $\boldsymbol{f}^{\text{ext}}$. The total force acting on the sphere is the sum of the drag force $\boldsymbol{F}_{d}$, evaluated via Eq. (17) or (28), respectively, and the buoyancy force $\boldsymbol{F}_{b}=\tfrac{\pi}{6}D^{3}\boldsymbol{a}$ due to external forcing [5]. The total dimensionless force in forcing direction is then given as

[TABLE]

with the domain average macroscopic velocity

[TABLE]

The simulation is carried out for $L=32$ and various values of $\nu$ with different coupling methods and is run until $C$ converges. For the MEM, the TRT collision operator is applied together with either the BB from Eq. (12), the CLI from Eq. (14), or the MR scheme from Eq. (15) for the no-slip boundary condition along the sphere’s surface. We choose $\Lambda_{\pm}=\tfrac{3}{16}$ which has been shown to yield good flow properties in combination with particles [12, 21]. For the PSM, the commonly used combinations M1B1, M2B2 and M3B2 of Eqs. (21)–(25) are applied.

The results are shown in Fig. 1 as relative error of $C$ to the reference value $C^{\text{(ref)}}=2.8402$ from [41]. The MEM variants yield forces that are independent of the particular choice of $\nu$. The results of CLI and MR are especially accurate whereas BB overestimates the force by less than $5\%$. In contrast, all PSM variants exhibit a strong dependency on the chosen viscosity and the obtained force gets smaller with increasing $\nu$. The behavior of M1B1 and M2B2 is similar with a nearly constant offset in between whereas the dependency is less pronounced for M3B2.

This study reveals one general drawback of the PSM which essentially relies on the BGK collision operator and thus cannot yield viscosity independent results. To the best of our knowledge, no attempts that aim to improve this behavior are reported which is most likely due to the lack of a thorough analysis of the PSM. As expected, this drawback is not present for the MEM when applying the TRT collision operator with a constant $\Lambda_{\pm}$ and suitable boundary conditions [12].

4 Motion of a single heavy sphere in ambient fluid

4.1 Description

In this benchmark scenario, a sphere that is heavier than the surrounding fluid and has a diameter $D$ is exposed to an external flow field. The simulation domain is rectangular and of size $5.34D\times 5.34D\times 16D$ and set up as periodic in the horizontal directions $x$ and $y$. A constant inflow from below with the velocity $\boldsymbol{u}_{\infty}$ is imposed and an outflow condition is applied at the top of the domain. The sphere is initially placed at the position $(2.67D,2.67D,5.34D)$. In order to keep the sphere positioned inside the simulation domain, the net force due to buoyancy and gravity is chosen to balance the acting drag force $F_{\text{drag}}$ in vertical direction. This setup corresponds to the one used in [48] for the simulation with the immersed boundary method from [47]. As a reference, the results from the spectral element method will be used which are also reported in [48].

More specifically, the ratio between solid and fluid density, $\rho_{p}/\rho_{f}$, is fixed to $1.5$ for all simulations. This density ratio allows the definition of an important characteristic quantity, the Galileo number

[TABLE]

with the magnitude of the gravitational acceleration $g=|\boldsymbol{g}|$. The Galileo number characterizes different regimes of the sphere’s motion which will be explored and investigated with this test case. For increasing $Ga$, the setup undergoes transitions from the steady axisymmetric regime to the steady oblique regime and to the oscillating oblique regime until it arrives at the chaotic regime.

However, the Galileo number cannot be set directly in this test case since it depends on the gravitational acceleration which is chosen with respect to the acting drag force as mentioned earlier. This force is mainly determined by the inflow velocity and the fluid viscosity, but is also influenced by the applied fluid-solid coupling algorithm. Therefore, the execution of the test case consists of two parts: At the beginning, the simulation is started with the sphere held fixed at its initial position. The viscosity is initially set to $0.01$ in lattice units and independent of the numerical resolution. The inflow velocity as well as the initial fluid velocity are determined as $\boldsymbol{u}_{\infty}=\left(0,0,Re_{||,\text{ref}}\nu/D\right)^{\top}$ where $Re_{||,\text{ref}}$ is the particle Reynolds number reported in Tabs. 2-5 of [48] for the reference simulations. The simulation is then run until the drag force on the sphere is sufficiently converged. Then, the gravitational acceleration is determined via

[TABLE]

and the Galileo number can be computed via Eq. (32). If its value is close enough to the desired one, the second part of the simulation can be started. In all presented simulations, the relative difference between the actual and the targeted Galileo number is below $10^{-4}$. Else, the viscosity is adjusted via $\nu\leftarrow\tfrac{Ga}{Ga_{\text{target}}}\nu$ and the simulation is continued, repeating the evaluation procedure with the updated viscosity. In our experiments, this routine typically converged after a few iterations and thus avoids tedious manual experimentation.

In the second part of the simulation, the sphere is released and allowed to move freely while maintaining the flow field from the first part. Its position and velocity are updated in each simulation step based on the acting forces via an explicit Euler integration. Unless stated otherwise, all simulations are carried out for $250$ dimensionless time steps, with $t_{\text{ref}}=D/u_{\text{ref}}$ and $u_{\text{ref}}=\sqrt{|\rho_{p}/\rho_{f}-1|gD}$.

As inflow condition, the common velocity bounce-back scheme like in Eq. (12) with the inflow velocity $\boldsymbol{u}_{\infty}$ is applied. The outflow is modeled following [11] with $\rho_{\text{out}}=1$ in lattice units to yield

[TABLE]

To obtain the forces on the sphere, two LBM subcycles are executed and the calculated forces are averaged over these two cycles.

For the actual test case, the four different regimes are simulated by choosing a corresponding Galileo number. Here we consider a range of the Galileo numbers from $144$ up to $250$ which corresponds to particle Reynolds numbers of approximately $185$ to $365$. The quantities of interest for each regime are evaluated, e.g. the magnitudes of the vertical and horizontal relative sphere velocity, $u_{pV}$ and $u_{pH}$, and of the vertical and horizontal rotational velocities, $\omega_{pV}$ and $\omega_{pH}$. The exact definition of all quantities and coordinate systems is given in [48] and the reported values here are made dimensionless by the same scaling. Then, the relative errors of these quantities with respect to the spectral element results are evaluated analogously to [48]. In particular, the errors in velocity components are always evaluated relative to the reference value of $u_{pV}$.

For the MEM from Sec. 2.2, the TRT collision operator with $\Lambda_{\pm}=\tfrac{3}{16}$ with extrapolation reconstruction technique of Eq. (19) is used in all cases. The variants only differ in the no-slip boundary condition along the sphere’s surface which are the BB from Eq. (12), the CLI from Eq. (14) and the MR scheme from Eq. (15). For the PSM from Sec. 2.3, the combinations M1B1, M2B2 and M3B2 of Eqs. (21)–(25) are applied. The respective labels are chosen to reflect the test regime (A-D), the coupling method, and the cells per diameter $D/\Delta x$, where $\Delta x$ is the cell size. Values between $18$ and $48$ are typically used for $D/\Delta x$, which then also determines the overall resolution of the simulation.

4.2 Steady axisymmetric regime

The first regime with a Galileo number of $144$ features a steady and axisymmetric flow around the sphere. The simulation has been carried out for four different resolutions $D/\Delta x\in\{18,24,36,48\}$. The contour plots of the vertical relative flow velocity for $D/\Delta=18$ and $48$ are shown in Fig. 2 for the CLI variant. The results for the vertical sphere velocity and the recirculation length are given in Tab. 1, together with their errors with respect to the reference values from the spectral element results in [48]. Generally, the values for $u_{pV}$ converge against the desired reference value in all cases. The differences in the MEM variants are overall rather small but especially at higher resolutions CLI and MR are more accurate. Regarding the PSM variants, clearly M1B1 performs worse than M2B2 and M3B2. In contrast to that, all LBM coupling methods converge to a recirculation length that is smaller than in the reference case.

In addition to those two values, reference results for the flow velocity profiles along $x_{p\parallel}\in\{-1,-3,-5,-7\}$ and the pressure coefficient $c_{p}$ along a circle around the sphere are reported in [48]. To evaluate the data in a compact form, they are compressed into single root mean square (RMS) values of the absolute deviation from these values. For a quantity $\phi$, the RMS error is here defined as

[TABLE]

The $N$ evaluation points along the axis $x_{\text{coord}}$ are prescribed by the available reference data. The RMS errors for $u_{r\parallel}$ and $u_{r\perp}$ along the $x_{p\perp}$ axis and at the positions $x_{p\parallel}=-1$ and $-7$, denoted by $|_{-1}$ and $|_{-7}$ respectively, are given in the graphs of Fig. 3. Additionally, the RMS errors of $c_{p}$ along a circle near the sphere’s surface and evaluated at different angles $\theta_{1}$ are shown in Fig. 4. The differences between the MEM variants are again rather small. The MR scheme shows overall slightly better results than CLI, particularly for the higher resolution cases of $u_{r\parallel}$. BB is less accurate than CLI and MR in almost all cases, especially for $c_{p}$. M2B2 is the most accurate PSM variant for $u_{r\parallel}$ and overall similar to M3B2 for $u_{r\perp}$ and $c_{p}$. The errors of M1B1 are up to four times larger for $u_{r\parallel}$ and for $u_{r\perp}$ at $x_{p\parallel}=-1$ in comparisons to the other PSM combinations.

Summarizing, all MEM variants are able to capture the characteristics of this regime very well. The outcomes of M2B2 and M3B2 are comparable to the ones from the MEM variants but M1B1 clearly falls behind the other two variants.

4.3 Steady oblique regime

When the Galileo number is increased to a nominal value of $178.46$, which corresponds to a particle Reynolds number of approximately $243$, the sphere eventually starts to move in the horizontal plane with a constant horizontal velocity. To trigger this movement stemming from an instability, the stabilizing numerical symmetry of this system is broken by slightly perturbing the sphere’s initial horizontal position. Note, however, that this does not change the physical setting due to the periodic boundaries in horizontal directions. In Fig. 5, the contour lines of the projected flow velocities are shown for two different planes. As a result of the horizontal movement, the horizontal and rotational sphere velocities, $u_{pH}$ and $\omega_{pH}$, can additionally be used for a comparison. The values and errors of these four quantities are given in Tab. 2. Following the reasoning given in [48], the relative errors for quantities triggered by instabilities are calculated with respect to the reference value of $u_{pV}$ in all cases. For our setup, the simulations become unstable for the lowest resolution of $D/\Delta x=18$ which are therefore excluded from the table. Regarding the vertical sphere velocity, the MEM variants with errors below $2\%$ for the finest resolution perform systematically better than the PSM variants. The opposite is the case for $u_{pH}$ where the PSM variants exhibit very small errors almost independent of the resolution. The results obtained by CLI and MR are still acceptable but BB shows significant deviations from the reference value. The horizontal rotational velocity is captured reasonably well by all methods with errors below $1\%$ for resolutions of $36$ and $48$ except for BB. Like for case A, all methods tend to have a shorter recirculation length. It has to be noted, however, that the results from MR show significant oscillations over the simulation time and cannot be regarded as steady state results. This explains the fluctuations in the results for the different resolutions.

Due to the steady nature of the flow, profiles of various velocity components can again be compared and will be expressed by RMS errors with the definition from Eq. (35). Graphs of some representative errors over the resolution are shown in Fig. 6. Once more, it can be seen that BB for $\mathcal{E}_{\text{RMS}}^{x_{p\perp}}(u_{r\perp})$ and M1B1 for $\mathcal{E}_{\text{RMS}}^{x_{pHz\perp}}(u_{r\parallel})$ fall behind the accuracy of the other methods. M2B2 and M3B2 yield the overall best results followed by CLI. The large errors in MR for the lowest resolution can again be explained by the unsteady results there. When analyzing the RMS errors of $c_{p}$ along $\theta_{1}$ and $\theta_{2}$, as done in Fig. 7, BB performs significantly worse than the rest. The results of all other methods are comparable with slightly more accurate results for PSM at low resolutions.

To summarize the findings for this regime, again M1B1 and now also BB show significantly worse results compared to the other variants. The MEM variants seem to be better suited for predicting the vertical sphere velocity whereas the PSM variants work very well for the horizontal sphere velocity. The flow properties are best represented by M2B2 followed by CLI and M3B2. The higher order MR scheme seems to destabilize the system as the results show oscillations which are expected to occur only at larger Galileo numbers like the one used for the next regime.

4.4 Oscillating oblique regime

For a Galileo number of $190$ ($Re_{\parallel}\approx 263$), the sphere is expected to exhibit an oscillating motion with a constant frequency $f$. As explained in [48], this behavior is triggered by a secondary bifurcation and is only present in a rather narrow range of Galileo numbers. Numerical inaccuracies can lead to a shift of the bifurcation point and the respective method will then fail to capture this motion at $Ga=190$. In Fig. 8, a phase-space diagram of the results for the different coupling methods together with the reference data is shown for the two resolutions $D/\Delta x=36$ and $48$. The expected time-periodic behavior is a closed curve around a fixed midpoint. Even for the finer resolution, only CLI and MR are able to capture this oscillating motion accurately. Oscillations can also be found for BB but the amplitude in $u_{pH}$ is too large and the value for $u_{pV}$ around which the curve oscillates is slightly changing in time. On the other hand, all PSM variants yield exponentially decaying oscillations and thus fail to capture this instability. It is worth to note that CLI is also able to reproduce the time-periodic oscillations with a resolution of $D/\Delta x=36$, whereas MR shows strong deviations from a closed curve. This motion can be analyzed in more detail by calculating the time average and fluctuation values of the different sphere velocities. These values are given in Tab. 3, where $\overline{\phi}$ denotes the average and $\phi^{\prime}$ the fluctuation part of a quantity $\phi$. Their exact definitions can be found in [48]. Table 3 also shows the frequency of the oscillation calculated with the help of a discrete Fourier transformation. It can be seen that the average of $u_{pV}$ is captured well by the MEM variants with errors well below $2\%$ for the fine resolution. In contrast to that, the errors of the PSM variants are larger and always above $2\%$. They, however, predict the average of $u_{pH}$ significantly better than the MEM counterparts with errors below $1\%$ and approximately independent of the resolution. The MEM variants tend to underestimate this quantity. Regarding $\overline{\omega}_{pH}$, all methods yield comparable results except for BB being off by a factor of three and M3B2 which is particularly close to the reference results. Having a look at the fluctuation parts, all PSM variants show large deviations as expected from their inability to capture the oscillating motion already discussed in the phase-space diagram in Fig. 8. BB shows too large fluctuations, whereas CLI and MR offer overall good predictions, especially for the finer resolution. Comparing the obtained oscillation frequencies $f$, CLI and MR can be considered again acceptable whereas BB underestimates the frequency considerably. All PSM variants predict the frequency more accurately than the MEM variants.

Summarizing, the PSM variants cannot capture the expected oscillating motion in this regime. However, they perform well in predicting the frequency and the temporal averages of the horizontal quantities, with M2B2 and M3B2 being the best choice. The MEM variants show two distinct behaviors. BB is too inaccurate and, except for $\overline{u}_{pV}$, all its values deviate significantly from the reference values. CLI and MR, on the other hand, are the methods that overall represent the expected trajectories best and additionally yield acceptable results for the rotational velocity and the frequency.

4.5 Chaotic regime

The sphere sets into a chaotic motion when the Galileo number is increased to $250$, corresponding to a particle Reynolds number of around $365$. In order to obtain reliable data for a comparison with the reference values, the simulation for this regime is run seven times. Each time, the initial sphere position is slightly perturbed at random in the horizontal plane resulting in identical physical conditions for each run. All samples combined amount to a total dimensionless simulation time of $3570$ units. These seven samples are then used to calculate average and fluctuation quantities as defined in [48]. Fig. 9 shows the phase-space diagrams of the horizontal velocity components of the seven samples. The BB variant, like the MR, shows only a few chaotic trajectories at the beginning of each sample which quickly settle into a rather time-periodic motion. In contrast to that, the motion obtained by CLI can be considered chaotic. All PSM variants exhibit an essentially circular structure in the phase-space with distinct chaotic oscillations inside this structure. More insight can be gained from the sample averages and the average amplitudes of the instantaneous fluctuations, given in Tab. 4. The only non-zero sample average, $\langle u_{pV}\rangle$, is captured reasonably well by all methods with CLI and MR showing the best results. The fluctuation of $u_{pV}$ is best reproduced by CLI, closely followed by the PSM variants. For $u_{pr}$ fluctuations, M3B2 gives the best result followed by M1B1 and CLI. Regarding the fluctuation values for the rotational velocities, CLI is again the closest to the reference value for $\omega_{pV}$. Remarkably, M3B2’s value is off by a factor of $10$ whereas the other two PSM variants perform quite well. The last quantity, the fluctuation of $\omega_{px}$, seems to be especially difficult to capture since only M3B2 comes close to the reference value. BB and MR show large deviations for almost all fluctuation quantities due to the essentially non-chaotic movement.

Apart from these statistical quantities, also the probability density functions (pdfs) of the translational and rotational sphere velocities, Figs. 10 and 11, as well as their temporal auto-correlations, Figs. 12 and 13, can be compared to the reference values. Only the cases D-CLI-36 and D-M3B2-36 are shown as BB and MR could not capture the chaotic motion and the results of the PSM variants do not differ much.

Analyzing the graphs in Fig. 10 first, CLI shows a clearly visible peak in the vertical component which creates a deviation from the assumed Gaussian shape. The plateau-like shape of $u_{pr}$ is present with a less steep descent around the value for the standard deviation $\sigma$. The PSM methods, on the other hand, show this plateau very distinctly and feature a sharp drop like the reference data. The vertical component is more pronounced towards negative fluctuations. Regarding the pdfs of the rotational velocities, the CLI results capture the Gaussian distribution quite well and also show the expected mild peak in the middle for the vertical component. The PSM methods feature a slightly narrower distribution for both, the vertical and horizontal component, and also exhibit the mild peak in $\omega_{pV}$.

Finally, the auto-correlations for both coupling techniques are compared. CLI captures the decay and dominant frequencies of the vertical translational sphere velocity very well, but shows deviations for the horizontal component with a too fast decay. As opposed to this, the M3B2 results reproduce the horizontal component accurately but contain a higher frequency in the vertical component. In Fig. 13, CLI shows the same frequency in the vertical rotational sphere velocity as the reference case with less fluctuations. Again, the horizontal component’s auto-correlation deviates from the reference in a way that the decay is slower and the frequency lower. The PSM method exhibits a very slow decay in both components but the frequency of the small visible fluctuations matches the reference ones.

Putting the findings for the chaotic regime together, CLI seems to be a better choice to capture the characteristics of vertical velocity components more accurately. The PSM variants have their strength in the prediction of the horizontal components where again M2B2 and M3B2 perform better than M1B1.

4.6 Discussion

Based on all results from the different tested regimes, substantial differences between the applied fluid-solid coupling methods can be seen. The momentum exchange method consistently yields good approximations of the sphere’s streamwise velocity component whereas the partially saturated cells method seems to have its strength in predicting the spanwise component. Apart from those general differences of the two approaches, also their variants show large variations in some cases.

The very common bounce back scheme shows good agreement with the reference data for the small Reynolds number regime but then clearly falls behind the other methods for larger Reynolds numbers. It cannot reproduce the desired oscillating and chaotic motion of cases C and D even at high resolutions. When applying the central linear interpolation scheme, which makes use of the exact surface position for achieving better accuracy, the results improve drastically. The CLI scheme can reproduce the expected characteristics of all regimes with good accuracy and yields acceptable results already at lower resolutions. This comes at the cost of a slightly increased computational effort compared to BB but the results clearly demonstrate that this extra cost is worth the effort. A further improvement by the multi-reflection scheme can only be observed at lower Reynolds numbers which is, however, not significant. The MR scheme seems to destabilize the system as the instability that triggers the oscillatory motion is already present at lower Galileo numbers. This behavior could potentially be improved by applying even finer resolutions. But here, this higher order scheme does not yet pay off at the coarser resolutions that are relevant for current computational practice. Another noticeable drawback of MR is the correction term, see Eq. (15), that requires old PDF values which have to be stored additionally. This leads to a noticeable slow down of the method. As CLI seems to work better and more robust than BB or MR for this benchmark scenario, its usage is recommended when applying the momentum exchange coupling for particulate flows.

Comparing the variants of the partially saturated cells method, the M1B1 version clearly falls behind the other two schemes in terms of accuracy. As it also does not offer any additional benefits, this version is not considered favorable for coupled simulations. On the other hand, M2B2 and M3B2 both are viable choices that yield reasonable accuracy for the tested regimes. This could mean that it is advisable to apply Eq. (22) instead of Eq. (21) in these methods This is an interesting finding also for possible future developments and improvements of the PSM formulation. All PSM methods have in common that they seemingly introduce more dissipation into the system. They shift the threshold for instabilities to occur towards higher Galileo numbers and they produce damped oscillations and less chaotic motion. Such a behavior is also reported for the immersed boundary method applied in [48] in the context of finite difference computations which could hint towards a common feature of the boundary treatment through immersed interface methods. In fact, when comparing the PSM collision operator, Eq. (20), to the BGK collision operator, Eq. (5), one could argue that the PSM method is a BGK model with modified relaxation time $\tau^{\text{PSM}}=\tfrac{\tau}{(1-B)}$ and an additional external forcing operator given by the second term of Eq. (20), i.e. $\sum_{s}B_{s}\mathcal{C}_{q,s}^{\text{solid}}$. As a consequence, $\tau^{\text{PSM}}$ increases in the vicinity of a particle and thus the apparent fluid viscosity, Eq. (6), locally increases as well.

It remains to note that all applied LBM coupling methods converge to values for the recirculation length that are smaller than the reference ones. The exact reason for this behavior is unclear since all other investigated quantities converge to the reference values. Issues with the current implementation are therefore unlikely. Possible influences originating from the applied inflow and outflow boundary conditions can also be excluded as the results remain unchanged when extending the domain by $16D$ in positive and negative $z$ direction and thus increasing the vertical domain size threefold. We note, however, that evaluating the recirculation length requires extensive interpolations of the velocity values and here differences in the implementation (also for the reference data of [48]) can easily result in substantially different values. In the present work, the $x_{pH}-z-$plane is first triangulated in the vicinity of the assumed recirculation point and then the velocity values are interpolated in a cubic trilinear fashion to obtain the contour lines using VTK 6.3.0. functionalities.

5 Conclusion

In this work, two popular, yet substantially different fluid-particle coupling approaches for the lattice Boltzmann method have been revisited and three variants of each method have been compared in detail. The applied test case of a single solid sphere settling in an ambient fluid has been found to be a suitable benchmark scenario for the evaluation of coupling methods that fully resolve the submerged particles. In contrast to other commonly applied test cases in this area, the particle Reynolds number is at the order of $100$ and the sphere is moving throughout the simulation. We find that the momentum exchange method can provide accurate predictions of the streamwise sphere velocity whereas Noble and Torczynski’s partially saturated cells method yields more accurate spanwise velocities. Even though no unique explicit recommendation for one of the methods can be given, we find that the central linear interpolation boundary scheme [12] from the family of the momentum exchange methods yields highly accurate and reliable predictions for all of the tested flow regimes. Compared to this scheme, both the bounce back and the multi-reflection scheme perform worse in the case of higher Reynolds numbers.

The partially saturated cells methods are found to seemingly add numerical diffusion which prevents the prediction of certain flow and coupling characteristics that appear at higher Reynolds numbers. Even though this affects the quality of the results, it provides a partial reason why this method is sometimes the preferred coupling method for high Reynolds number particulate flows. The extra dissipation of this method is likely to numerically stabilize the computation and thus to permit stable simulations even at lower resolutions.

To reach errors below 5% for the velocity components of the sphere, essentially all methods require a resolution of $24$ computational cells per diameter already for Reynolds numbers of $185$ and $243$. Depending on the applied method, $36$ or even $48$ cells must be used for similarly accurate predictions of the regimes with $Re=263$ and $365$. These resolutions are significantly higher than those commonly reported and used in current computational practice where often a single particle is resolved with around ten cells only. These values usually find their justification in drag force evaluations when using fixed spheres in Stokes flow. Even though such a low resolution might be enough to capture the relevant length scales in the flow, the choice for a certain resolution should also take into account the accuracy of the fluid-particle coupling at the respective flow regime. These findings agree with the arguments made by [45] that the boundary layers along the particles must be resolved sufficiently well for flows at higher Reynolds numbers. The evaluation in this article shows that a genuine direct numerical simulation of particles in inertial flows may require such an increased resolution. Future work should apply this test case also to other available fluid-particle coupling approaches not covered in this work like the immersed boundary method for LBM to obtain a complete overview of their properties.

These remarks must be seen in the light that one cannot directly relate the accuracy of a single particle system to that achieved for a system with several thousand particles, as they may be needed to simulate industrially relevant applications. In such cases, systematic evaluations must consider additional aspects like particle-particle interaction models and lubrication forces since they constitute additional numerical errors. Moreover, the evaluation of the quantities of interest when simulating such systems will often require statistical averaging. Because of such effects it is yet unclear how the findings of this article will carry over to simulations with many embedded particles. Therefore, future work must focus on developing meaningful test and benchmark cases for the direct numerical simulations of particulate flow scenarios with several particles.

Acknowledgements

The authors would like to thank Simon Bogner and Ehsan Fattahi for valuable discussions. They gratefully acknowledge the Gauss Centre for Supercomputing e.V. (www.gauss-centre.eu) for funding this project by providing computing time on the GCS Supercomputer SuperMUC at Leibniz Supercomputing Centre (LRZ, www.lrz.de).