A Multiresolution Mesoscale Approach for Microscale Hydrodynamics
Andrea Montessori, Adriano Tiribocchi, Marco Lauricella, Fabio Bonaccorso, Sauro Succi
AA multi-resolution mesoscale approach for microscale hydrodynamics
Andrea Montessori ∗1 , Adriano Tiribocchi , Marco Lauricella , Fabio Bonaccorso , and SauroSucci Istituto per le Applicazioni del Calcolo CNR, via dei Taurini 19, Rome, Italy Center for Life Nano Science@La Sapienza, Istituto Italiano di Tecnologia, 00161 Roma, Italy Institute for Applied Computational Science, Harvard John A. Paulson School ofEngineering And Applied Sciences, Cambridge, MA 02138, United StatesJune 29, 2020
Abstract
We present a new class of multiscale schemes for micro-hydrodynamic problems based on a dual repre-sentation of the fluid observables.The hybrid model is first tested against the classical flow between two parallel plates and then appliedto a plug flow within a micron-sized striction and a shear flow within a micro-cavity.Both cases demonstrate the capability of the multiscale approach to reproduce the correct macroscopichydrodynamics also in the presence of refined grids (one and two levels), while retaining the correct thermalfluctuations, embedded in the multiparticle collision method.This work provides the first step towards a novel class of fully mesoscale hybrid approaches able to capturethe physics of fluids at the micro and nano scales whenever a continuum representation of the fluid fallsshort of providing the complete physical information, due to a lack of resolution and thermal fluctuations.
Complex flow simulations across scales represent a grand-challenge in computational physics and set a ceaselessquest for novel and efficient computational methods.The main computational challenge is represented by the spatial and time scale separations involved in aplethora of complex systems, from biology to soft matter, which can be addressed by developing efficient hybridapproaches able to bridge these spatiotemporal gaps.These methods, generally referred to as multiscale approaches, are based on composite computationalschemes relying on different levels of description of the system at hand.Many attempts have been deployed so far in order to tackle this problem, and a number of concurrent hybridapproaches have been proposed to investigate a vast spectrum of physics problems, including soft matter andmolecular fluids, fluctuating hydrodynamics, as well as both dilute and dense hydrodynamics [7, 8, 33, 28, 29, 10].In concurrent approaches, most of the physical domain is solved by employing a macroscopic (continuum)description, while only small selected portions are investigated through particle-based methods.The coupling and communication across the different regions take place within a handshaking region wherethe hydrodynamic information is exchanged between the two representations.In this work, we propose a hybrid numerical scheme based on the lattice Boltzmann method (LB) [22, 31,20, 23] and the multiparticle collision dynamics (MPCD) [12, 16], capable of predicting the correct macroscopichydrodynamics also in the presence of multilevel grids (one and two levels), while retaining the correct thermalfluctuations, embedded by default in the multiparticle collision method.This work is intended as a very first step towards the development of a novel class of hybrid and fullymesoscale approaches for complex flows across scales, designed in such a way as to capture the physics at thesmallest scales, whenever a continuum description alone falls short of providing the correct physical informationdue to a lack of resolution and thermal fluctuations. ∗ Electronic address: [email protected] ; Corresponding author a r X i v : . [ phy s i c s . c o m p - ph ] J un Methods
The lattice Boltzmann method [3] (LB) is based on a minimal version of the Bathnagar-Gross-Krook equation,in which a discrete set of probability distribution functions, flies freely along the links of a cartesian lattice andcollide on each node, relaxing towards a Maxwell-Boltzmann equilibrium. The LB equation reads as follows: f i ( ~x + ~c i ∆ t, t + ∆ t ) = f i ( ~x, t ) + ∆ tτ [ f eqi ( ρ, ~u ) − f i ( ~x, t )] (1)where f i ( ~x, t ) is the discrete distribution function representing the probability of finding a particle at position ~x and time t , streaming along the i − th lattice direction with a (discrete) velocity ~c i , i running over the latticedirections (i=0,...,b) [31, 24]. The lattice time step ∆ t is usually set to unity as well as the lattice spacing∆ x = 1.The left hand side of the eq.1 codes for the streaming process, namely the free-flight of particles along thelattice directions hopping from a site to the neighboring one.The right hand side implements the relaxation of the set of discrete distributions towards a local thermody-namic equilibrium f eqi , namely a truncated, low Mach number expansion of the Maxwell-Boltzmann distribution[31, 19].In this work,a second-order isotropic nine speed two dimensional lattice ( b = 8, D Q f eqi = w i ρ (cid:20) ~c i · ~uc s + ( ~c i · ~u ) c s − ~u · ~u c s (cid:21) (2)In the above equation, w i are the weights of the discrete equilibrium distribution functions, c s = P i w i ~c ix = 1 / ~u the fluid velocity. The relaxation parameter τ in equation (1) is linkedto the fluid kinematic viscosity via the linear relation ν = c s ( τ − ∆ t/
2) [31].The local (time-space) hydrodynamic moments, up to the order supported by the lattice in use are thedensity ( ρ ), linear momentum ( ρ~u ) and momentum flux tensor (Π) [6] and are evaluated by computing thestatistical moments of the discrete set of distribution as ρ = P i f i , ρ~u = P i f i ~c i and Π = P i f i ( ~c i ~c i − c s I),being I the identity matrix.
In the multiparticle collision dynamics, the fluid is modeled via a large number of pointlike particles of mass m , typically of the order of 10 − in two dimensions [12, 14]. During each time step, the system evolves viathe subsequent application of a streaming and a collision step. In the streaming step, the particles positions areupdated via a forward Euler step: ~r k ( t + δt ) = ~r k ( t ) + ~v k ( t ) δt (3)being ~r k the vector position and ~v k the velocity of the k − th particle and δt the value of the MPCD discretetime step.As per the collision, several approaches are available [26, 15, 2].A possible choice is to perform stochastic rotations of the particle velocities relative to the center-of-massmomentum. In this case, the multiparticle collision model is usually referred to as Stochastic Rotation Dynamics(SRD) [21].In the SRD, the domain is divided into cells of side a [12]. Multi-particle collisions are then performedwithin each cell, by rotating the velocity, ~υ k , of each k -th particle with respect to the velocity of the cell centerof mass, ~υ cm , of all particles in the cell: ~v k ( t + δt ) = ~v cm ( t ) + R [ ~v k ( t ) − v cm ( t )] (4)In the above equation, R is the rotation matrix, which rotates the particle velocities of a given cell by anangle ± α with uniform probability [ / , − / ].The kinematic viscosity of the SRD fluid [18] which, in two dimensions, reads as follows: ν = N k B T δta (cid:20) N ( N − e − N )(1 − cos (2 α )) (cid:21) + m (1 − cos ( α ))12 δt ( N − e − N ) (5)2eing N is the average number of particles in a collision cell.A stochastic rotation of the particle velocities is not the only possibility to perform multi-particle collisions.In particular, MPCD simulations can be performed directly in the canonical ensemble by employing an Andersenthermostat (MPCD-AT) [12, 26, 2].In the MPCD-AT, new relative velocities are generated during each computational step in each cell and thecollision step writes as [26]: ~v k ( t + δt ) = ~v cm ( t ) + δ~v rank = ~v cm ( t ) + ~v rank − N c X j ∈ cell ~v ranj (6)where ~v rank are random numbers sampled from a Gaussian distribution with variance p k B T /m and N c is theactual number of particles in the collision cell. In the above equation, the sum runs over all particles in a givencell.In the case of the MPCD with Andersen Thermostat, the kinematic viscosity is given by [12]: ν = k B T ∆ t (cid:20) NN − e − N − (cid:21) + a t (cid:20) N − e − N N (cid:21) (7)In this work we employed an MPCD with the Andersen thermostat since, although computationally moreexpensive than SRD, it features shorter relaxation times. This implies smaller number of time steps requiredfor transport coefficients to reach their asymptotic values [16]. As a consequence the restricition on maximumnumber of particles per cell (5 − N ) − [12].We wish to point out that the hydrodynamic limit of the multiparticle collision dynamics approach cor-responds to the fluctuating Navier Stokes equation for athermal and weakly-compressible fluids. Indeed, ithas been previously shown that the MPCD defines a discrete-time dynamics which has been demonstrated toreproduce the correct long-time hydrodynamics [12]. In addition to the conservation of mass and momentum,MPCD locally conserves energy as well, which enables simulations in the canonical ensemble, i.e. it fully incor-porates both thermal fluctuations and hydrodynamic interactions. The temperature of the fluid is then constant(apart from fluctuations) within the bulk and, consequently, the present MCPD model cannot describe thermaltransport phenomena. In this sense, the LB approach and MPCD simulation share some features in common.In particular, in both approaches, the system works at a well defined temperature ( c s in the lattice Boltzmannand k B T in the MPCD) but the MPCD allows the correct incorporation of thermal fluctuations while in theLB, being a pre-averaged (noise free) version of the lattice gas automata, thermal fluctuations are not included,unless one forces them into the scheme via a stochastic noise [1]. In this subsection, we detail the coupling strategy employed to exchange the hydrodynamic information betweenthe lattice Boltzmann and the multiparticle collision dynamics.First, we start from the simplest coupling case, namely ∆ x LB = ∆ x MP CD = a . In Figure 1, we report asketch of the hybrid simulation domain. The right region is the pure LB region (green area), while the left onerepresents the MPCD region (light blue area); the red area in the middle is the handshaking region, where theparticles velocities are generated from the LB values of the (local) linear momenta.It is worth noting that, in this work, we employ a one-way coupling, LB to MPCD, between the twoapproaches [13]. The LB runs across the entire domain, the information exchange between the two models isimplemented in the coupling region and the MPCD evolves only in selected regions of the domain.In other words, the information is transferred one-way only, from the LB to the MPCD domain and, conse-quently, the LB simulation is unaffected by the noise, which is built-in within the MPCD method.Moreover, in all the simulations performed, the coupling region covers up to two lattice nodes, which hasproved to be sufficient to obtain smooth transitions of the solutions between the two grids.The coupling algorithm proceeds as follows:1) LB step . The LB model runs across the whole domain over a single computational step ( full streamingand collision process ).2)
Information exchange step . In the coupling region the particles velocities are re-generated by drawingthem from a normal distribution, generated via the Box-Muller algorithm [5].3)
MPCD step . A full streaming and collision step of the MPCD-AT is performed within the MPCDregion.For the sake of clarity, a pseudo-code is reported in Algorithm 1.3
PCD-AT LBCoupling D x a= D x Figure 1: Representation of the hybrid domain. On the right region, the pure LB region (green area), andon the left, the high resolution (MPCD, light blue area) region where the streaming collision process of themultiparticle collision dynamics is performed. In the middle, the buffer region (red area) is drawn, where theLB → MPCD step is performed via the generation of particles velocity from the LB values of the (local) linearmomenta.
MPCD-AT LB Coupling a= D x/2 D x D x Figure 2: Representation of the hybrid multigrid (2 level) domain.
Algorithm 1
Coupling procedure : pseudo-code procedure LB-MPCD coupling LB step: call
LB collision call LB streaming
Information exchange step: for particles ∈ coupling region, draw ~v ( ~x, t ) from p π kBTm e − ( ~v − u )22 kBTm MPCD step call
MPCD streaming + boundary conditions call
MPCD collision 4e then extended the hybrid MPCD-LB approach to run on multigrid domains.Specifically, the LB code runs on the coarse grid, which covers the entire domain, while the MPCD runs ona finer grid with a = 1 / a = 1 / a = 1 / Algorithm 2
Coupling procedure : pseudo-code procedure LB-MPCD two level coupling LB step: call
LB collision call LB streaming
Bilinear interpolation step: for i, j ∈ coupling region do ~u coarse → ~u fine Information exchange step: for particles in the coupling region, draw ~v ( ~x, t ) from p π kBTm e − ( ~v − u )22 kBTm MPCD step call
MPCD streaming + boundary conditions call
MPCD collisionTo this purpose, we employed a simple bilinear interpolation.In the case of the refined simulation the integration time has been chosen so as to avoid sub-cycling inthe MPCD region. This is possible since the MPCD allows to set both the integration time step and the griddiscretization independently. Thus, in the case a = 1 / t = 2, insuch a way to guarantee synchronization between the LB and the MPCD simulation step.As pointed out above, the coupling procedure is performed via the generation of the particle velocities in theoverlapping region drawn from a Maxwell-Boltzmann distribution. This implies that, in the coupling region,the non-equilibrium part of the velocity distribution is set to zero. For the cases investigated in this work,this particular choice turned out to be effective, not compromising the accuracy of the coupling procedure, asdetailed in the following.It is worth noting that this coupling may become inaccurate in the presence of strong velocity gradients,where the non-equilibrium information of the velocity distribution cannot be ignored.In these cases, different approaches, for example a sampling from an Enskog distribution [11, 8], would bemore appropriate. It is of interest to compare the computational times associated with the LB and MPCD simulations.First, we note that a particle update (streaming and collision) takes about ∼ . µs/particle/step , comparableto the time needed to update a lattice unit in a single component LB code, 10 M LU P S (i.e. Mega LatticeUpdates per Second), which is a typical performance of a serial implementation of a single phase, LB code [19].In the case of the MPCD model, the running time scales roughly linearly with the total number of particlesand, in our simulations, it varies between 0 . s/step ÷ . s/step with the particle densities ranging between60 ÷
600 particles per bin.It is clear that an efficient parallelization able to exploit the computational capabilities of the latest super-computers is mandatory in order to make the coupling approach amenable to large scale simulations.The performance data reported above refer to a serial implementation of the code, run on a Intel XeonPlatinum 8176 CPU based on the Skylake microarchitecture (2 sockets each one with 28 cores with an installedDRAM of 720 GB ). The code was compiled by Intel Fortran Compiler version 18.0.2 with the recent AVX-512instruction set supported on Skylake architectures. As per the data structure employed, we opted for a Arrayof Structure for both the LB and MPCD impolementation (see [30]). A parallel version (openMP-MPI) of thehybrid approach is currently under development. 5 oupling Region LB MPCD-AT D x D x LB MPCD-AT D x D x LB MPCD-AT D x D x LB MPCD-AT D x D x u ( l u / s t e p ) w (lu) k B T=0.167 v x , p lbx v y , p lb y (a) (b)(c) (d) (e) Figure 3: (a-d) Hybrid approach, Poiseuille flow velocity field. The coupling region is denoted by dottedregion. The cell density of particles was varied between 60 ÷ N . Inset: time-space averaged velocitydistribution, comparison between LB(spider plot) and MPCD (2 d histogram). We start by testing the coupling procedure for the case of flow between two parallel plates, the main resultsbeing shown in figure 3.As reported above, the LB runs over a grid covering the entire fluid domain with the MPCD restricted tothe lower half of the fluid domain, the depth of the coupling zone extending over two grid nodes, as denoted bythe region included within the dashed line in figure 3.Four separate cases have been run for different values of the density of MPCD particles per cell, N = 60 ÷ ν = 0 . lu /step ( lu standing for lattice units).As evidenced in figure 3 (a-d), the LB and MPCD flow fields smoothly connect in the coupling region,even for the lowest value of particle density per cell ( N = 60). By increasing the cell density, the statisticalfluctuations associated to the flow field decrease accordingly, as shown in panel (e) of figure 3 which reports theLB and the MPCD velocity profiles.In the same plot,a comparison between the continuous (MPCD) and discrete (LB) velocity distributionfunctions at the steady state is reported (inset of figure 3(e)).The inset displays the time-space averaged velocity distribution. The MPCD velocity distribution (2 d colored field) presents its typical (shifted) Maxwell-Boltzmann shape while, the LB distributions (spider plot),is skewed, its skewness increasing for larger values of the mean channel velocity (i.e. as the average Machnumber increases).This departure from the Maxwell Boltzmann behaviour is due to the fact that the set of lattice distributionsis not allowed to shift, as instead occurs in the continuous case, but it occurs via a positive bias in the co-flowingdistributions. For this reason, the larger the non-equilibrium (represented herein by the large values of the meanchannel velocity) the more apparent is the skewness of the discrete set of distributions.As a further test, we have checked the velocity evolution in two sections of the channel (at L/ L/ L the height of the channel) during the transient comparing the results against the analytical solution[25] and we observed an excellent agreement between the numerical and analytical solutions.We further tested the capability of the hybrid approach to handle multi-resolution grids by performingsimulations of the Poiseuille test with near-wall coupling and grid refinement of the MPCD region .Specifically, we run two separate simulations with the following features:1) A reference case where the LB and the MPCD run on a grid with the same discretization (figure 4(a)).In this case the LB runs on a 40 ×
60 nodes grid and is coupled to the MPCD near the wall which runs on a40 ×
14 bins grid.2) In the multi-resolution implementation, the LB runs on a 20 ×
30 nodes grid ( halved with respect tothe previous case), and is coupled near the wall to the MPCD which runs on 39 ×
13 grid, with a twice finerdiscretization ( a = 1 /
2) with respect to the LB grid (figure 4(b)).6 B MPCD-AT
Coupling Region LB MPCD-AT D x D x D x D x/2 u(lu/step)width(lu) Figure 4: (left) Coupled LB-MPCD approach with the same grid discretization (fine grid simulation). (right)Hybrid approach with grid refinement at the wall. Inset : Comparison between the analytical solution of thePoiseuille flow near the wall and the MPCD solution.The first case is used as a a benchmark to test the ability of the coupling procedure to correctly reproducethe Poiseuille solution and to handle grid-refined domains.In this case, we set N = 1000 and, as before, the overlapping zone extends to two lattice units.The main results are reported in figure 4.Even in this case, the velocity profiles smoothly reconnect in the overlapping region, which is also quantita-tively confirmed by comparing the analytical solution near the wall against the MPCD averaged velocity profile,as shown in the inset of fig. 4, thus proving the effectiveness of the hybrid approach to handle multi-resolutiongrids. The multi-resolution algorithm was then employed to simulate the plug flow in a channel with a micrometricstriction.As an example, such geometries, embedded in more complex microfluidic chips, are widely used for biomi-crofluidics applications to perform deformability measurements on red blood cells [27]. The domain geometryis reported in figure 5(a). The fluid evolves under the effect of a constant body force ( ~g ) directed from left toright. In the LB simulation, the domain is periodic along the flow direction, while no-slip boundary conditionsare applied both at the upper and lower boundaries and on the striction walls. In the MPCD sub-domain,denoted by the white dotted rectangular area in figure 5, no-slip conditions are applied to the particles hittingwith the top and bottom walls while, on the right(outlet) and left (inlet) boundaries, particles which wouldescape from the MPCD domain are periodically re-injected at the opposite side, in order to exactly conservethe mass within the particle domain.As shown by the colorbar in figure 5, the maximum flow speed attained within the channel is ∼ · − (inlattice unit per time step). The channel is discretized with 10 lattice nodes and the fluid viscosity (in latticeunits) is ν = 0 . Re = Uhν ∼ µm height channel, with a maximum speed of ∼ cm/s , typicalvalues in microfluidic channels, the Reynolds is O (1), thus of the same order as the simulation at hand.The LB and MPCD domains are coupled at the inlet and outlet sections of the striction and the overlappingzone extends, as before, over two grid nodes. Within the coupling regions, the particle velocities are initializedby employing the procedure highlighted in Algorithm 1.Figures 5 (b) and (c) report a visual comparison between the flow fields within the striction obtained withthe full LB simulation and with the hybrid approach, while the plot below reports the velocity profiles taken atthree different sections of the microchannel. These results confirm that the hybrid approach is able to captureboth the dynamical features along the entrance and exit lengths of the striction inlet (see caption in figure 5)and the bulk dynamics within the channel.Moreover, the LB and MPCD velocity profiles show a fairly good agreement, thus proving the effectivenessof the coupling approach also in the presence of inlet and outlet velocity gradients.The same simulations were then performed by discretizing the micrometric channel in the MPCD simulationon finer grids. To this purpose, two sets of simulations have been performed by running the MPCD on a two-fold( a = 1 /
2) and four-fold ( a = 1 /
4) grid resolution. The main results are reported in figure 6, showing the flowfields ((a) and (b)) and the velocity profiles along the red-dashed vertical sections of the striction. In order tocompare the MPCD with the LB results, two full resolution LB simulations have been run as a reference.These results further demonstrate the ability of the multi-resolution approach to correctly reproduce thephysics of fluid within the micron-sized channel.Moreover, as stated in the introduction, the hybrid mesoscale approach not only allows to locally refine thegrid to obtain more accurate solutions of the hydrodynamic fields, but also naturally preserves the fluctuating7 u h h LB MPCD g wallwall Figure 5: (top panel) Sketch of the flow domain. (center panel) LB (left) and MPCD (right) solutions of theflow field within the microchannel. (bottom panel) comparison between LB and MPCD velocity profiles alongthe red-dashed vertical lines. Symbols stands for MPCD solutions. Open circles leftmost, crosses central andtriangles rightmost section. h u u h MPCD D x /2 MPCD D x /4 Figure 6: (left and right panels) MPCD velocity fields on grids with increased resolutions with the plot of thevelocity along vertical sections denoted by the dashed lines. Symbols denotes MPCD results.8 all u w Figure 7: Sketch of the channel with the micro-cavity.statistics, typical of confined fluids the micro- and nano-scales, which is absent in the non-fluctuating LB models[1, 9].
As a final case, we present a shear flow in a periodic channel with a microcavity placed at the bottom. Thiskind of geometry is typical of micropatterned surfaces [32] and rough walls in microdevices.Similar systems are widely employed in biological and medical applications for selective cell sorting andentrapment, enabling for tunable size-based cell capture in microvortices [17].A sketch of the domain is reported in figure 7, in which a shear is applied to the fluid by imposing amomemtum exchange boundary condition to the top boundary [4]. Periodic boundary conditions are appliedalong the flow directions, while no-slip condition is implemented on the bottom boundary and on the cavitywalls.The main simulation and geometrical parameters are the following:i) the LB domain is discretized with 40 ×
30 grid nodes and the square cavity has a side L = 20. Thekinematic viscosity is set to ν = 0 .
167 and the shear velocity imposed at the uppermost boundary is U = 0 . N = 1000.iii) As per the boundary conditions employed in the MPCD, the bottom and side walls are no slip, whilethe uppermost boundary is treated as a free slip wall as the velocity is fixed during the information exchangebetween the LB and MPCD grids.It is again of interest to compute the cavity Reynolds number Re = ULν ∼
6, typical of micro-cavities undershear [32, 17]. The MPCD run in the region inside the dashed rectangular perimeter (see figure 7)and thecoupling region coincides with the separation line between the microcavity and the main channel. As a firstcase, we run the coupled approach by using the same grid discretization both in the LB and in the MPCDregions. The first set of results is shown in figure 8. As shown in panel (a) and (b), the main vortex inside themicrocavity is accurately captured both in the LB and in the MPCD simulation. Further, the vortex positionsand sizes in both the representaion also show a good agreement.We then plotted the horizontal component and the magnitude of the velocity ( √ u + v ) along the height ( h )of the cavity (panels (c) and (d)). The agreement is evident and confirmes the smooth reconnection between theMPCD and LB solutions within the coupling region. Furthermore,the resolution in the MPCD region has beenincreased by a factor two and four( a = 1 / a = 1 / In this work a concurrent lattice Boltzmann-multiparticle collision dynamics multiscale/multigrid approach hasbeen developed and validated against the classical laminar flow through parallel plates and employed to simulatefluid flows in confined geometries at the micro-scale.In particular, for the Poiseuille flow across two parallel plates, several cases have been investigated, namely9
10 20 30-0.100.10.20.30.4 uu h h h LB MPCD
Figure 8: Comparison between the flow fields inside the micro-cavity computed with the LB(left) and theMPCD(right) solver. The plot reports a comparison of the horizontal component of the velocity (open circlesand dashed line) and of its magnitude (crosses and solid line),along the vertical axes of the cavity, between theLB and MPCD solutions. Symbols stand for the MPCD solutions. D x D x/2 LB MPCD
20 30 -0.1 -0.05 u h LB MPCD
Figure 9: Comparison between the flow fields inside the micro-cavity computed with the LB(left) and the refinedMPCD(right). The plot reports a comparison of the horizontal component of the velocity (open circles anddashed line) and its magnitude (crosses and solid line),along the vertical axes of the cavity, between the LB andMPCD solutions. Symbols stand for the MPCD solutions.10 B, D x MPCD, D x MPCD, D x/2 MPCD, D x/4 Figure 10: Flow fields inside the micro-cavity. Comparison between LB and MPCD with grid refinement.1) hybrid LB-MPCD with the coupling region positioned at the center of the channel2) Coupled LB-MPCD with one level near-wall grid refinementFurthermore, the multi-resolution approach has been used to simulate a plug flow in a microchannel with amicrometric striction and a shear flow in a channel with a micro-cavity placed on the bottom.The results presented are aimed at demonstrating the effectiveness of the hybrid approach to simulatefluctuating Navier-Stokes fluid dynamics even in the presence of complex geometries and velocity gradients nearthe coupling region.As stated in the introduction this work represents a first step towards a fully coupled (MPCD-LB) algorithm,which is currently under development. Nevertheless, a one-way coupling, like the one proposed in this work,can be employed in a number of applications where the effect of thermal fluctuations needs to be taken intoaccount only in small portions of the fluid domain. In this respect, an example is represented by the flow inthe microcavity. Indeed, within the microcavity, both hydrodynamic interactions and thermal fluctuations playan important role while outside, in the sheared channel, fluctuations are not essential and a lattice Boltzmannapproach can safely be applied.Further developments of the hybrid approach, currently ongoing, will aim at the simulation of complexfluid phenomena in confined systems, also in the presence of hydrodynamic interactions with nanoparticles, andwithin thin films at the fluid-fluid interface, wherein increased resolution and thermal fluctuations are bothneeded to capture the relevant microscale physics.
Data Availability
The code is available upon request.
Acknowledgements
A. M., M. L., A. T. and S. S. acknowledge funding from the European Research Council under the European196 Union’s Horizon 2020 Framework Programme (No. FP/2014-2020) ERC Grant Agreement No.739964(COPMAT).
References [1] R Adhikari, K Stratford, ME Cates, and AJ Wagner. Fluctuating lattice boltzmann.
EPL (EurophysicsLetters) , 71(3):473, 2005.[2] E Allahyarov and G Gompper. Mesoscopic solvent simulations: Multiparticle-collision dynamics of three-dimensional flows.
Physical Review E , 66(3):036702, 2002.[3] R. Benzi, S. Succi, and M Vergassola. The lattice boltzmann equation: theory and applications.
Phys.Rep. , 222(3):145–197, 1992. 114] M’hamed Bouzidi, Mouaouia Firdaouss, and Pierre Lallemand. Momentum transfer of a boltzmann-latticefluid with boundaries.
Physics of fluids , 13(11):3452–3459, 2001.[5] George EP Box. A note on the generation of random normal deviates.
Ann. Math. Stat. , 29:610–611, 1958.[6] Hudong Chen, Isaac Goldhirsch, and Steven A Orszag. Discrete rotational symmetry, moment isotropy,and higher order lattice boltzmann models.
Journal of Scientific Computing , 34(1):87–112, 2008.[7] Rafael Delgado-Buscalioni and Peter V Coveney. Hybrid molecular–continuum fluid dynamics.
Philo-sophical Transactions of the Royal Society of London. Series A: Mathematical, Physical and EngineeringSciences , 362(1821):1639–1654, 2004.[8] G. Di Staso, H. J. H. Clercx, S. Succi, and F. Toschi. Dsmc-lbm mapping scheme for rarefied and non-rarefied gas flows.
J. Comp. Sci. , 17:357–369, 2016.[9] Burkhard Dünweg, Ulf D Schiller, and Anthony JC Ladd. Statistical mechanics of the fluctuating latticeboltzmann equation.
Physical Review E , 76(3):036704, 2007.[10] Dmitry A Fedosov, Huan Lei, Bruce Caswell, Subra Suresh, and George E Karniadakis. Multiscale modelingof red blood cell mechanics and blood flow in malaria.
PLoS computational biology , 7(12):e1002270, 2011.[11] Alejandro L Garcia and Berni J Alder. Generation of the chapman–enskog distribution.
Journal of com-putational physics , 140(1):66–70, 1998.[12] G Gompper, T Ihle, DM Kroll, and RG Winkler. Multi-particle collision dynamics: A particle-basedmesoscale simulation approach to the hydrodynamics of complex fluids.
Advances in polymer science ,221(PreJuSER-3284):1–87, 2009.[13] D Hash and H Hassan. A decoupled dsmc/navier-stokes analysis of a transitional flow experiment. In , page 353, 1996.[14] T Ihle and DM Kroll. Stochastic rotation dynamics: A galilean-invariant mesoscopic model for fluid flow.
Physical Review E , 63(2):020201, 2001.[15] Thomas Ihle, E Tüzel, and Daniel M Kroll. Consistent particle-based algorithm with a non-ideal equationof state.
EPL (Europhysics Letters) , 73(5):664, 2006.[16] Raymond Kapral. Multiparticle collision dynamics: Simulation of complex systems on mesoscales.
Advancesin Chemical Physics , 140:89, 2008.[17] Reem Khojah, Ryan Stoutamore, and Dino Di Carlo. Size-tunable microvortex capture of rare cells.
Labon a Chip , 17(15):2542–2549, 2017.[18] N Kikuchi, CM Pooley, JF Ryder, and JM Yeomans. Transport coefficients of a mesoscopic fluid dynamicsmodel.
The Journal of chemical physics , 119(12):6388–6395, 2003.[19] Timm Krüger, Halim Kusumaatmaja, Alexandr Kuzmin, Orest Shardt, Goncalo Silva, and Erlend MagnusViggen. The lattice boltzmann method.
Springer International Publishing , 10:978–3, 2017.[20] Marco Lauricella, Simone Melchionna, Andrea Montessori, Dario Pisignano, Giuseppe Pontrelli, and SauroSucci. Entropic lattice boltzmann model for charged leaky dielectric multiphase fluids in electrified jets.
Physical Review E , 97(3):033308, 2018.[21] Anatoly Malevanets and Raymond Kapral. Mesoscopic model for solvent dynamics.
The Journal of chemicalphysics , 110(17):8605–8613, 1999.[22] A Montessori, M Lauricella, N Tirelli, and S Succi. Mesoscale modelling of near-contact interactions forcomplex flowing interfaces.
Journal of Fluid Mechanics , 872:327–347, 2019.[23] A Montessori, P Prestininzi, M La Rocca, and S Succi. Lattice boltzmann approach for complex nonequi-librium flows.
Physical Review E , 92(4):043308, 2015.[24] Andrea Montessori and Giacomo Falcucci.
Lattice Boltzmann modeling of complex flows for engineeringapplications . Morgan & Claypool Publishers, 2018.[25] Joseph P. Morris, Patrick J. Fox, and Yi Zhu. Modeling low reynolds number incompressible flows usingsph.
Journal of Computational Physics , 136(1):214 – 226, 1997.1226] Hiroshi Noguchi, Norio Kikuchi, and Gerhard Gompper. Particle-based mesoscale hydrodynamic tech-niques.
EPL (Europhysics Letters) , 78(1):10005, 2007.[27] Diana Pinho, Tomoko Yaginuma, and Rui Lima. A microfluidic device for partial cell separation anddeformability assessment.
BioChip Journal , 7(4):367–374, 2013.[28] Matej Praprotnik, Luigi Delle Site, and Kurt Kremer. Adaptive resolution scheme for efficient hybridatomistic-mesoscale molecular dynamics simulations of dense liquids.
Physical Review E , 73(6):066701,2006.[29] Matej Praprotnik, Luigi Delle Site, and Kurt Kremer. Multiscale simulation of soft matter: From scalebridging to adaptive resolution.
Annu. Rev. Phys. Chem. , 59:545–571, 2008.[30] Aniruddha G Shet, Shahajhan H Sorathiya, Siddharth Krithivasan, Anand M Deshpande, Bharat Kaul,Sunil D Sherlekar, and Santosh Ansumali. Data structure and movement for lattice-based simulations.
Physical Review E , 88(1):013314, 2013.[31] S. Succi.
The Lattice Boltzmann Equation: For Complex States of Flowing Matter . Oxford UniversityPress, 2018.[32] Hiroshi Suzuki, Ruri Hidema, and Yoshiyuki Komoda. Flow characteristics in a micro-cavity swept by avisco-elastic fluid.
Experimental Thermal and Fluid Science , 67:96–101, 2015.[33] HS Wijesinghe, RD Hornung, AL Garcia, and NG Hadjiconstantinou. Three-dimensional hybrid continuum-atomistic simulations for multiscale hydrodynamics.