1 Introduction

Complex flow phenomena involving dispersions of particles moving in viscous fluids are of interest for their theoretical relevance in the framework of non-equilibrium statistical mechanics [1, 2]. Such phenomena are also relevant in a variety of applications, ranging from large [3] to small scales [4]. The corresponding theoretical description at the large scales hinges on the deterministic Navier–Stokes equations [5, 6], suitably coupled to the surface of the particles via hydrodynamic boundary conditions; these, in turn, account for the affinity of the particle toward the fluid and result in macroscopic properties, such as slip and wettability. The deterministic dynamics of the Navier–Stokes equations is however unsuitable for the description at smaller scales, where the assumptions of negligible fluctuations cease to be valid. Consistently, fluctuations need to be taken into account, see [7, 8] for some reference textbooks and [9,10,11] (plus references therein) for some recent works on the topic. In these conditions, mesoscale methods represent methods of choice [12]. By definition, mesoscale modeling is constructed at scales which are intermediate between the large scales and the small scales; hence, a suitable coarse-graining allows to recover the hydrodynamical description based on the Navier–Stokes equations. Additionally, one can enrich the modeling with nanoscale features like thermal fluctuations [7, 8]. Among all mesoscale methods, we are interested in the lattice Boltzmann models (LBM). Over the last decades, LBM have been successfully used to model complex hydrodynamic phenomena at large scales, such as particle suspensions [4, 13, 14], non-ideal fluids with phase transition and/or phase segregation [15,16,17,18,19,20], polymer flows [21,22,23], active matter [24] just to cite some prominent examples. Especially in the last decade, there has been a boost to push the applicability of LBM simulations toward nano-scales via the inclusion of thermal fluctuations [25,26,27,28,29, 31], designing the so-called fluctuating lattice Boltzmann methodology (FLBM). This methodology has been recently applied to the study of multicomponent fluids in the presence of thermal fluctuations [32, 33] and also to study the effects of thermally excited capillary waves on the break-up properties of a thin liquid ligament [34]. In this paper, we couple the FLBM with a wetted finite-size particle model [35]. This prompts the need of understanding how to choose the various tunable parameters in this complex system to obtain a stable numerical simulation; these parameters include the particle’s resolution, particle’s wettabilities, Shan-Chen forcing coupling coefficient, collision model relaxation time, fluids density ratio, gravity force, confinement ratio, and thermal energy. Some preliminary data on the particle diffusivity and mean square displacement without external driving forcing are presented in [30]. Therefore, here we focus on the quantitative assessment of the potentiality of such coupled methodology in modeling fluctuations of finite-size particles in the presence of confinement and external driving forces. Specifically, we quantitatively characterize the motion of a spherical colloidal particle driven by a constant body force in a confined nanofluidic channel. The channel is filled with a fluctuating multicomponent mixture of two fluids. The results of numerical simulations are compared with the expectations of a simplified hydrodynamical Langevin model for a finite-size particle comprising Gaussian noise and effective friction that accounts for the effects of confinement. We observe that numerical simulations capture very well the theoretical expectations at changing the various free parameters in the problem, i.e., the particle radius, the thermal energy, the driving force, and the particle wettability. In particular, the increasing importance of particle velocity fluctuations (with respect to its mean settling velocity) at increasing confinement is correctly modeled by the simulations. These numerical observations bear non-straightforward methodological importance, in view of the fact that simulations with FLBM cannot be granted a priori the “hydrodynamical limit” [32, 33], hence one has to verify a posteriori if the outcome of simulations is well-captured by hydrodynamical models.

The paper is organized as follows. In Sect. 2, we summarize the essential methodological aspects of the FLBM for multicomponent fluids and the coupling between particles and the multicomponent fluid. In Sect. 3, we will present the set-up for the numerical simulations, and we will present validation studies in the absence of thermal fluctuations, by comparing the settling velocity with previous experimental and numerical data. Results in the presence of thermal fluctuations are presented in Sect. 4. Conclusions will follow in Sect. 5.

2 Methodology

To model the bulk fluid, we consider LBM that allow for the simulations of multicomponent mixture of two components in the presence of thermal fluctuations [32]. We additionally introduce finite-size particles via a suitable coupling between the particle and the multicomponent fluid [35,36,37]. The essential technical details of the LBM used here are briefly summarized, the interested reader can refer to the reference works [32, 35,36,37] for more extensive technical coverage.

The multicomponent LBM considers the evolution equation of probability distribution functions, \(f_{l i}({\mathbf {x}},t)\), representing the probability density to find a particle of fluid component \(l=A,B\) with kinetic velocity \({{\varvec{c}}}_i\) in the space-time location \(({\mathbf {x}},t)\). Lattice velocities are discretized (\(i=0,1, \ldots , Q-1\)), and we employ the D3Q19 model, with \(Q=19\) velocity directions. The density of each component and the mixture velocity can be obtained via a proper coarse-graining in the kinetic velocity

$$\begin{aligned} \rho _{l}({\mathbf {x}}, t)= & {} \sum _{i=0}^{Q-1} f_{l i}({\mathbf {x}}, t),\nonumber \\ \rho _{\text {tot}}({\mathbf {x}}, t)\mathbf {v}({\mathbf {x}}, t)= & {} \sum _{i=0}^{Q-1}\sum _{l=A, B} f_{l i}({\mathbf {x}}, t){\mathbf {c}}_{i} \end{aligned}$$
(1)

being \(\rho _{\text {tot}}({\mathbf {x}},t)=\sum _{l=A, B} \rho _{l}({\mathbf {x}},t)\) the total density. The evolution equation for the distribution functions over a unitary time step is given by

$$\begin{aligned}&f_{l i}({\mathbf {x}}+{\mathbf {c}}_{i},t+1)-f_{l i}({\mathbf {x}}, t)={\mathfrak {L}}\left[ f_{l i}({\mathbf {x}},t )-f^{(eq)}_{l i}({\mathbf {x}},t )\right] \nonumber \\&\quad +S^{(F)}_{l i}({\mathbf {x}},t) + \xi _{l i}({\mathbf {x}},t) \quad l=A,B. \end{aligned}$$
(2)

The collision operator \({\mathfrak {L}}\) is designed in such a way that it expresses the relaxation of the whole system toward a local Maxwellian distribution function \(f^{(eq)}_{l i}({\mathbf {x}},t)\) [12, 38]. Technically, we make use of the MRT (multiple relaxation time) scheme [28, 39, 40]: the distribution functions are decomposed in modes (density, momentum, stress, etc.) and the action of \({\mathfrak {L}}\) consists in relaxing the different modes with different relaxation times [39]. The relaxation time of the momentum modes will determine the species diffusivity [28], whereas the relaxation time of the stress modes will determine the fluid viscosity [28, 39]. The term \(S^{(F)}_{l i}({\mathbf {x}},t)\) is a deterministic source term, accounting for the external body forces and the interactions between the two components. For the modeling of non-ideal interactions, we adopt the Shan-Chen formulation for multicomponent mixtures [41,42,43,44,45], where the force experienced by the fluid component l due to the surrounding fluid component \(l'\) can be written as

$$\begin{aligned} F_{l}({\mathbf {x}},t) = - {{{\mathcal {G}}}} \rho _{l}({\mathbf {x}},t) \sum _{l' \ne l}\sum _{i=0}^{Q-1} \omega _{i} \rho _{l'}({\mathbf {x}}+{\mathbf {c}}_{i},t) {\mathbf {c}}_{i} \quad \end{aligned}$$
(3)

where \({{{\mathcal {G}}}}\) is a strength coefficient and \(\omega _{i}\) a suitable weight needed to impose the isotropy in the interactions [41, 42, 44]. In all the simulations performed, we consider a non-ideal mixture with \(\mathcal{G}= 1.5\) (lattice Boltzmann units, lbu), and we simulate a bulk fluid with a majority of component A and fluid densities \(\rho _A=2.21\) lbu (majority component) and \(\rho _B=0.09\) lbu (minority component). The term \(\xi _{l i}({\mathbf {x}},t)\) is a stochastic force, which adds to the deterministic evolution a stochastic term. The stochastic terms are chosen in such a way that the conserved mass densities do not receive any stochastic force, while non-conserved modes receive a stochastic force in compliance with the fluctuation–dissipation relation [32]. The FLBM equations Eq. 2 imply evolution equations for macroscopic density and velocity. If we apply a Chapman–Enskog procedure [28, 40] by treating the stochastic terms as “generic” forcing terms, the macroscopic equations of a binary mixture in the presence of thermal fluctuations are recovered for the fluid densities and the hydrodynamical velocity \(\mathbf{v}^{(H)}=\mathbf{v}+(\mathbf{F}_{A}+\mathbf{F}_{B})/2 \rho _{\text {tot}}\) (superscript T means transposition) [8, 32]Footnote 1

$$\begin{aligned} \partial _t \rho _{\text {tot}} + {\varvec{\nabla }} \cdot (\rho _{\text {tot}} \mathbf{v}^{(H)})=0 \end{aligned}$$
(4)
$$\begin{aligned}&\partial _t (\rho _{\text {tot}} \mathbf{v}^{(H)})+{\varvec{\nabla }} (\rho _{\text {tot}} \mathbf{v}^{(H)} \mathbf{v}^{(H)}) \nonumber \\&\quad = -{\varvec{\nabla }} P_b + {\varvec{\nabla }} \cdot [\eta ({\varvec{\nabla }} \mathbf{v}^{(H)} + ({\varvec{\nabla }} \mathbf{v}^{(H)})^T )+{\varvec{\Sigma }}_{{\text {tens}}}] + \rho _{\text {tot}} \mathrm {g}\nonumber \\ \end{aligned}$$
(5)
$$\begin{aligned} \partial _t \rho _A+{\varvec{\nabla }} \cdot (\rho _A \mathbf{v}^{(H)})={\varvec{\nabla }} \cdot \left[ D {\varvec{\nabla }} \mu + {\varvec{\Psi }}_{{\text {vec}}} \right] \end{aligned}$$
(6)

where the bulk pressure \(P_b\) and the chemical potential \(\mu \) assume the form \(P_b=\frac{1}{3} \rho _{\text {tot}} +\frac{\mathcal{G}}{3} \rho _A \rho _B\) and \(\mu =\frac{1}{3} \log \rho _A - \frac{1}{3} \log \rho _B+\frac{1}{3} {{{\mathcal {G}}}} (\rho _A -\rho _B)\) [32], \(\mathrm {g}\) is external body force density acting on the fluid. The transport coefficients D and \(\eta \) are related to the relaxation times of the fluid. These will be fixed to \(D=1/6\) lbu and \(\eta =0.383\) lbu in all the simulations performed. The capital Greek symbols identify the stochastic stress (\({\varvec{\Sigma }}_{{\text {tens}}}\)) and the stochastic diffusion (\({\varvec{\Psi }}_{{\text {vec}}}\)) contributions to the equations of hydrodynamics

Fig. 1
figure 1

Sketch of the setup for the particle settling numerical simulations. The computational box is a rectangular parallelepiped of height \(L_z\) and square base \(L \times L\). The solvent fluid is a fluctuating mixture of two non-ideal components, A and B, with majority of the component A in the bulk phase. The whole system is under the effect of the body force density, \(\mathrm {g}\), acting in the z direction. The mixture and the particle are simulated using lattice Boltzmann models (LBM) on a regular three dimensional lattice (cfr. Sect. 2). The LBM implementation is further equipped with thermal fluctuations (fluctuating LBM, FLBM [32]) to mimic the effect of noise at the small scales. The particle’s velocity is tracked as a function of time and we quantify its statistical properties (average \(\langle U^{\mathrm {(z)}}_{\mathrm {conf}}\rangle \) and fluctuations \(\Delta U^{\mathrm {(z)}}_{\mathrm {conf}}\)) in the statistically steady state

Full size image
$$\begin{aligned} {\varvec{\Sigma }}_{{\text {tens}}}=\sqrt{\eta {\mathrm {k_B T}}}(\mathbf{W_{{\text {tens}}}}+\mathbf{W}_{{\text {tens}}}^T) \quad {\varvec{\Psi }}_{{\text {vec}}}=\sqrt{2 D {\mathrm {k_B T}}} \mathbf{W_{{\text {vec}}}}\nonumber \\ \end{aligned}$$
(7)

where \({\mathrm {k_B T}}\) is the thermal energy, while \(\mathbf{W}_{{\text {tens}}}\) and \(\mathbf{W}_{{\text {vec}}}\) are a Gaussian tensor and a Gaussian vector with independent and uncorrelated components and variance equal to unity. In both homogeneous and heterogeneous systems, validations for the stochastic term have been done in previous studies [32, 33]. Specifically, for homogeneous systems, it was demonstrated that the correlations of the hydrodynamical fields are exactly those that can be predicted from the fluctuating hydrodynamic equations that we write above. For heterogeneous systems, the model has been successfully benchmarked against capillary fluctuations at non-ideal interfaces. In a more recent study [33], the model has also been shown to reproduce quantitative details of non-equilibrium fluctuations. We remark that the hydrodynamical equations reported above are obtained via a Chapman–Enskog procedure. This requires fields that slowly vary in time and space, hence in the presence of fluctuations it may become questionable. We will discuss more in details these issues while presenting the results of numerical simulations.

For the LBM modeling of the particle, we follow References [35,36,37]. The particle is modeled on the lattice, by declaring the fluid nodes belonging to the particle (“particle nodes”), as sketched in Fig. 1. The motion of the particle is determined by Newton’s equation [37], and the evolution of the finite-size particle is solved with the leap-frog algorithm [46]. The integration of the leap-frog algorithm has been validated in previous studies for the finite-size particles in turbulent channel flows [47, 48]. The bounce back boundary condition is implemented at the interface between the particle and the fluid [36]. During the bounce back procedure, the particle exchanges the momentum with the surrounding fluid. Due to the particle movement in the fluid, there will be the creation of new particle nodes which originally were fluid nodes (cover-nodes behavior). Analogously, the movement of the particle can delete the particle nodes and create new fluid nodes (uncover-nodes behavior). In order to impose the total mass conservation, we implement the mass correction algorithm described in [35]. Also, we introduce a virtual fluid layer [35] at the interface between the particle and the fluid to be able to tune the particle’s wettability. In such a layer, the fluid densities are set equal to the average densities of the neighboring fluid nodes, plus a correction \(\Delta \rho \) that is instrumental to model the affinity of the particle toward the two components. The wettability properties described in the following (i.e., hydrophobic, neutral, hydrophilic) refer to the affinity of the particle toward the majority component in the bulk phase.

3 Numerical set-up and validation

The set-up for the numerical simulations is sketched in Fig. 1. A particle with diameter d is placed in a long channel with square cross section \(L \times L\). The particle is initially placed with its center of mass lying in the center of the square cross section and is driven by a body force density, \(\mathrm {g}\), acting in the z direction. The resulting force on the particle is

$$\begin{aligned} \mathrm {F_p}=\frac{4}{3}\pi \left( \frac{d}{2}\right) ^3(\rho _{\mathrm {p}}-\rho _{\text {tot}})\mathrm {g}\end{aligned}$$
(8)

where \(\rho _{\mathrm {p}}\) is the particle density which is set to \(\rho _{\mathrm {p}}=2\rho _{\text {tot}}\). The channel is resolved with \(L \times L \times L_z\) = \(60 \times 60 \times 900\) lbu. The channel is closed with walls in all directions; this choice is instrumental to fully appreciate the effects of confinement. A neutral wettability boundary condition is chosen for all the bounding walls, while three different wettabilities are considered at the interface between the particle and the fluid: these correspond to wetting angles \(\theta =120.5\), \(\theta =90^{\circ }\), \(\theta =55.0\) and will be denoted hereafter as “hydrophobic,” “neutral” and “hydrophilic.” Due to the wide range of parameters, it is highly challenging to fully validate the theoretical prediction with numerical simulations. In this work, we spent around 1.4 Million computing hours to accomplish the validation. The square cross section is kept fixed in all numerical simulations, while different particle’s diameters are considered. Different values of the thermal energy and driving force are also simulated. The simulations parameters are chosen in the following ranges: \(d/L \in [0.133:0.67]\), \({\mathrm {k_B T}}\in [1\cdot 10^{-5}:0.45 \cdot 10^{-3}]\) lbu, \(\mathrm {g}\in [5 \cdot 10^{-7}: 5 \cdot 10^{-5}]\) lbu.

We have first validated the numerical set-up without thermal fluctuations. To this aim, we measured the steady settling velocity of the particle at changing the particle diameter. The steady settling velocity in the confined (conf) channel will be proportional to the driving force and inversely proportional to the friction \(\gamma _{\mathrm {conf}}\)

$$\begin{aligned} U^{\mathrm {(z)}}_{\mathrm {conf}}=\frac{\mathrm {F_p}}{\gamma _{\mathrm {conf}}}. \end{aligned}$$
(9)

In unconfined (unconf) domains, one would expect the Stokes law for the friction \(\gamma _{\mathrm {unconf}}=3 \pi \eta d\). However, it is known from the literature that confinement enhances friction and reduces the settling velocity in comparison with the unconfined cases [49,50,51,52,53]. We therefore focused the attention on the ratio \(c_m\) as a function of d/L. \(c_m\) is defined as the ratio between the particle’s settling velocity under confinement and the Stokes’ prediction for an unconfined particle driven by the same body force

$$\begin{aligned} c_m=\frac{U^{\mathrm {(z)}}_{\mathrm {conf}}}{U^{\mathrm {(z)}}_{\mathrm {unconf}}}. \end{aligned}$$
(10)

Notice that \(c_m\) is a function of the aspect ratio d/L, with the property that \(\lim _{d/L\rightarrow 0} c_m=1\). Results are in good agreement with the experimental observations. Discrepancies which may be observed with coarser grids (results from [37]) become essentially negligible with our finer grids. We also observe that there is almost no dependency on the wettability condition.

Fig. 2
figure 2

We report the ratio, \(c_m\), between the particle’s settling velocity under confinement and the Stokes’ prediction for an unconfined particle driven by the same body force (cfr. Eq. 10). The ratio \(c_m\) is considered as a function of the degree of confinement d/L (cfr. Fig. 1) and for different wettability boundary conditions at the particle’s surface: hydrophilic (blue squares), neutral (green triangles), hydrophobic (red circles). Results are compared with the experimental results in [53] and with the numerical results in [37]. Our numerical investigation agrees well with the previous numerical and experimental results

Full size image

4 Results and discussions

After the validation of the results in the absence of thermal fluctuations, we switched on the thermal noise in the LBM simulations and studied the corresponding fluctuations in the particle’s velocity \(U^{\mathrm {(z)}}_{\mathrm {conf}}\) in the confined environment. As we have seen in the previous section, the friction acting on the particle is clearly affected by confinement, and it increases with respect to the unconfined case. This increase in friction is quantitatively well-reproduced by the simulations (cfr. Fig. 2). Fluctuations are added in the LBM in compliance with the fluctuation–dissipation balance [32, 33]; thus—as a first guess—one could invoke a simplified picture based on a Langevin equation for the particle’s velocity in the direction of the body force [54]

$$\begin{aligned} m_{\mathrm {p}}\frac{dU^{\mathrm {(z)}}_{\mathrm {conf}}(t)}{dt} + \gamma _{\mathrm {conf}}U^{\mathrm {(z)}}_{\mathrm {conf}}(t) = F_{\mathrm {p}} + \zeta (t), \end{aligned}$$
(11)

where \(m_{\mathrm {p}}=\pi d^3 \rho _{\mathrm {p}}/6\) represents the particle’s mass. The scalar term \(\zeta (t)\) stands for the stochastic noise in compliance with the fluctuation dissipation theorem [54], i.e.,

$$\begin{aligned} \langle \zeta (t) \zeta (t') \rangle = 2 \gamma _{\mathrm {conf}}{\mathrm {k_B T}}\delta (t-t'). \end{aligned}$$

Establishing the correspondence between the mesoscale FLBM dynamics (cfr. Sect. 2) and Eq. 11 is not simple from the methodological point-of-view. Indeed, interpretations based on hydrodynamical equations for lattice Boltzmann simulations rely on a coarse-graining view in the kinetic velocity space and invoke some multi-scale expansion technique (e.g., Chapman–Enskog [12, 32]) to find the corresponding hydrodynamical equations. By treating the stochastic source terms as “generic” one can surely carry out the detailed expansion calculations and derive fluctuating hydrodynamics equations (cfr. Sect. 2). It has to be noted, however, that such a procedure typically requires that fields under study slowly vary in time and space; thus, the “equivalence” between the FLBM simulations and the simple model Eq. 11 could well fail. One is therefore left with the need of assessing a posteriori the correctness of numerical simulations and if they match the predictions of Eq. 11 without fitting parameters. Based on this view, we started to analyze the statistical properties in the particle’s velocity.

Fig. 3
figure 3

Panel a: Standardized PDFs of particle’s settling velocity at fixed body force density \(\mathrm {g}=5\cdot 10^{-5}\) (lbu) and fixed thermal energy \({\mathrm {k_B T}}=0.45\cdot 10^{-3}\) (lbu). Three different wettabilities were chosen: hydrophilic, neutral and hydrophobic. To make the PDFs comparable with a standard Gaussian distribution, we have considered rescaled variables with zero mean and unitary variance. Data come from different values of the particle diameter \(d/L=0.13, 0.47, 0.67\). Results are matched well with standard Gaussian distribution. Panel b: we report \(c_m\) as function of d/L in three cases: hydrophilic, neutral and hydrophobic. The quantity \(c_m\) is computed as the ratio between the average particle’s settling velocity under confinement and the Stokes’ prediction for an unconfined particle driven by the same body force. Error bars are estimated from standard deviation of the particle’s settling velocity fluctuations

Full size image

The steady state predictions from Eq. 11 imply a Gaussian distribution for the velocity fluctuations

$$\begin{aligned} P(U^{\mathrm {(z)}}_{\mathrm {conf}})=\sqrt{\frac{1}{2\pi \sigma _{\mathrm {p}}^2}} e^{-\frac{(U^{\mathrm {(z)}}_{\mathrm {conf}}-\langle U^{\mathrm {(z)}}_{\mathrm {conf}}\rangle )^2}{2 \sigma _{\mathrm {p}}^2}} \quad \end{aligned}$$
(12)

where

$$\begin{aligned} \langle U^{\mathrm {(z)}}_{\mathrm {conf}}\rangle =\frac{F_{\mathrm {p}}}{\gamma _{\mathrm {conf}}} \quad \sigma _{\mathrm {p}}^2=\frac{{\mathrm {k_B T}}}{m_{\mathrm {p}}}. \end{aligned}$$
(13)

First of all, we checked that \(\langle U^{\mathrm {(z)}}_{\mathrm {conf}}\rangle =\frac{F_{\mathrm {p}}}{\gamma _{\mathrm {conf}}}\) holds and that the results are compatible with a Gaussian shape. We report in Fig. 3 some representative results for different d/L and different wettabilities, while keeping the body force density and the thermal energy fixed to \(\mathrm {g}=5\cdot 10^{-5}\) lbu and \({\mathrm {k_B T}}=0.45\cdot 10^{-3}\) lbu. To check for the Gaussian shape, we report the PDF of the quantity \(x=(U^{\mathrm {(z)}}_{\mathrm {conf}}-\langle U^{\mathrm {(z)}}_{\mathrm {conf}}\rangle )/\sigma _{\mathrm {p}}\). As can be seen, \(\langle U^{\mathrm {(z)}}_{\mathrm {conf}}\rangle =\frac{F_{\mathrm {p}}}{\gamma _{\mathrm {conf}}}\) holds and the numerical results collapse well on the Gaussian shape \(f(x)=e^{-x^2/2}/\sqrt{2 \pi }\). Then, we proceeded in characterizing the dependency of the particle’s velocity fluctuations \(\Delta U^{\mathrm {(z)}}_{\mathrm {conf}}=\sqrt{\sigma _{\mathrm {p}}^2}\) on the three parameters d/L, \({\mathrm {k_B T}}\) and \(\mathrm {F_p}\). From Eq. 13 and \(m_{\mathrm {p}}=\pi d^3 \rho _{\mathrm {p}}/6\), one gets

$$\begin{aligned} \Delta U^{\mathrm {(z)}}_{\mathrm {conf}}= \sqrt{\frac{6}{\pi \rho _{\mathrm {p}}L^3}}\,({\mathrm {k_B T}})^{1/2}\,(d/L)^{-3/2}. \end{aligned}$$
(14)

The behavior of the velocity fluctuations at changing \(\mathrm {g}\), d/L and \({\mathrm {k_B T}}\), is analyzed in Figs. 4, 5, 6.

Fig. 4
figure 4

Particle’s velocity fluctuations \(\Delta U^{\mathrm {(z)}}_{\mathrm {conf}}\) as a function of the body force density \(\mathrm {g}\), at changing d/L and thermal energy \({\mathrm {k_B T}}\). Hydrophilic (blue squares), neutral (green triangles) and hydrophobic (red circles) cases have been shown. Results are compared with the theoretical prediction given in Eq. 14. Panel a shows the results under the degree of the confinement \(d/L=0.13\), and Panel b presents the results at \(d/L=0.47\). Our simulation data fit well with the theory at all \(\mathrm {g}\) for Panel a and Panel b. Also, particle’s velocity fluctuations show no dependency on the body force density \(\mathrm {g}\). When the particle reaches the stationary state, we equally split the data set in five time intervals. Error bars are the standard deviations from different groups of the configurations

Full size image
Fig. 5
figure 5

Particle’s velocity fluctuations as a function of d/L at changing the body force density \(\mathrm {g}\) and the thermal energy \({\mathrm {k_B T}}\). Hydrophilic (blue squares), neutral (green triangles) and hydrophobic (red circles) cases have been shown. Results are compared with the theoretical prediction given in Eq. 14. Panel a shows results at \(\mathrm {g}=5\cdot 10^{-7}\mathrm {(lbu)}\), and Panel b shows results at the largest body force density \(\mathrm {g}=5\cdot 10^{-5}\mathrm {(lbu)}\). Three different lines are the theoretical predictions from Eq. 14 at \({\mathrm {k_B T}}=1\cdot 10^{-5}, 1\cdot 10^{-4}, 0.45\cdot 10^{-3}\mathrm {(lbu)}\). Our simulation data fit well with the theory at all d/L. When the particle reaches the stationary state, we equally split the data set in five time intervals. Error bars are the standard deviations from different groups of the configurations

Full size image
Fig. 6
figure 6

Particle’s velocity fluctuations as a function of the thermal energy \({\mathrm {k_B T}}\) at changing the body force density \(\mathrm {g}\) and d/L. Hydrophilic (blue squares), neutral (green triangles) and hydrophobic (red circles) cases have been shown. The back dotted line is the theoretical predictions given in Eq. 14. Panel a and b show that simulation data match well with the theory at all \({\mathrm {k_B T}}\) under the body force density \(\mathrm {g}=5\cdot 10^{-5}, 5\cdot 10^{-6}, 5\cdot 10^{-7}\mathrm {(lbu)}\). When the particle reaches the stationary state, we equally split the data set in five time intervals. Error bars are the standard deviations from different groups of the configurations. The subplots of Panel a and b checks the \((\Delta U^{\mathrm {(z)}}_{\mathrm {conf}})^2/{\mathrm {k_B T}}\), at changing d/L the value remain constant which is equal to \(1/m_{\mathrm {p}}\)

Full size image
Fig. 7
figure 7

Particle’s velocity fluctuations normalized to the average velocity in both confined (conf) and unconfined (unconf) environments, as a function of the normalized particle’s diameter d/L. Hydrophilic (blue squares), neutral (green triangles) and hydrophobic (red circles) cases have been presented. We change both the body force density \(\mathrm {g}\) and the thermal energy \({\mathrm {k_B T}}\): \(\mathrm {g}=5 \cdot 10^{-7}\) lbu, \({\mathrm {k_B T}}=1\cdot 10^{-5}\) lbu (Panel a), \(\mathrm {g}=5 \cdot 10^{-5}\) lbu, \({\mathrm {k_B T}}=1\cdot 10^{-5}\) lbu (Panel b), \(\mathrm {g}=5 \cdot 10^{-7}\) lbu, \({\mathrm {k_B T}}=0.45 \cdot 10^{-3}\) lbu (Panel c), \(\mathrm {g}=5 \cdot 10^{-5}\) lbu, \({\mathrm {k_B T}}=0.45 \cdot 10^{-3}\) lbu (Panel (d)). Theoretical prediction for the confined cases is given in Eq. 15. Theoretical prediction for unconfined cases is obtained using Eq. 15 with \(c_m=1\). When the particle reaches the stationary state, we equally split the data set in five time intervals. Error bars are the standard deviations from different groups of the configurations

Full size image

The results are also compared with the prediction of Eq. 14. As predicted by Eq. 14, the velocity fluctuations are independent of the body force for fixed d/L and \({\mathrm {k_B T}}\) (cfr. Fig. 4); we also observe the scaling \(\sim (d/L)^{-3/2}\) for fixed \(\mathrm {g}\) and \({\mathrm {k_B T}}\) (cfr. Fig. 5) and the scaling \(\sim ({\mathrm {k_B T}})^{1/2}\) for fixed d/L and \(\mathrm {g}\) (cfr. Fig. 6). To be noticed that not only the scaling laws, but also the pre-factor \(\sqrt{6/(\pi \rho _{\mathrm {p}}L^3)}\) in Eq. 14 matches very well with the numerical observations. Overall, hydrophilic (blue squares), neutral (green triangles), and hydrophobic (red circles) results are well overlapping. We observe little dependency on the particle’s wettability.

Finally, in Fig. 7, we consider velocity fluctuations normalized to the mean settling velocity in hydrophilic, neutral and hydrophobic cases, to highlight the impact of the fluctuations with respect to the characteristic order magnitude of the velocity. Based on Eqs. 14, 10 and \(U^{\mathrm {(z)}}_{\mathrm {unconf}}=\mathrm {F_p}/\gamma _{\mathrm {unconf}}\), \(\gamma _{\mathrm {unconf}}=3 \pi \eta d\), we obtain

$$\begin{aligned} \frac{\Delta U^{\mathrm {(z)}}_{\mathrm {conf}}}{U^{\mathrm {(z)}}_{\mathrm {conf}}}= & {} 18 \sqrt{\frac{6 \eta ^2}{(\rho _A+\rho _B)^2 \pi \rho _{\mathrm {p}}L^{7}}}\nonumber \\&\frac{ \,({\mathrm {k_B T}})^{1/2} \,(d/L)^{-7/2}\,\mathrm {g}^{-1}}{c_m} \end{aligned}$$
(15)

where we have related the particle’s velocity to the unconfined velocity via the ratio \(c_m\) (cfr. Eq. 10). To gain insight on the importance of confinement, we also compared the present results with the unconfined predictions \(\Delta U^{\mathrm {(z)}}_{\mathrm {unconf}}/U^{\mathrm {(z)}}_{\mathrm {unconf}}\) obtained by setting \(c_m=1\) in Eq. 15.

In Fig. 7, we report \(\Delta U^{\mathrm {(z)}}_{\mathrm {conf}}/U^{\mathrm {(z)}}_{\mathrm {conf}}\), \(\Delta U^{\mathrm {(z)}}_{\mathrm {unconf}}/U^{\mathrm {(z)}}_{\mathrm {unconf}}\) at changing d/L for selected values of \(\mathrm {g}\) and \({\mathrm {k_B T}}\) in hydrophilic, neutral and hydrophobic cases. We also compare with the theoretical predictions obtained from Eq. 15. The numerical data are well in agreement with the theory for all values of d/L. The unconfined theory is well-reproduced only at small d/L, as expected. Notice that for the largest d/L we observe a dramatic enhance of the importance of confinement, which is about one magnitude higher than the unconfined theory. This is expected based on the solution of Eq. 13.

Before closing this section, we stress once more that all the theoretical predictions that we have verified in the numerical simulations are the natural consequence of the simplified Langevin equation Eq. 11, where the noise follows the fluctuation dissipation theorem [54] and the friction accounts for confinement. The theoretical outcomes of such scenario predict that the mean settling velocity is reduced by confinement (cfr. Fig. 2) and the fluctuations around the mean settling velocity are unchanged and equal to \(\frac{{\mathrm {k_B T}}}{m_{\mathrm {p}}}\) (cfr. Eq. 13). What we target here is not the solution of such equation, but the verification that the results of the adopted methodology can be explained in the hydrodynamical limit with such an equation.

5 Conclusions

We studied the settling of a spherical particle with diameter d in a fluctuating multicomponent fluid. The system is driven by a constant body force in a confined channel with a square cross-sectional area \(L \times L\). Our simulations hinge on the fluctuating lattice Boltzmann methodology (FLBM) coupled with a finite-size particle model with tunable wettability [35]. This methodological coupling has never been tested in the literature: due to the fluctuating nature of the lattice Boltzmann populations, it requires careful numerical verification in the assessment of its hydrodynamical properties. We have first validated the numerical set-up in the absence of thermal fluctuations. In agreement with earlier numerical studies [37] on single component LBM, our numerical simulations with the multicomponent LBM well-reproduce the frictional properties of a confined particle [53]. We have then switched on thermal fluctuations, and we have systematically characterized the steady-state statistical properties (i.e., average and fluctuations) of the particle’s velocity at changing the thermal energy \({\mathrm {k_B T}}\), the degree of confinement d/L, the body force. The results of the numerical simulations show a neat matching with the predictions of a simplified Langevin-type scenario, accounting for the motion of a particle subject to the linear frictional law induced by confinement and in the presence of a stochastic force satisfying the fluctuation-dissipation theorem [55]. We think this is a non-straightforward result since the coarse-grained description of FLBM requires some “hydrodynamical assumption,” and this could be well-violated by the presence of mesoscale fields that do not vary smoothly in space and time. On a quantitative basis, results in the presence of confinement show that the numerical tool is quite versatile in handling quantitative changes in frictional properties across orders of magnitude. Correspondingly, the measured ratio between the velocity fluctuations and the mean velocity comes out to be dramatically increased in the presence of confinement. Taken all together, the “ensemble” of simulations here proposed underscore the robustness and versatility of the proposed methodology in concrete applications involving the motion of colloidal particles in the presence of confinement and multicomponent fluids.

In future work, it would be interesting to explore regimes where noise effects produce a Brownian time larger than the Stokes’ time. This would allow also to study lubrication effects coming from particle/wall interactions. Another follow-up could be represented by the numerical simulations of the particle motion settling at the interface separating two immiscible fluids [56], where a change of wettability is expected to lead to more sizeable effects than those observed in the present study. This makes the presented numerical results particularly relevant on the future perspective of achieving a further upgrade of the FLBM simulations as quantitative tools for the study of complex fluids with colloidal particles.