A Multiresolution Mesoscale Approach for Microscale Hydrodynamics

A new class of multiscale scheme is presented for micro-hydrodynamic problems based on a dual representation of the fluid observables. The hybrid model is first tested against the classical flow between two parallel plates and then applied to a plug flow within a micrometer-sized striction and a shear flow within a microcavity. Both cases demonstrate the capability of the multiscale approach to reproduce the correct macroscopic hydrodynamics also in the presence of refined grids (one and two levels), while retaining the correct thermal fluctuations, embedded in the multiparticle collision method. This provides the first step toward a novel class of fully mesoscale hybrid approaches able to capture the physics of fluids at the micro- and nanoscales whenever a continuum representation of the fluid falls short of providing the complete physical information, due to a lack of resolution and thermal fluctuations.


Introduction
Complex flow simulations across scales represent a grand-challenge in computational physics and set a ceaseless quest for novel and efficient computational methods.
The main computational challenge is represented by the spatial and time scale separations involved in a plethora of complex systems, from biology to soft matter, which can be addressed by developing efficient hybrid approaches able to bridge these spatiotemporal gaps.
These methods, generally referred to as multiscale approaches, are based on composite computational schemes 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 hybrid approaches have been proposed to investigate a vast spectrum of physics problems, including soft matter and molecular 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 where the 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 macroscopic hydrodynamics also in the presence of multilevel grids (one and two levels), while retaining the correct thermal fluctuations, 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 fully mesoscale approaches for complex flows across scales, designed in such a way as to capture the physics at the smallest scales, whenever a continuum description alone falls short of providing the correct physical information due to a lack of resolution and thermal fluctuations.

Lattice Boltzmann method
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 and collide on each node, relaxing towards a Maxwell-Boltzmann equilibrium. The LB equation reads as follows: 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 lattice directions (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 the lattice 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 thermodynamic equilibrium f eq i , 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, D2Q9 lattice) has been employed, equipped with a set of second-order expansion of the Maxwell-Boltzmann distribution function, which read as: In the above equation, w i are the weights of the discrete equilibrium distribution functions, c 2 s = i w i c 2 ix = 1/3 the squared lattice speed of sound and u the fluid velocity. The relaxation parameter τ in equation (1) is linked to the fluid kinematic viscosity via the linear relation ν = c 2 s (τ − ∆t/2) [31]. The local (time-space) hydrodynamic moments, up to the order supported by the lattice in use are the density (ρ), linear momentum (ρ u) and momentum flux tensor (Π) [6] and are evaluated by computing the statistical moments of the discrete set of distribution as s I), being I the identity matrix.

Multiparticle collision dynamics with Andersen Thermostat
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 3 − 10 5 in two dimensions [12,14]. During each time step, the system evolves via the subsequent application of a streaming and a collision step. In the streaming step, the particles positions are updated via a forward Euler step: being r k the vector position and v k the velocity of the k − th particle and δt the value of the MPCD discrete time 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-mass momentum. 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 performed within each cell, by rotating the velocity, υ k , of each k-th particle with respect to the velocity of the cell center of mass, υ cm , of all particles in the cell: In the above equation, R is the rotation matrix, which rotates the particle velocities of a given cell by an angle ±α with uniform probability The kinematic viscosity of the SRD fluid [18] which, in two dimensions, reads as follows: being 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 Andersen thermostat (MPCD-AT) [12,26,2].
In the MPCD-AT, new relative velocities are generated during each computational step in each cell and the collision step writes as [26]: where v ran k are random numbers sampled from a Gaussian distribution with variance k B T /m and N c is the actual number of particles in the collision cell. In the above equation, the sum runs over all particles in a given cell.
In the case of the MPCD with Andersen Thermostat, the kinematic viscosity is given by [12]: In this work we employed an MPCD with the Andersen thermostat since, although computationally more expensive than SRD, it features shorter relaxation times. This implies smaller number of time steps required for transport coefficients to reach their asymptotic values [16]. As a consequence the restricition on maximum number of particles per cell (5 − 20), which commonly applies for the SRD, does not apply to the MPCD-AT, where the relaxation times scale as (ln N ) −1 [12].
We wish to point out that the hydrodynamic limit of the multiparticle collision dynamics approach corresponds to the fluctuating Navier Stokes equation for athermal and weakly-compressible fluids. Indeed, it has been previously shown that the MPCD defines a discrete-time dynamics which has been demonstrated to reproduce 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 incorporates 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 thermal transport 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 Boltzmann and k B T in the MPCD) but the MPCD allows the correct incorporation of thermal fluctuations while in the LB, 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].

Coupling Procedure
In this subsection, we detail the coupling strategy employed to exchange the hydrodynamic information between the lattice Boltzmann and the multiparticle collision dynamics.
First, we start from the simplest coupling case, namely ∆x LB = ∆x M P CD = a. In Figure 1, we report a sketch of the hybrid simulation domain. The right region is the pure LB region (green area), while the left one represents the MPCD region (light blue area); the red area in the middle is the handshaking region, where the particles 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 two approaches [13]. The LB runs across the entire domain, the information exchange between the two models is implemented 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, consequently, 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 has proved 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 streaming and collision process ).
2) Information exchange step. In the coupling region the particles velocities are re-generated by drawing them 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 MPCD region.
For the sake of clarity, a pseudo-code is reported in Algorithm 1.

MPCD-AT LB Coupling
Dx a=Dx call LB collision 5: call LB streaming 6: 7: Information exchange step: k B T m 9: 10: 11: call MPCD streaming + boundary conditions 12: call MPCD collision We 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 on a finer grid with a = 1/2 or a = 1/4 (see fig. 2 for a visual sketch of the multigrid hybrid domain with a = 1/2).

MPCD step
The implementation follows the same procedure outlined in Algorithm 1 with an additional intermediate interpolation step between the LB and the MPCD steps, where the linear momentum is interpolated from the coarse to the fine grid (see Algorithm 2).
In the case of the refined simulation the integration time has been chosen so as to avoid sub-cycling in the MPCD region. This is possible since the MPCD allows to set both the integration time step and the grid discretization independently. Thus, in the case a = 1/2 the MPCD, integration time step is set to ∆t = 2, in such 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 the overlapping 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, as detailed 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 be more appropriate.

A note on the computational performances
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 ∼ 0.1µs/particle/step, comparable to the time needed to update a lattice unit in a single component LB code, 10 M LU P S (i.e. Mega Lattice Updates 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 particles and, in our simulations, it varies between 0.01s/step ÷ 0.1s/step with the particle densities ranging between 60 ÷ 600 particles per bin.
It is clear that an efficient parallelization able to exploit the computational capabilities of the latest supercomputers 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 Xeon Platinum 8176 CPU based on the Skylake microarchitecture (2 sockets each one with 28 cores with an installed DRAM of 720 GB ). The code was compiled by Intel Fortran Compiler version 18.0.2 with the recent AVX-512 instruction set supported on Skylake architectures. As per the data structure employed, we opted for a Array of Structure for both the LB and MPCD impolementation (see [30]). A parallel version (openMP-MPI) of the hybrid approach is currently under development.

Coupled LB-MPCD-AT: center line coupling and near-wall coupling with grid refinement
We start by testing the coupling procedure for the case of flow between two parallel plates, the main results being shown in figure 3.
As reported above, the LB runs over a grid covering the entire fluid domain with the MPCD restricted to the lower half of the fluid domain, the depth of the coupling zone extending over two grid nodes, as denoted by the 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÷1000 ((a) to (d)).
The grids share the same spatial discretization and the same kinematic viscosity, ν = 0.167 lu 2 /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 statistical fluctuations associated to the flow field decrease accordingly, as shown in panel (e) of figure 3 which reports the LB and the MPCD velocity profiles.
In the same plot,a comparison between the continuous (MPCD) and discrete (LB) velocity distribution functions at the steady state is reported (inset of figure 3(e)).
The inset displays the time-space averaged velocity distribution. The MPCD velocity distribution (2d 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 Mach number increases).
This departure from the Maxwell Boltzmann behaviour is due to the fact that the set of lattice distributions is not allowed to shift, as instead occurs in the continuous case, but it occurs via a positive bias in the co-flowing distributions. For this reason, the larger the non-equilibrium (represented herein by the large values of the mean channel 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/2 and L/4, being 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 performing simulations 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 a 40 × 14 bins grid.
2) In the multi-resolution implementation, the LB runs on a 20 × 30 nodes grid ( halved with respect to the previous case), and is coupled near the wall to the MPCD which runs on 39 × 13 grid, with a twice finer discretization (a = 1/2) with respect to the LB grid ( figure 4(b)). The first case is used as a a benchmark to test the ability of the coupling procedure to correctly reproduce the 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 quantitatively 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-resolution grids.

Plug Flow in a microchannel with a striction
The multi-resolution algorithm was then employed to simulate the plug flow in a channel with a micrometric striction.
As an example, such geometries, embedded in more complex microfluidic chips, are widely used for biomicrofluidics applications to perform deformability measurements on red blood cells [27]. The domain geometry is reported in figure 5(a). The fluid evolves under the effect of a constant body force ( g) directed from left to right. In the LB simulation, the domain is periodic along the flow direction, while no-slip boundary conditions are 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 hitting with the top and bottom walls while, on the right(outlet) and left (inlet) boundaries, particles which would escape from the MPCD domain are periodically re-injected at the opposite side, in order to exactly conserve the mass within the particle domain.
As shown by the colorbar in figure 5, the maximum flow speed attained within the channel is ∼ 5 · 10 −2 (in lattice unit per time step). The channel is discretized with 10 lattice nodes and the fluid viscosity (in lattice units) is ν = 0.167, corresponding to a maximum Reynolds number within the channel Re = U h ν ∼ 3. By considering a water flow within a 10µm height channel, with a maximum speed of ∼ 10cm/s, typical values 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 overlapping zone extends, as before, over two grid nodes. Within the coupling regions, the particle velocities are initialized by employing the procedure highlighted in Algorithm 1. Figures 5 (b) and (c) report a visual comparison between the flow fields within the striction obtained with the full LB simulation and with the hybrid approach, while the plot below reports the velocity profiles taken at three different sections of the microchannel. These results confirm that the hybrid approach is able to capture both 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 effectiveness of 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 simulation on 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 flow fields ((a) and (b)) and the velocity profiles along the red-dashed vertical sections of the striction. In order to compare 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 the physics of fluid within the micron-sized channel.
Moreover, as stated in the introduction, the hybrid mesoscale approach not only allows to locally refine the grid to obtain more accurate solutions of the hydrodynamic fields, but also naturally preserves the fluctuating statistics, typical of confined fluids the micro-and nano-scales, which is absent in the non-fluctuating LB models [1,9].

Shear flow in a channel with a micro-cavity
As a final case, we present a shear flow in a periodic channel with a microcavity placed at the bottom. This kind 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 and entrapment, 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 a momemtum exchange boundary condition to the top boundary [4]. Periodic boundary conditions are applied along the flow directions, while no-slip condition is implemented on the bottom boundary and on the cavity walls.
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. The kinematic viscosity is set to ν = 0.167 and the shear velocity imposed at the uppermost boundary is U = 0.4 in lattice units.
ii) The cell density in the MPCD simulations is set to N = 1000. iii) As per the boundary conditions employed in the MPCD, the bottom and side walls are no slip, while the uppermost boundary is treated as a free slip wall as the velocity is fixed during the information exchange between the LB and MPCD grids.
It is again of interest to compute the cavity Reynolds number Re = U L ν ∼ 6, typical of micro-cavities under shear [32,17]. The MPCD run in the region inside the dashed rectangular perimeter (see figure 7)and the coupling region coincides with the separation line between the microcavity and the main channel. As a first case, we run the coupled approach by using the same grid discretization both in the LB and in the MPCD regions. The first set of results is shown in figure 8. As shown in panel (a) and (b), the main vortex inside the microcavity is accurately captured both in the LB and in the MPCD simulation. Further, the vortex positions and sizes in both the representaion also show a good agreement.
We then plotted the horizontal component and the magnitude of the velocity ( √ u 2 + v 2 ) along the height (h) of the cavity (panels (c) and (d)). The agreement is evident and confirmes the smooth reconnection between the MPCD and LB solutions within the coupling region. Furthermore,the resolution in the MPCD region has been increased by a factor two and four(a = 1/2 a = 1/4 ). The main results are reported in figures 9 and 10, which report, respectively, a comparison between the LB and MPCD solutions and the flow fields within the cavity for different values of the grid size. Here too, the agreement between the LB and the MPCD solution is remarkable, as is also confirmed by the subplot in panel (c), showing the LB and MPCD velocity profiles, taken along the height of the cavity, in close match with each other. Also in this case, the multi-resolution approach proves capable of providing the correct macroscopic hydrodynamic information with thermal fluctuations directly incorporated in the MPCD solver.

Conclusions
In this work a concurrent lattice Boltzmann-multiparticle collision dynamics multiscale/multigrid approach has been developed and validated against the classical laminar flow through parallel plates and employed to simulate fluid flows in confined geometries at the micro-scale.
In particular, for the Poiseuille flow across two parallel plates, several cases have been investigated, namely  1) hybrid LB-MPCD with the coupling region positioned at the center of the channel 2) Coupled LB-MPCD with one level near-wall grid refinement Furthermore, the multi-resolution approach has been used to simulate a plug flow in a microchannel with a micrometric 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 simulate fluctuating Navier-Stokes fluid dynamics even in the presence of complex geometries and velocity gradients near the 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 into account only in small portions of the fluid domain. In this respect, an example is represented by the flow in the microcavity. Indeed, within the microcavity, both hydrodynamic interactions and thermal fluctuations play an important role while outside, in the sheared channel, fluctuations are not essential and a lattice Boltzmann approach can safely be applied.
Further developments of the hybrid approach, currently ongoing, will aim at the simulation of complex fluid phenomena in confined systems, also in the presence of hydrodynamic interactions with nanoparticles, and within thin films at the fluid-fluid interface, wherein increased resolution and thermal fluctuations are both needed to capture the relevant microscale physics.

Data Availability
The code is available upon request.