A lattice Boltzmann study of particle settling in a fluctuating multicomponent fluid under confinement

We present mesoscale numerical simulations based on the coupling of the fluctuating lattice Boltzmann method for multicomponent systems with a wetted finite-size particle model. This newly coupled methodologies are used to study the motion of a spherical particle driven by a constant body force in a confined channel with a fixed square cross section. The channel is filled with a mixture of two liquids under the effect of thermal fluctuations. After some validations steps in the absence of fluctuations, we study the fluctuations in the particle’s velocity at changing thermal energy, applied force, particle size, and particle wettability. The importance of fluctuations with respect to the mean settling velocity is quantitatively assessed, especially in comparison with unconfined situations. Results show that the expected effects of confinement are very well captured by the numerical simulations, wherein the confinement strongly enhances the importance of velocity fluctuations, which can be one order of magnitude larger than what expected in unconfined domains. The observed findings underscore the versatility of the proposed methodology in highlighting the effects of confinement on the motion of particles in the presence of thermal fluctuations.


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 nonstraightforward 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.

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 li (x, t), representing the probability density to find a particle of fluid component l = A, B with kinetic velocity c i in the space-time location (x, t). Lattice velocities are discretized (i = 0, 1, . . . , 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 being ρ tot (x, t) = l=A,B ρ l (x, t) the total density. The evolution equation for the distribution functions over a unitary time step is given by The collision operator L is designed in such a way that it expresses the relaxation of the whole system toward a local Maxwellian distribution function f li (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 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 ) li (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 where G is a strength coefficient and ω 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 G = 1.5 (lattice Boltzmann units, lbu), and we simulate a bulk fluid with a majority of component A and fluid densities ρ A = 2.21 lbu (majority component) and ρ B = 0.09 lbu (minority component). The term ξ li (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 fluctuationdissipation 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 hydrodynam- where the bulk pressure P b and the chemical potential μ assume the form , g is external body force density acting on the fluid. The transport coefficients D and η are related to the relaxation times of the fluid. These will be fixed to D = 1/6 lbu and η = 0.383 lbu in all the simulations performed. The capital Greek symbols identify the stochastic stress (Σ tens ) and the stochastic diffusion (Ψ vec ) contributions to the equations of hydrodynamics where k B T is the thermal energy, while W tens and W 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 nonequilibrium fluctuations. We remark that the hydrodynamical equations reported above are obtained via a 1 When applying the Chapman-Enskog procedure to obtain the momentum equation, we can first consider the sum of all populations l f li (x) as an "effective" 1 component population; then, we can rely on the known results for onecomponent systems [28,40].
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 Δρ 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.

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 × 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, g, acting in the z direction. The resulting force on the particle is where ρ p is the particle density which is set to ρ p = 2ρ tot . The channel is resolved with L × L × L z = 60 × 60 × 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 θ = 120. 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 γ conf In unconfined (unconf) domains, one would expect the Stokes law for the friction γ unconf = 3πη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 veloc-ity under confinement and the Stokes' prediction for an unconfined particle driven by the same body force Notice that c m is a function of the aspect ratio d/L, with the property that lim d/L→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.

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 (z) 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 Fig. 2 We report the ratio, cm, 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 cm 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 with the fluctuation-dissipation balance [32,33]; thusas 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] m p dU (z) where m p = πd 3 ρ p /6 represents the particle's mass. The scalar term ζ(t) stands for the stochastic noise in compliance with the fluctuation dissipation theorem [54], i.e., 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 coarsegraining 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. The steady state predictions from Eq. 11 imply a Gaussian distribution for the velocity fluctuations First of all, we checked that U (z) conf = Fp γ 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 g = 5 · 10 −5 lbu and k B T = 0.45 · 10 −3 lbu. To check for the Gaussian shape, we report the PDF of the quantity x = (U The behavior of the velocity fluctuations at changing g, d/L and k B T, is analyzed in Figs. 4, 5, 6. 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 k B T (cfr. Fig. 4); we also observe the scaling ∼ (d/L) −3/2 for fixed g and k B T (cfr. Fig. 5) and the scaling ∼ (k B T) 1/2 for fixed d/L and g (cfr. Fig. 6). To be noticed that not only the scaling laws, but also the pre-factor 6/(πρ 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. The quantity cm 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 conf as a function of the body force density g, at changing d/L and thermal energy kBT. 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 g for Panel a and Panel b. Also, particle's velocity fluctuations show no dependency on the body force density 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  (Panel (d)). Theoretical prediction for the confined cases is given in Eq. 15. Theoretical prediction for unconfined cases is obtained using Eq. 15 with cm = 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 and U 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 ΔU unconf obtained by setting c m = 1 in Eq. 15.
In Fig. 7, we report ΔU 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 kBT mp (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.

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 × 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 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.