2.1. Lattice Boltzmann method for gas–liquid two-phase flows

The LBM is a mesoscopic numerical approach based on the simplified kinetic models, which solves a specific discrete Boltzmann equation for the density distribution functions (DDFs) ${f_i}$, designed to reproduce the macroscopic Navier–Stokes equations in the low-Mach-number limit. In the standard LBM algorithm, an iteration of ‘collision and streaming’ is repeatedly executed until convergence (for steady problems) or the target time (for unsteady problems), where the collision step represents the relaxation towards the local equilibrium state due to molecular collisions and streaming represents molecular free streaming afterward. The streaming step is trivial, as it always transports the post-collision DDFs $(\,f_i^*)$ at the present lattice node $(\boldsymbol {x},t)$ to their neighbours along directions of a set of discrete velocities ${\boldsymbol {e}_i}$ within a time step $\Delta t$, i.e.

(2.1)\begin{equation} {f_i}(\boldsymbol{x} + {\boldsymbol{e}_i}\Delta t,t + \Delta t) = f_i^*(\boldsymbol{x},t). \end{equation}

To execute the next collision step, the equilibrium DDFs $f_i^{eq}$ (or their moment representations) are needed based on the updated hydrodynamic variables (density $\rho$ and velocity $\boldsymbol {u}$):

(2.2a,b)\begin{equation} {} \rho = \sum\limits_i {{f_i}} ,\quad{{ }}\rho {\boldsymbol{u}} = \sum\limits_i {{f_i}} {{\boldsymbol{e}}_i} + \Delta t{\boldsymbol{F}}/2, \end{equation}

where ${\boldsymbol {F}}$ is the total force exerted on the fluid.

In the lattice Boltzmann community, various collision schemes, such as the single-relaxation-time scheme (Qian, d'Humières & Lallemand Reference Qian, d'Humières and Lallemand1992), multiple- relaxation-time (MRT) scheme (Lallemand & Luo Reference Lallemand and Luo2000), cascaded or central-moment scheme (Geier, Greiner & Korvink Reference Geier, Greiner and Korvink2006) and the entropic-MRT scheme (Karlin, Bösch & Chikatamarla Reference Karlin, Bösch and Chikatamarla2014), can be employed to suit the problems under investigation. These schemes have been discussed in detail and integrated into a unified framework (Fei & Luo Reference Fei and Luo2017; Fei, Luo & Li Reference Fei, Luo and Li2018). In this work, the cascaded scheme is used as it possesses excellent numerical stability (Geier et al. Reference Geier, Greiner and Korvink2006), while it allows achieving non-slip boundary conditions (Fei & Luo Reference Fei and Luo2017; Fei et al. Reference Fei, Luo and Li2018). In the CLBM, the DDFs are first projected onto the central moment space, then the central moments of different orders are relaxed separately and finally the post-collision DDFs are reconstructed. To make the implementation more efficient, the collision step can be described as follows (Fei & Luo Reference Fei and Luo2017; Fei et al. Reference Fei, Luo and Li2018):

(2.3)\begin{align} f_i^*\left( {\boldsymbol{x},t} \right) &= {f_i}\left( {\boldsymbol{x},t} \right) - ({{\boldsymbol{M}}^{ - 1}}{\boldsymbol{N}^{ - 1}}{\boldsymbol{SNM}})\left[ {{f_i}\left( {\boldsymbol{x},t} \right) - f_i^{eq}\left( {\boldsymbol{x},t} \right)} \right]\nonumber\\ &\quad + \Delta t{{\boldsymbol{M}}^{ - 1}}{\boldsymbol{N}^{ - 1}}({\boldsymbol{I}} - {\boldsymbol{S}}/2){\boldsymbol{NM}}{R_i}\left( {\boldsymbol{x},t} \right), \end{align}

where $f_i^{eq}$ is the equilibrium DDF, $\boldsymbol {M}$ is a transformation matrix, $\boldsymbol {N}$ is a lower triangular shift matrix, $\boldsymbol {I}$ is the unit matrix, $\boldsymbol {S}$ is a (block) diagonal relaxation matrix and ${R_i}$ represents the forcing term in the discrete velocity space. The explicit formulations of $\boldsymbol {M}$, $\boldsymbol {N}$ and their inverses depend on the discrete velocity model and the adopted central moment set. In this work, we use the two-dimensional nine-velocity lattice (Qian et al. Reference Qian, d'Humières and Lallemand1992) with the discrete velocity set:

(2.4)\begin{equation} \left.\begin{gathered} \left| {{e_{ix}}} \right\rangle = {[0,{{ }}1,{{ }}0, - 1,{\rm{ }}0,{{ }}1, - 1, - 1,{{ }}1]^{\rm T}},\\ \left| {{e_{iy}}} \right\rangle = {[0,{{ }}0,{{ }}1,{\rm{ }}0, - 1,{{ }}1,{{ }}1, - 1, - 1]^{\rm T}}, \end{gathered}\right\} \end{equation}

where $i = 0,1, \ldots,8$, $| {\cdot } \rangle$ denotes a nine-column vector and $\textrm {T}$ denotes the transposition. To construct the CLBM, the raw and central moments of ${f_i}$ are introduced:

(2.5)\begin{equation} {k_{mn}} = \left\langle {{\,f_i}|e_{ix}^me_{iy}^n} \right\rangle {\rm{, }}~{\tilde k_{mn}} = \left\langle {{\,f_i}|{{\left( {{e_{ix}} - {u_x}} \right)}^m}{{\left( {{e_{iy}} - {u_y}} \right)}^n}} \right\rangle, \end{equation}

where m and n are positive integers. The equilibrium raw and central moments $k_{mn}^{k,eq}$ and $\tilde k_{mn}^{k,eq}$ are defined analogously by replacing ${f_i}$ with their equilibrium counterparts $f_i^{eq}$. Then, one needs to choose an appropriate moment set. The following natural raw moment set is used in the present work:

(2.6)\begin{equation} \left| {T_i^k} \right\rangle = {\left[ {k_{00}^k,k_{10}^k,k_{01}^k,k_{20}^k + k_{02}^k,k_{20}^k - k_{02}^k,k_{11}^k,k_{21}^k,k_{12}^k,k_{22}^k} \right]^{\rm T}}, \end{equation}

in which the first three moments represent the component density and momentum, and the middle three and last three stand for the viscous stress and high-order non-hydrodynamic moments, respectively. Such a choice gives very concise and sparse expressions for ${\boldsymbol {M}}$ and ${\boldsymbol {N}}$ according to the relation (Fei & Luo Reference Fei and Luo2017)

(2.7)\begin{equation} \left| {\tilde T_i^k} \right\rangle = {\boldsymbol{N}}\left| {T_i^k} \right\rangle = {\boldsymbol{NM}}\left| {\,f_i^k} \right\rangle. \end{equation}

The explicit expressions for ${\boldsymbol {M}}$ and ${\boldsymbol {N}}$ are given in Appendix A. In the implementation of (2.3), the explicit formulations of $f_i^{eq}$ and ${R_i}$ and the related matrix calculations are not needed because their central moments can be consistently defined by matching the continuous integration of the Maxwellian distribution (Fei & Luo Reference Fei and Luo2017). More specifically, for the considered moment set in (2.6), we have

(2.8) \begin{equation} \left.\begin{gathered} \left| {\tilde T_i^{eq}} \right\rangle = {\boldsymbol{NM}}\left| {\,f_i^{eq}} \right\rangle = {\left[ {\rho ,0,0,2{\rho _k}c_s^2,0,0,0,0,{\rho _k}c_s^4} \right]^{\rm T}}, \\ \left| {{C_i}} \right\rangle = {\boldsymbol{NM}}\left| {{R_i}} \right\rangle = {\left[ {0,{F_x},{F_y},0,0,0,{F_y}c_s^2,{F_x}c_s^2,0} \right]^{\rm T}}. \end{gathered}\right\} \end{equation}

To mimic the interaction between the liquid phase and gas phase, the pseudopotential- based interaction force model (Shan & Chen Reference Shan and Chen1993, Reference Shan and Chen1994), which has been successfully applied to various complex multiphase flows (Li et al. Reference Li, Luo, Kang, He, Chen and Liu2016a), is adopted:

(2.9)\begin{equation} {{\boldsymbol{F}}_{int}} =- G\psi ({\boldsymbol{x}})\sum_{i \ne 0} {{w_i}} \psi ({\boldsymbol{x}} + {{\boldsymbol{e}}_i}\Delta t){{\boldsymbol{e}}_i}, \end{equation}

where $w$ is the weight with ${w_{1 - 4}} = 1/3$ and ${w_{5 - 8}} = 1/12$, and $\psi ({\boldsymbol {x}})$ is the pseudopotential function. The square-root form pseudopotential $\psi = \sqrt {2({\,p_{EOS}} - \rho c_s^2)/G{c^2}}$ is used (Yuan & Schaefer Reference Yuan and Schaefer2006), which incorporates the non-ideal equation of state (EOS) ${p_{EOS}}$ into the system. For such a choice, the interaction strength is fixed to be $G = - 1$, and $c = 1$ and ${c_s} = \sqrt {1/3}$ are the lattice speed and sound speed, respectively. In the present work, Peng–Robinson non-ideal gas EOS is used (Peng & Robinson Reference Peng and Robinson1976):

(2.10)\begin{equation} {p_{EOS}} = \frac{{\rho \bar RT}}{{1 - b\rho }} - \frac{{a\varphi (T){\rho ^2}}}{{1 + 2b\rho - {b^2}{\rho ^2}}}, \end{equation}

where $\varphi (T) = {[1 + (0.37464 + 1.54226\varpi - 0.26992{\varpi ^2})(1 - \sqrt {T/{T_c}} )]^2}$, and we set the gas constant $\bar R = 1$ and the acentric factor $\varpi = 0.344$. The critical pressure ${p_c}$ and temperature ${T_c}$ are determined by $a = 0.4572{R^2}T_c^2/{p_c}$ and $b = 0.0778R{T_c}/{p_c}$. Unless specified, in the present work, the saturation temperature is set as ${T_{sat}} = T/{T_{cr}} = 0.775$, and $a = 1/40$ and $b = 2/21$ are chosen, which leads to a density ratio ${\rho _l}/{\rho _g} = 7.452/0.147 \approx 50$ and interface thickness $W \approx 5\Delta x$.

When such a square-root-form pseudopotential is used, the system suffers from the so-called thermodynamic inconsistency in that the liquid–vapour coexistence densities by the mechanical stability solution deviate from those of the Maxwell construction. To solve the problem, Li et al. (Reference Li, Luo and Li2012, Reference Li, Luo and Li2013) proposed to adjust the mechanical stability condition to be consistent with the Maxwell construction. Recently, such an adjustment method has been extended to the CLBM for single-component thermal multiphase systems, such as two-dimensional convective boiling (Saito et al. Reference Saito, De Rosis, Fei, Luo, Ebihara, Kaneko and Abe2021) and three-dimensional pool boiling (Fei et al. Reference Fei, Yang, Chen, Mo and Luo2020). More specifically, the central-moment forcing terms in (2.8) are slightly modified as (Saito et al. Reference Saito, De Rosis, Fei, Luo, Ebihara, Kaneko and Abe2021)

(2.11)\begin{equation} \left| {{C_i}} \right\rangle = {\left[ {0,{F_x},{F_y},\eta ,0,0,{F_y}c_s^2,{F_x}c_s^2,\eta c_s^2} \right]^{\rm T}}, \end{equation}

where $\eta = 4\sigma {| {{{\boldsymbol {F}}_{int}}} |^2}/[{\psi ^2}\Delta t(1/{s_b} - 0.5)]$ and the tuning parameter is chosen as $\sigma = 0.117$.

Using the Chapmann–Enskog analysis, the above CLBM with the pseudopotential model recovers the following macroscopic equations:

(2.12)\begin{align} \left.\begin{aligned} {\partial _t}\rho + \boldsymbol{\nabla} \boldsymbol{\cdot} (\rho {\boldsymbol{u}}) & = 0,\\ {\partial _t}(\rho {\boldsymbol{u}}) + \boldsymbol{\nabla}\boldsymbol{\cdot} (\rho {\boldsymbol{uu}}) & =- \boldsymbol{\nabla} (\rho c_s^2) + \boldsymbol{\nabla} \boldsymbol{\cdot} \left[ {\rho \nu (\boldsymbol{\nabla}{\boldsymbol{u}} + {{( \boldsymbol{\nabla}{\boldsymbol{u}})}^{\rm{T}}}) + \rho ({\nu _b} - \nu )( \boldsymbol{\nabla}\boldsymbol{\cdot}{\boldsymbol{u}}){\boldsymbol{I}}} \right] \\ & \quad + {\boldsymbol{F}} - 2{G^2}{c^4}\sigma \boldsymbol{\nabla}\boldsymbol{\cdot} ( {{{\left| {\boldsymbol{\nabla} \psi } \right|}^2}{\boldsymbol{I}}} ), \end{aligned}\right\} \end{align}

where $\nu = c_s^2\Delta t(1/{s_v} - 0.5)$ and ${\nu _b} = c_s^2\Delta t(1/{s_b} - 0.5)$ are the kinetic and bulk viscosities, respectively, and ${s_v}$ and ${s_b}$ are relaxation rates for second-order central moments in the diagonal relaxation matrix ${\boldsymbol {S}} = \textrm {diag}({s_0},{s_1},{s_1},{s_b},{s_v},{s_v},{s_3},{s_3},{s_4})$. The last term can be absorbed into the pressure tensor, leading to a thermodynamically consistent mechanical stability condition (Li et al. Reference Li, Luo and Li2013).

2.2. Discrete element method for particle dynamics

The DEM calculates the translational and angular velocities (${{\boldsymbol {U}}_a}$ and ${{\boldsymbol {\varOmega }}_a}$, respectively) of a classical multi-body system by applying Newton's second law to each grain to describe the deformation of the granular assembly:

(2.13)\begin{equation} {M_a} \cdot \frac{{{\rm d}{{{\boldsymbol{U}}_a}}}}{{{\rm d}t}} = {{{\boldsymbol{F}}_a}}{{, }}\quad {I_a} \cdot \frac{{{\rm d}{{\boldsymbol{\varOmega}}_a}}}{{{\rm d}t}} = {{\boldsymbol{T}}_a}\quad {{(}}a = 1,2,\ldots,N{{),}} \end{equation}

where N is the number of grains in the system and ${M_a}$ and ${I_a}$ are the mass and the moment of inertia of the a-th grain, respectively. The method consists of calculating the total force ${{\boldsymbol {F}}_a}$ and total torque ${{\boldsymbol {T}}_a}$ acting on each particle and the numerical solution of the ordinary differential equation (2.13) (Cundall & Strack Reference Cundall and Strack1979). In the present work, we consider each grain as a two-dimensional circular particle (disc) with radius ${r_a}$ and all grains are made of the same material with density ${\rho _p}$; therefore we have ${M_a} = {\rho _p}{\rm \pi} r_a^2$ and ${I_a} = {\rho _p}{\rm \pi} r_a^4/2$. The angular velocity ${{\boldsymbol {\varOmega }}_a} = \omega {\boldsymbol {k}}$ has only one non-zero component $\omega$ (${\boldsymbol {k}}$ is the unit vector perpendicular to the paper).

The total force in (2.13) is composed of the contact interaction forces ${\boldsymbol {F}}_a^b$ (between particles a and b), hydrodynamic force ${\boldsymbol {F}}_a^h$, capillary force ${\boldsymbol {F}}_a^{cap}$ and body force ${\boldsymbol {F}}_a^v$. Each force leads to a torque contribution, except for the body force that always acts along the centre of mass. Thus,

(2.14a,b)\begin{equation} {\boldsymbol{F}}_a = \sum\limits_b {\boldsymbol{F}}_a^b + {\boldsymbol{F}}_a^h + {\boldsymbol{F}}_a^{cap} + {\boldsymbol{F}}_a^v,\quad {\boldsymbol{T}}_a = \sum\limits_b {\boldsymbol{T}}_a^b + {\boldsymbol{T}}_a^h + {\boldsymbol{T}}_a^{cap}. \end{equation}

The hydrodynamic force, capillary force and their torques come from the interaction between the fluid and solid particles, and are discussed in the next subsection. The body force is most often the gravity force, ${\boldsymbol {F}}_a^v = - {M_a}g\,{\boldsymbol {j}}$, where g is the magnitude of gravitational acceleration and ${\boldsymbol {j}}$ is the unit vector in the vertical direction.

In the DEM, the solid particles are not deformable, but they can overlap in the sense they are considered to be rigid while their contact is assumed to be soft. Considering two particles a and b in contact, as shown in figure 1, the overlap zone at the contact is limited to be very small by introducing a repulsive normal force ${F_n}$. Related to the relative tangential motion, there is a tangential force in addition to the normal force ${F_t}$. Therefore, the contact interaction force can be decomposed into two components:

(2.15)\begin{equation} {\boldsymbol{F}}_a^b = {F_n}{\boldsymbol{n}} + {F_t}{\boldsymbol{t}}, \end{equation}

where ${\boldsymbol {n}}$ and ${\boldsymbol {t}}$ are unit vectors along the normal and tangential directions, respectively. The particle contact torque is

(2.16)\begin{equation} {\boldsymbol{T}}_a^b = ({{\boldsymbol{X}}_c} - {{\boldsymbol{X}}_a}) \times {F_t}{\boldsymbol{t}}, \end{equation}

where ${{\boldsymbol {X}}_a}$ and ${{\boldsymbol {X}}_b}$ are the vectors fixing the particle centre and contact point, respectively. Following Luding et al. (Reference Luding, Clément, Blumen, Rajchenbach and Duran1994), we consider the contact as a linear spring with damping. The repulsive normal force depends linearly on the mutual compression distance ${\delta _n}$, which is defined as

(2.17)\begin{equation} {\delta _n} = \left| {{{\boldsymbol{X}}_a} - {{\boldsymbol{X}}_b}} \right| - {r_a} - {r_b}. \end{equation}

The normal force appears when ${\delta _n}<0$ and is controlled by the stiffness of the grain particle ${k_n}$. Also, another damping term that opposes the relative velocity of two contact particles is added to consider energy dissipation. The normal force can be explicitly written as (Luding et al. Reference Luding, Clément, Blumen, Rajchenbach and Duran1994)

(2.18)\begin{equation} {F_n} = \left\{ \begin{array}{@{}ll} 0, & {{ }}{\delta _n} \geq 0\\ - {k_n}{\delta _n} - {\gamma _n}\dfrac{{{\rm d}{\delta _n}}}{{{\rm d}t}}, & {{ }}{\delta _n} < 0,{{ }} \end{array} \right. \end{equation}

where ${\gamma _n}$ is the normal damping constant and the relative normal velocity is calculated by $\textrm {d}{\delta _n}/\textrm {d}t = ({{\boldsymbol {U}}_a} - {{\boldsymbol {U}}_b}) \boldsymbol {\cdot } {\boldsymbol {n}}$. To avoid the appearance of an attractive force, one can simply set ${F_n} = 0$ when (2.18) produces ${F_n} > 0$ (Pöschel & Schwager Reference Pöschel and Schwager2005). The tangential force is defined similarly, and is assumed to obey Coulomb's friction law, i.e.

(2.19)\begin{equation} {F_t} = \left\{ \begin{array}{@{}ll} {{\rm sgn}} ({{\hat{F}}_t}) \cdot {\mu _f}{F_n}, & {{ }}\left| {{{\hat{F}}_t}} \right| \geq {\mu _f}{F_n}{{ }}\\ {{\hat{F}}_t}, & {{ }}\left| {{{\hat{F}}_t}} \right| < {\mu _f}{F_n}{{ }} \end{array} \right.\quad {{}}{\hat F_t} =- {k_t}{\delta _t} - {\gamma _t}\frac{{{\rm d}{\delta _t}}}{{{\rm d}t}}, \end{equation}

where ${k_t}$ is the stiffness of the tangential spring, ${\gamma _t}$ is the damping constant, ${\mu _f}$ is the friction coefficient and the relative tangential velocity is calculated by $\textrm {d}{\delta _t}/\textrm {d}t = ({U_a} - {{\boldsymbol {U}}_b}) \boldsymbol {\cdot } {\boldsymbol {t}} - ({r_a}{\varOmega _a} + {r_b}{\varOmega _b}) \times {\boldsymbol {n}}$. To save memory, we only consider a particle interacting with those particles within a cut-off radius, which is achieved using the Verlet list algorithm (Verlet Reference Verlet1967). By default, the DEM parameters in (2.18) and (2.19) are set according to the physical properties of granite sand (Soundararajan Reference Soundararajan2015).

Figure 1. Illustration of the interaction between particles a and b.

In the time evolution of the DEM, a suitable numerical scheme is needed to integrate (2.13) twice, i.e. to obtain the particle linear and angular velocities as well as the translational and angular positions. Here we use the Verlet integration algorithm (Verlet Reference Verlet1967) since it has second-order accuracy, is fast to implement and requires less storage memory. The Verlet scheme first evaluates the velocities at the half-integer time steps and computes the new positions:

(2.20)\begin{equation} \left.\begin{gathered} {{\boldsymbol{U}}_a}\left(t + \dfrac{{\Delta {t_{DEM}}}}{2}\right) = {{\boldsymbol{U}}_a}(t) + \frac{{\Delta {t_{DEM}}}}{2} \cdot \frac{{{{\boldsymbol{F}}_a}(t)}}{{{M_a}}},\\ {{}}{\omega _a}\left(t + \frac{{\Delta {t_{DEM}}}}{2}\right){\boldsymbol{k}} = {\omega _a}(t){\boldsymbol{k}} + \frac{{\Delta {t_{DEM}}}}{2} \cdot \frac{{{{\boldsymbol{T}}_a}(t)}}{{{I_a}}},\\ {{\boldsymbol{X}}_a}(t + \Delta {t_{DEM}}) = {{\boldsymbol{X}}_a}(t) + \Delta {t_{DEM}}{{\boldsymbol{U}}_a}\left(t + \frac{{\Delta {t_{DEM}}}}{2}\right),\\ {{ }}{\varphi _a}(t + \Delta {t_{DEM}}) = {\varphi _a}(t) + \Delta {t_{DEM}}\omega \left(t + \frac{{\Delta {t_{DEM}}}}{2}\right), \end{gathered}\right\} \end{equation}

where ${\varphi _a}$ is particle angular position and $\Delta {t_{DEM}}$ is the time step in the DEM, which is not necessarily equal to $\Delta t$ in the LBM, as is discussed later. Using the new position information, the particle velocities are updated as

(2.21)\begin{equation} \left.\begin{gathered} {{\boldsymbol{U}}_a}(t + \Delta {t_{DEM}}) = {{\boldsymbol{U}}_a} \left(t + \frac{{\Delta {t_{DEM}}}}{2}\right) + \frac{{\Delta {t_{DEM}}}}{2} \cdot \frac{{{{\boldsymbol{F}}_a}(t + \Delta {t_{DEM}})}}{{{M_a}}},{{ }}\\ {{\boldsymbol{\varOmega}}_a}(t + \Delta {t_{DEM}}) = {{\boldsymbol{\varOmega}}_a} \left(t + \frac{{\Delta {t_{DEM}}}}{2}\right) + \frac{{\Delta {t_{DEM}}}}{2} \cdot \frac{{{{\boldsymbol{T}}_a}(t + \Delta {t_{DEM}})}}{{{I_a}}}. \end{gathered}\right\} \end{equation}

2.3. Fluid–solid interactions

The hydrodynamic force is calculated based on the momentum-exchange method (Ladd Reference Ladd1994a). The basic consideration is as sketched in figure 2: at a boundary point ${{\boldsymbol {x}}_b}$, some post-collision DDFs $f_i^*({{\boldsymbol {x}}_b},t)$ stream into the solid particle and each of them contributes a momentum to the particle, i.e. $f_i^*({{\boldsymbol {x}}_b},t){{\boldsymbol {e}}_i}$, while at the solid point ${{\boldsymbol {x}}_s}$, the opposite post-collision DDFs $f_{\boldsymbol {i}}^*({{\boldsymbol {x}}_s},t)$ stream out of the solid particle and each of them contributes a momentum decrement $f_{\boldsymbol {i}}^*({{\boldsymbol {x}}_s},t){{\boldsymbol {e}}_{\boldsymbol {i}}}$. Wen et al. (Reference Wen, Zhang, Tu, Wang and Fang2014) proposed that the relative velocity (discrete velocity displaced by the wall velocity ${\boldsymbol {u}}_w$) should be used in the momentum computation and proposed the following Galilean invariant momentum-exchange method:

(2.22)\begin{equation} \left.\begin{gathered} {\boldsymbol{F}}_a^h =\sum {\left[ {({{\boldsymbol{e}}_i} - {{\boldsymbol{u}}_w})f_i^*({{\boldsymbol{x}}_b},t) - ({{\boldsymbol{e}}_{\boldsymbol{i}}} - {{\boldsymbol{u}}_w})f_{\boldsymbol{i}}^*({{\boldsymbol{x}}_s},t)} \right]},\\ {{}}{\boldsymbol{T}}_a^h = \sum {\left[ {({{\boldsymbol{e}}_i} - {{\boldsymbol{u}}_w})f_i^*({{\boldsymbol{u}}_b},t) - ({{\boldsymbol{e}}_{\boldsymbol{i}}} - {{\boldsymbol{u}}_w})f_{\boldsymbol{i}}^*({{\boldsymbol{x}}_s},t)} \right]} \times ({{\boldsymbol{x}}_w} - {{\boldsymbol{X}}_a}), \end{gathered}\right\} \end{equation}

where ${{\boldsymbol {X}}_a}$ is the particle centre and the summation runs over all the fluid–solid links. Similar to (2.5), the Galilean invariant momentum-exchange method in (2.22) is based on central moments of the post-collision distributions, which is therefore very consistent with the central-moment-based CLBM used in § 2.1. The velocity of each wall point ${{\boldsymbol {x}}_w}$ is determined by

(2.23)\begin{equation} {{\boldsymbol{u}}_w} = {{\boldsymbol{U}}_a} + \omega {\boldsymbol{k}} \times ({{\boldsymbol{x}}_w} - {{\boldsymbol{X}}_a}). \end{equation}

Figure 2. Illustration of fluid–solid interactions and the curved boundary scheme. The solid line connected by wall points represents the boundary of particle a.

The capillary force calculation depends on the contact angle scheme. In this work, the improved virtual-density contact angle scheme by Li, Yu & Luo (Reference Li, Yu and Luo2019) is used due to its simplicity and locality. The basic idea is that the solid points in the neighbourhood of the liquid are allocated virtual densities by

(2.24)\begin{equation} {\rho _s} = \left\{\begin{array}{@{}ll} \phi {\rho _{ave}}({{\boldsymbol{x}}_s}), & {{ }}\theta \leq {90^ \circ },\\ {\rho _{ave}}({{\boldsymbol{x}}_s}) - \delta \rho , & {{}}\theta > {90^ \circ }, \end{array} \right. \end{equation}

where ${\rho _{ave}}({{\boldsymbol {x}}_s})$ is calculated based on the weighted average of the density at the neighbouring fluid points (Li et al. Reference Li, Yu and Luo2019). The parameters $\phi$ and $\delta \rho$ are used to tune the contact angle $\theta$. Then, the pseudopotential at the solid point is substituted into (2.9) to calculate the interaction force at the boundary point. The capillary force and torque are obtained by adding up the adhesive force contribution along all the fluid–solid links:

(2.25)\begin{equation} \left.\begin{gathered} {\boldsymbol{F}}_a^c = \sum_{{{\boldsymbol{x}}_b}} {G\psi ({{\boldsymbol{x}}_b})\sum_{{{\boldsymbol{x}}_b} + {{\boldsymbol{e}}_i}\Delta t = {{\boldsymbol{x}}_s}} {{w_i}} \psi ({{\boldsymbol{x}}_b} + {{\boldsymbol{e}}_i} \Delta t){{\boldsymbol{e}}_i}} ,\\ {\boldsymbol{T}}_a^c = \sum_{{{\boldsymbol{x}}_b}} ({{\boldsymbol{x}}_w} - {{\boldsymbol{X}}_a}) \times \left[G\psi ({{\boldsymbol{x}}_b})\sum_{{{\boldsymbol{x}}_b} + {{\boldsymbol{e}}_i}\Delta t = {{\boldsymbol{x}}_s}} {{w_i}} \psi ({{\boldsymbol{x}}_b} + {{\boldsymbol{e}}_i}\Delta t){{\boldsymbol{e}}_i}\right]. \end{gathered}\right\} \end{equation}

2.4. Curved boundary treatment

In (2.22), the opposite post-collision DDFs $f_{\boldsymbol {i}}^*({{\boldsymbol {x}}_s},t)$ at the solid point node are needed. However, inside the solid particle, the lattice Boltzmann solver is not implemented. This requires boundary condition treatment: how does one obtain $f_{\boldsymbol {i}}^*({{\boldsymbol {x}}_s},t)$ or ${f_{\boldsymbol {i}}}({{\boldsymbol {x}}_b},t + \Delta t)$ after the streaming step? A simple and straightforward strategy is the half-way bounce-back scheme, i.e. ${f_{\boldsymbol {i}}}({{\boldsymbol {x}}_b},t + \Delta t) = f_i^*({{\boldsymbol {x}}_b},t)$. Since it is first-order-accurate, more lattice nodes are needed to resolve one solid particle using the half-way bounce-back scheme, which limits the total number of particles in our system. Many second- and higher-order curved boundary schemes have also been proposed in the literature (Yu et al. Reference Yu, Mei, Luo and Shyy2003a), among which the method originally proposed by Filippova & Hänel (Reference Filippova and Hänel1998) is used in the present work.

As sketched in figure 2, the physical boundary intersects the fluid–solid link between the boundary node ${{\boldsymbol {x}}_b}$ and the solid node ${{\boldsymbol {x}}_s}$ at the wall node ${{\boldsymbol {x}}_w}$. The fraction of the intersected link in the fluid region is defined as $\delta = | {{{\boldsymbol {x}}_b} - {{\boldsymbol {x}}_w}} |/| {{{\boldsymbol {x}}_b} - {{\boldsymbol {x}}_s}} |$. To obtain a second-order scheme, Filippova & Hanel proposed to obtain $f_{\boldsymbol {i}}^*({{\boldsymbol {x}}_s},t)$ using the following linear interpolation:

(2.26)\begin{equation} f_{\boldsymbol{i}}^*({x_s},t) = (1 - \chi )\,f_{_i}^ * ({x_b},t) + \chi f_i^{eq}({x_s},t) + 2{\omega _i}\rho ({x_b})\frac{{{e_{\boldsymbol{i}}} \cdot {u_w}}}{{c_s^2}}, \end{equation}

where $\chi$ is a weighting factor. The main feature of the scheme by Filippova & Hanel lies in that it introduces a fictitious equilibrium distribution function given by

(2.27)\begin{equation} f_{_i}^{eq}({x_s},t) = {\omega _i}\rho ({x_b},t)\left[ {1 + \frac{{{e_i} \cdot {u_s}}}{{c_s^2}} + \frac{{{{({e_i} \cdot {u_b})}^2}}}{{2c_s^4}} - \frac{{{u_b} \cdot {u_b}}}{{2c_s^2}}} \right], \end{equation}

where ${{\boldsymbol {u}}_s}$ is the fictitious velocity at ${{\boldsymbol {x}}_s}$, and the weights ${\omega _i}$ are defined as ${\omega _0} = 4/9$, ${\omega _{1 - 4}} = 1/9$ and ${\omega _{5 - 8}} = 1/36$. The weighting factor $\chi$ depends on $\delta$ and ${{\boldsymbol {u}}_s}$. In this work, we use the original choice by Filippova & Hanel for $\delta \geq 1/2$ and the improved scheme by Mei et al. (Reference Mei, Luo and Shyy1999) for $\delta < 1/2$, namely

(2.28)\begin{equation} \left.\begin{gathered} \delta > 1/2:{{ }}\quad{{\boldsymbol{u}}_s} = \frac{{(\delta - 1)}}{\delta }{{\boldsymbol{u}}_f} + \frac{1}{\delta }{{\boldsymbol{u}}_w},\quad {{ }}\chi = (2\delta - 1){s_v},\\ \delta \leq 1/2:{{ }}\quad{{\boldsymbol{u}}_s} = {{\boldsymbol{u}}_f},\quad {{ }}\chi = \frac{{(2\delta - 1){s_v}}}{{(1 - {s_v})}}, \end{gathered}\right\} \end{equation}

where ${{\boldsymbol {u}}_f}$ is the fluid velocity at the node next to the boundary node ${{\boldsymbol {x}}_f} = {{\boldsymbol {x}}_b} - {{\boldsymbol {e}}_i}\Delta t$.

The above scheme performs well for the fixed curved boundary in single-phase flows. However, in multiphase systems, the scheme leads to mass leakage due to the unbalance of outgoing and incoming mass at the boundary nodes. As a result, for a static droplet resting on a disc or sphere, the droplet shrinks or expands spontaneously (Yu et al. Reference Yu, Li and Wen2020; Yao et al. Reference Yao, Liu, Zhong and Wen2022). To overcome mass leakage, Yu et al. (Reference Yu, Li and Wen2020) proposed a simple strategy by adding a compensation term in the rest DDFs:

(2.29)\begin{equation} {f_0}({{\boldsymbol{x}}_b},t + \Delta t) = f_0^ * ({{\boldsymbol{x}}_b},t) + \left[ {\sum_{outgoing} {f_i^ * ({{\boldsymbol{x}}_b},t)} - \sum_{ingoing} {{f_{\boldsymbol{i}}}({{\boldsymbol{x}}_b},t + \Delta t)} } \right]. \end{equation}

Yao et al. (Reference Yao, Liu, Zhong and Wen2022) proposed another scheme which also uses the above compensation term but slightly modified, ${f_{\boldsymbol {i}}}({{\boldsymbol {x}}_b},t + \Delta t)$, by considering the non-ideal force effect. It may be noted that both schemes are developed for a fixed curved boundary, in the sense that the particle cannot move. For a moving curved boundary, the third term has been added in (2.26) to include the boundary velocity effect. Based on simple tests, we find (2.29) still leads to sensible mass leakage for moving boundary problems, and we propose to use the following alternative scheme by removing the velocity correction term in ${f_{\boldsymbol {i}}}({{\boldsymbol {x}}_b},t + \Delta t)$, i.e.

(2.30)\begin{align} {f_0}({{\boldsymbol{x}}_b},t + \Delta t) &= f_0^ * ({{\boldsymbol{x}}_b},t) \nonumber\\ &\quad + \left( {\sum_{outgoing} {f_i^ * ({{\boldsymbol{x}}_b},t)} - \sum_{ingoing} {\left[\,{{f_{\boldsymbol{i}}}({{\boldsymbol{x}}_b},t + \Delta t) - 2{\omega _i}\rho ({{\boldsymbol{x}}_b})\frac{{{{\boldsymbol{e}}_{\boldsymbol{i}}} \boldsymbol{\cdot} {{\boldsymbol{u}}_w}}}{{c_s^2}}} \right]} } \right). \end{align}

2.5. Refilling scheme and effective hydraulic particle radius

For our considered problems, some nodes may emerge from the solid region to the fluid domain due to the movement of particles in the time interval t to $t + \Delta t$, as sketched in figure 3. To initialize the DDFs at these points, we use the following non-equilibrium refilling scheme:

(2.31)\begin{equation} {f_i}({{\boldsymbol{x}}_b},t + \Delta t) = f_i^{eq}\left[ {{\rho _{ave}}({{\boldsymbol{x}}_b}),{\boldsymbol{u}}({{\boldsymbol{x}}_b})} \right] + f_i^{neq}({{\boldsymbol{x}}_b},t + \Delta t), \end{equation}

where ${\rho _{ave}}({{\boldsymbol {x}}_b})$ is calculated based on the weighted average of the density at the neighbouring fluid nodes, and ${\boldsymbol {u}}({{\boldsymbol {x}}_b})$ is calculated according to the particle velocity:

(2.32)\begin{equation} {\boldsymbol{u}}({{\boldsymbol{x}}_b}) = {{\boldsymbol{U}}_a} + {\omega _a}({{\boldsymbol{x}}_b} - {{\boldsymbol{X}}_a}). \end{equation}

Previously, Caiazzo (Reference Caiazzo2008) simply copied the non-equilibrium part from a neighbour of the new fluid node. However, such a neighbour of the present new fluid node along a certain discrete velocity direction could as well be a new fluid node or even a solid node, which, for either, the non-equilibrium part is unknown. To avoid ambiguity and improve stability, we use a weighted average method to approximate the non-equilibrium part:

(2.33)\begin{equation} f_i^{neq}({{\boldsymbol{x}}_b},t + \Delta t) \approx \frac{{\displaystyle\sum {{\omega _i}\left[ {{f_i}({\boldsymbol{x}} + {{\boldsymbol{e}}_i}\Delta t,t + \Delta t) - f_i^{eq}({\boldsymbol{x}} + {{\boldsymbol{e}}_i}\Delta t,t + \Delta t)} \right]} }}{{\displaystyle\sum {{\omega _i}} }}, \end{equation}

where the summation runs over all the old neighbouring fluid nodes. Taking ${{\boldsymbol {x}}_b}$ in figure 3 as an example, the summation includes ${{\boldsymbol {x}}_c}$ , ${{\boldsymbol {x}}_d}$ and ${{\boldsymbol {x}}_e}$.

Figure 3. Illustration of refilling scheme.

As discussed by Jiang et al. (Reference Jiang, Liu, Chen and Tsuji2022), two-dimensional simulations may encounter the problem that there are no connected paths for fluid to flow through densely packed sediments, while in three dimensions pore spaces may actually be interconnected. In addition, when two particles are in contact with each other, it brings difficulties in treating curved boundary conditions using the scheme discussed in above. To solve this problem, we employ a method similar to that proposed by Boutt, Cook & Williams (Reference Boutt, Cook and Williams2011), by introducing an effective hydraulic radius ${r_{a,h}}$ for each particle a in LBM simulations. By default, the effective hydraulic radius in the LBM is $2.5\Delta x$ smaller than the particle radius in the DEM, i.e. ${r_{a,h}} = {r_a} - 2.5\Delta x$.

2.6. Unit conversion and time-step matching

Usually, lattice units are used in LBM simulations while physical units are used in DEM simulations. In the coupled LBM–DEM model, information exchange is needed at the interface of the LBM and DEM solvers. For example, we need to convert the lattice unit hydrodynamic force and capillary forces calculated based on DDFs and pseudopotential functions into physics units for solving particle dynamics. We also need to convert the physical unit particle locations and velocities into lattice units to deal with the boundary conditions and refilling in the LBM solver. To this aim, unit conversion is needed. The variables in lattice units can be converted into physical units based on the reference variables, and we consider the following three primary variables, namely reference length $\Delta {{x}_r}$, reference time $\Delta {t_r}$ and reference mass ${m_r}$, i.e.

(2.34a–c)\begin{equation} \Delta {x_r} = \frac{{{L_p}}}{{{L_l}}},\quad{{ }}\Delta {t_r} = \frac{{{\nu _l}}}{{{\nu _p}}}\Delta x_r^2,\quad{{}}{m_r} = \frac{{{\rho_l}}}{{{\rho _p}}}\Delta x_r^2, \end{equation}

where L is the system size and the indexes p and l indicate a variable in physical and lattice units, respectively. Other reference variables, such as reference acceleration, reference force, reference torque, reference translation velocity and reference angular velocity, are derived from the primary variables:

(2.35a–e)\begin{equation} {a_r} = \frac{{\Delta {x_r}}}{{\Delta t_r^2}},\quad{{ }}{F_r} = {m_r}\frac{{\Delta {x_r}}}{{\Delta t_r^2}},\quad{{ }}{T_r} = {m_r}\frac{{\Delta x_r^2}}{{\Delta t_r^2}},\quad{{}}{V_r} = \frac{{\Delta {x_r}}}{{\Delta {t_r}}},\quad{{}}{w_r} = \frac{1}{{\Delta {t_r}}}. \end{equation}

As noted above, the time steps used in the LBM and in the DEM do not need to be identical. Considering both in physical units, the LBM time step is determined by fluid viscosity ${v_{phys}}$ and the mesh size $\Delta {x_{phys}}$, i.e. $\Delta {t_{LBM}} = (1/{s_v} - 0.5)\Delta x_{phys}^2/3{v_{phys}}$, while the DEM time step $\Delta {t_{DEM}}$ should be smaller than a critical oscillation length scale by

(2.36)\begin{equation} \Delta t_{DEM}^{cr} = \lambda {\rm \pi}{r_{min }}\sqrt {{\rm \pi} {\rho _p}/{k_g}}, \end{equation}

where ${r_{min}}$ is the smallest particle radius and the step factor $\lambda$ is chosen to be around 0.1 to ensure both stability and accuracy (Soundararajan Reference Soundararajan2015). Usually, $\Delta {t_{DEM}}$ is smaller than $\Delta {t_{LBM}}$; therefore a subcycling time integration with ${N_{sub}}$ steps in the DEM is required to match one LBM step. To this aim, the DEM time step is chosen as

(2.37a,b)\begin{equation} \Delta {t_{DEM}} = {{\Delta {t_{LBM}}} / {{N_{sub}},\quad {{ }}{N_{sub}} = \left\lceil {\Delta t/\Delta t_{DEM}^{cr}} \right\rceil + 1,}} \end{equation}

where $\lceil {\cdot } \rceil$ is an integer round-off operator.

3.1. Single-disc sedimentation in single-phase fluid

The considered problem is sketched in figure 4. A disc with diameter D is initially released away from the channel centre, which is filled with static single-phase fluid. Due to its higher density than that of the fluid, the disc sinks and rotates under gravity and the hydrodynamic force. The disc accelerates in the early stage and finally reaches a steady state with a constant descending velocity along the centreline since the drag force and buoyancy exactly offset gravity. To compare with the arbitrary Lagrangian–Eulerian technique (ALE) simulation (Hu, Patankar & Zhu Reference Hu, Patankar and Zhu2001) and previous LBM simulations by Wen et al. (Reference Wen, Li, Zhang and Fang2012, Reference Wen, Zhang, Tu, Wang and Fang2014), the channel width is set as $4 \times {10^{ - 3}}$ m and the disc diameter is ${10^{ - 3}}$ m. The fluid density and kinematic viscosity are ${10^{3}}\ \textrm {kg}\ \textrm {m}^{-3}$ and ${10^{ -6}}\ \textrm {m}^2\ \textrm {s}^{-1}$. The disc is released at $7.6 \times {10^{ - 4}}$ m away from the left wall, and the gravitational acceleration is $g = - 9.8\ \textrm {m} \textrm {s}^{-2}$. The channel width is resolved by $120{{}}\Delta x$, and the length-to-width ratio is 12.5. The terminal vertical velocity ${v_p}$ depends on the particle density, and we consider two cases with ${\rho _p} = 1.01 \times {10^3}$ and $1.03 \times {10^3}\ \textrm {kg}\ \textrm {m}^{-3}$, leading to a terminal Reynolds number (defined as ${Re} = D{v_p}/\nu$) $Re=3.23$ and 8.33, respectively. In the LBM solver, the kinetic viscosity is set as 0.033, leading to ${s_v} = 1.667$. The time steps in LBM and DEM are $\Delta {t_{LBM}} = 3.704 \times {10^{ - 5}}$ s and $\Delta {t_{DEM}} = 6.173 \times {10^{ - 6}}$ s, respectively, with ${N_{sub}} = 6$.

Figure 4. Illustration of single-disc sedimentation problem.

The time-dependent trajectories, angular velocities, horizontal velocities and vertical velocities for the two cases are shown in figures 5 and 6 together with a comparison with the results by ALE (Hu et al. Reference Hu, Patankar and Zhu2001) and Wen et al. (Reference Wen, Li, Zhang and Fang2012, Reference Wen, Zhang, Tu, Wang and Fang2014). In terms of trajectories, angular velocities and horizontal velocities, the present simulations are in excellent agreement with the two reference solutions. The present model gives slightly lower terminal vertical velocities than Wen et al., deviates slightly from ALE for case $Re=3.23$ (figure 5d) but agrees better with ALE for case $Re=8.33$ (figure 6d). The difference between the present model and Wen et al. lies in that we use the Filippova & Hanel scheme to deal with curved boundary conditions and the non-equilibrium scheme to initialize the newborn fluid nodes, while the quadratic interpolation scheme and second-order extrapolation scheme have been used by Wen et al. As we have discussed above, the interpolation/extrapolation schemes work well for suspended particles in single-phase fluids, but suffer from numerical instabilities in two-phase fluids. Nevertheless, the differences between our method and Wen et al. are negligible, which confirms the accuracy of the curved boundary scheme and refilling scheme used in the present model.

Figure 5. Single-disc sedimentation at $Re=3.23$. Time-dependent (a) trajectories, (b) angular velocities, (c) horizontal velocities and (d) vertical velocities.

Figure 6. Single-disc sedimentation at $Re=8.33$. Time-dependent (a) trajectories, (b) angular velocities, (c) horizontal velocities and (d) vertical velocities.

In the above simulations, the resolution is chosen the same as that of Wen et al. (Reference Wen, Zhang, Tu, Wang and Fang2014) by setting $D = 30\Delta x$. Here, we want to test the dependence of the results on the grid size by considering another two settings, i.e. $D = 20\Delta x$ and $D = 25\Delta x$. For the trajectories, horizontal and vertical velocities, the differences among the results based on different resolutions are almost negligible, as shown in figures 6(a), 6(c) and 6(d), respectively. For the angular velocities, the differences are also very small, but the results based on coarser meshes (especially for $D = 20\Delta x$) present some oscillations, as seen in the inset of figure 6(b). We have also conducted more simulations with different resolutions at different $Re$; the trend is similar and is not shown here. Generally, we see that accurate and smooth results can be obtained at the condition $D \geq 25\Delta x$. Such a criterion will be taken into account in more complex cases in the following.

3.2. Single particle floating on the liquid–gas interface

We now simulate a particle with a diameter D floating at the liquid–gas interface, as sketched in figure 7. The system has a size of $2D \times 2D$ with non-slip boundary conditions on all the walls and is divided into two equal parts, filled with gas phase and liquid phase at the top and bottom, respectively. The particle is initially located in the centre of the system. To make the particle's steady location solely adjusted by the contact angle, gravity is neglected in this case (Zhang et al. Reference Zhang, Liu and Zhang2020). Due to the capillary force, the particle will move adaptively to a steady location, where the contact angle equals the equilibrium contact angle. In the simulations, the particle diameter is 5 mm, resolved by 100 $\Delta x$. The particle density, liquid density and liquid viscosity are chosen as ${\rho _p} = 2000\ \textrm {kg}\ \textrm {m}^{-3}$, ${\rho _l} = 1000\ \textrm {kg}\ \textrm {m}^{-3}$ and ${\nu _l} = {10^{ - 4}}\ \textrm {m}^2\ \textrm {s}^{-1}$, respectively. In the LBM, the viscosity is set as ${\nu _l} = 0.05$. For this problem, the time steps in LBM and DEM are identical, $\Delta {t_{LBM}} = \Delta {t_{DEM}} = 1.25 \times {10^{ - 6}}\ \textrm {s}$, since $\Delta {t_{LBM}} < \Delta {t_{cr}}$ and the subcycling step ${N_{sub}}$ in (2.37a,b) is set to be unity.

Figure 7. Illustration of single particle floating on the liquid–gas interface.

We consider different cases by changing $\phi$ and $\delta \rho$ in (2.24) to achieve various contact angles. Snapshots for three typical cases are shown in figure 8. According to (2.24), the case with $\phi = 1.25$ and $\delta \rho = 0.0$ leads to a hydrophilic wetting condition; therefore the particle moves downward with some transient oscillation and finally reaches a steady position (seen in figure 8a–d), with a contact angle measured using ImageJ Fiji software of ${\theta _1} = {56.2^ \circ }$. For the case in figure 8(e–h) with $\phi = 1.0$ and $\delta \rho = 0.0$, a neutral wetting condition is expected and the measured steady contact angle is ${\theta _1} = {94.2^ \circ }$. Although not exactly at ${90^ \circ }$, it is within the acceptable range compared with the results reported in Li et al. (Reference Li, Yu and Luo2019). The case with $\phi = 1.0$ and $\delta \rho = 0.3$ corresponds to a hydrophobic wetting condition; accordingly the particle moves upward and finally reaches the steady position with ${\theta _1} = {120.8^ \circ }$, as seen in figure 8(i–l). Among the three cases, the transient time is the smallest ($\sim$0.75 s) for the case in figure 8(e–h) since its initial location is closest to the steady state and capillary oscillation is barely present. Such an argument is further confirmed by the transient vertical locations for different cases shown in figure 9(a), where for the cases with larger deviation from the neutral case, either hydrophilic or hydrophobic, the oscillation is more significant. For a consistent contact angle scheme, one may expect that the contact angle on a moving boundary agrees with the contact angle measured on a fixed boundary. To this aim, a static droplet sitting on a fixed particle is also simulated, and the measured contact angle ${\theta _2}$ is compared with ${\theta _1}$ for different settings in figure 9(b). It is clearly seen that the two contact angles agree well with each other, which confirms the consistency of the contact angle implementation in the present model.

Figure 8. Evolution of a floating solid particle (in white) at the liquid–gas interface (liquid phase: in red; gas phase: in blue): (a–d) $\phi = 1.25$, $\delta \rho = 0.0$; (e–h) $\phi = 1.0$, $\delta \rho = 0.0$; (i–l) $\phi = 1.0$, $\delta \rho = 0.3$.

Figure 9. (a) Transient vertical location of a solid particle floating at the liquid–gas interface. (b) Comparison between the contact angles measured from a floating particle (${\theta _1}$: top left inset) and a static droplet sitting on a fixed particle (${\theta _2}$: bottom right inset).

3.3. Sinking of a horizontal cylinder

In order to further test the performance of our proposed LBM–DEM for the gas–liquid–solid interaction problems in the presence of moving contact lines, the sinking of a horizontal cylinder is investigated. The problem has been studied experimentally by Vella, Lee & Kim (Reference Vella, Lee and Kim2006) and numerically by He et al. (Reference He, Li, Huang, Hu and Wang2019) and Zhang et al. (Reference Zhang, Liu and Zhang2020). As sketched in figure 10, a circular cylinder with a diameter $D = 5\ \textrm {mm}$ (long enough in the axial direction) lies horizontally at the liquid–gas interface. The density of the cylinder is chosen as ${\rho _d} = 1920\ \textrm {kg}\ \textrm {m}^{-3}$ such that the capillary force and the weight of the displaced liquid (of density ${\rho _l} = 1000\ \textrm {kg}\ \textrm {m}^{-3}$) are not sufficient for the cylinder to remain afloat. Therefore, the cylinder sinks rapidly to become completely immersed in the bulk fluid. Length $h$ is the distance between the cylinder's centre and the undeformed free surface; $\beta$ is the angle between the position of the contact line and the vertical direction. In the simulation, the cylinder is initially half-immersed at the liquid–gas interface. The computational domain is chosen as $12D \times 4.8D$, which is sufficiently wide so that the capillary waves reflected by the sidewalls will not affect the sinking dynamics, and the cylinder diameter is resolved by $100\Delta x$. The four boundaries are set as solid walls with a contact angle of ${90^ \circ }$. The initial liquid–gas interface is beneath the top wall, with a distance $D$. The cylinder contact angle is set as the averaged value of the contact angles in the experiments (Vella et al. Reference Vella, Lee and Kim2006), i.e. $\theta = {111^ \circ }$ ($\phi = 1.0$ and $\delta \rho = 0.21$). To match the experimental conditions, we choose the Weber number as ${We} = {\rho _l}g{D^2}/\gamma = 3.7$ and the Ohnesorge number as ${Oh} = {\rho _l}{\nu _l}/\sqrt {{\rho _l}D\gamma } = 0.017$, with ${\rho _l} = 7.452$, gravity $g = 4.9 \times {10^{ - 6}}$, surface tension $\gamma = 0.098$ and ${\nu _l} = 0.02$ in lattice units.

Figure 10. Geometry set-up of the sinking of a cylinder from the liquid–gas interface.

Figure 11 shows snapshots of the sinking process, where the simulation results (figure 11b) are compared with the experimental results (figure 11a). A good agreement is observed at different instants. To be more quantitative, we also plot the evolutions of the cylinder's centre and the angle of the contact line relative to the centreline in figure 12. According to Vella et al. (Reference Vella, Lee and Kim2006) and He et al. (Reference He, Li, Huang, Hu and Wang2019), $h$ and $t$ have been normalized by the characteristic length $\sqrt {\gamma /{\rho _l}g}$ and the characteristic time ${(\gamma /{\rho _l}{g^3})^{1/4}}$, respectively. It is seen that the cylinder's centre experiences acceleration (due to gravity) in the beginning stage, then goes to a steady stage showing a linear drop of $h^*$ with $t^*$ and finally is subject to a slight deceleration (maybe due to the snap-off). Our present simulations agree very well with the experiments, as shown in figure 11(a). The angle $\beta$ decreases slightly in the beginning, which is as expected due to the hydrophobic condition, although the experimental data are missing at ${t^*} < 0.7$. Afterward, $\beta$ increases approximately linearly with time in both experiments and our simulations. The small differences between simulations and experiments may be attributed to the fact that the dynamic contact angle ${\theta _d} \in [{97^ \circ },{125^ \circ }]$ in experiments (Vella et al. Reference Vella, Lee and Kim2006) has been replaced by its mean value $\theta = {111^ \circ }$ in simulations. Nevertheless, the accuracy of the present simulation is within an acceptable range.

Figure 11. Snapshots of the sinking process at different time frames. (a) Experimental results from Vella et al. (Reference Vella, Lee and Kim2006) and (b) the corresponding simulation results of the present model.

Figure 12. Quantitative comparison between experiments (Vella et al. Reference Vella, Lee and Kim2006) and the present LBM–DEM simulations. (a) Time evolution of the cylinder's centre and (b) time evolution of the angle between the contact line and the vertical direction.

3.4. Self-assembly of three particles on a liquid–gas interface

In this subsection, we test the ability of our model for multi-particle liquid–gas–solid interactions by simulating the self-assembly of three hydrophilic particles on a liquid–gas interface. Following the previous simulations by the diffuse-interface immerse-boundary method (Liu, Gao & Ding Reference Liu, Gao and Ding2017) and LBM (Zhang et al. Reference Zhang, Liu and Zhang2020), three equal-size hydrophilic particles, marked as A, B and C, are initially half immersed at the liquid–gas interface with different distances between each other, as sketched in figure 13. The particle diameter is ${D}=5$ mm, resolved by 100 $\Delta x$, the contact angle is $\theta = {56.2^ \circ }$ and the system size is $30\ \textrm {mm} \times 20\ \textrm {mm}$. Non-slip velocity and neutral wetting conditions are set on all the wall boundaries. The liquid density, particle density and liquid viscosity are chosen as ${\rho _l} = 1000\ \textrm {kg}\ \textrm {m}^{-3}$, ${\rho _p} = 400\ \textrm {kg} \textrm {m}^{-3}$ and ${\nu _l} = 5 \times {10^{ - 6}}\ \textrm {m}^{2}\ \textrm {s}^{-1}$, respectively. Gravity is included with $g = - 9.8\ \textrm {m}\ \textrm {s}^{-2}$. In the LBM, the viscosity is set as ${\nu _l} = 0.04$. Based on the unit conversion, the time step in the LBM is $\Delta {t_{LBM}} = 5 \times {10^{ - 6}}\ \textrm {s}$, and the time step in the DEM is $\Delta {t_{DEM}} = 2.5 \times {10^{ - 6}}\ \textrm {s}$ by setting ${N_{sub}} = 2$. The surface tension based on the Laplace test in lattice units is $\gamma = 0.098$, equivalent to $0.075\ \textrm {N}\ \textrm {m}^{-1}$. The corresponding Weber and Ohnesorge numbers are ${We} = {\rho _l}g{D^2}/\gamma = 3.3$ and ${Oh} = {\rho _l}{\nu _l}/\sqrt {{\rho _l}D\gamma } = 0.033$, respectively.

Figure 13. Geometry set-up of the self-assembly of three hydrophilic particles on a liquid–gas interface.

Figure 14 shows the self-assembly process of the three particles at different times. From the figure, it is clear that there are three stages in the process: at the first stage, particle A moves towards B which remains almost stationary, while particle C slightly moves close to the right wall; when particle A meets with B, they move together towards the right until they collide with particle C; finally, the three particles assemble and adhere to the right wall. The physical origin of such a self-assembly behaviour can be explained with a force analysis, as done by Vella & Mahadevan (Reference Vella and Mahadevan2005). For a quantitative description, the time evolutions of the particle horizontal centres and vertical centres are plotted in figure 15. Figure 15(a) displays the above-mentioned three stages (marked as I, II and III), and the stage transitions (I/II and II/III) that happen at $t \approx 0.5\ \textrm {s}$ and $t \approx 0.85\ \textrm {s}$, respectively. During the assembly process, the particles float up and down, as evidenced by the oscillation behaviours in figure 15(b). In the end, particle A lies lower than B and C, which is due to the fact that the liquid–gas interface is lower to its left than its right and particle A must adaptively locate lower to obtain sufficient buoyancy. The present simulation results are in good agreement with the results reported in the literature (Liu et al. Reference Liu, Gao and Ding2017; Zhang et al. Reference Zhang, Liu and Zhang2020), which confirms the ability of our model for simulating complex liquid–gas–solid problems involving capillary force, gravity, moving boundary and multi-particle interactions.

Figure 14. Snapshots of the three-particle self-assembly process on a liquid–gas interface.

Figure 15. Evolution of normalized three-particle centres in the horizontal direction (a) and vertical direction (b) during their self-assembly process on a liquid–gas interface.