## Abstract

Hybrid particle–continuum computational frameworks permit the simulation of gas flows by locally adjusting the resolution to the degree of non-equilibrium displayed by the flow in different regions of space and time. In this work, we present a new scheme that couples the direct simulation Monte Carlo (DSMC) with the lattice Boltzmann (LB) method in the limit of isothermal flows. The former handles strong non-equilibrium effects, as they typically occur in the vicinity of solid boundaries, whereas the latter is in charge of the bulk flow, where non-equilibrium can be dealt with perturbatively, i.e. according to Navier–Stokes hydrodynamics. The proposed concurrent multiscale method is applied to the dilute gas Couette flow, showing major computational gains when compared with the full DSMC scenarios. In addition, it is shown that the coupling with LB in the bulk flow can speed up the DSMC treatment of the Knudsen layer with respect to the full DSMC case. In other words, LB acts as a DSMC accelerator.

This article is part of the themed issue ‘Multiscale modelling at the physics–chemistry–biology interface’.

## 1. Introduction

The need for the development of fast and accurate numerical methods to simulate gas flows characterized by a broad range of time and space scales represents a general challenge to statistical physics modelling. Moreover, despite the success demonstrated by computer simulations at the microscopic level as a powerful tool of investigation, the computational requirements of the microscopic methods, such as molecular dynamics (liquids) and direct simulation Monte Carlo (DSMC; gases), are still prohibitive for many practical applications. It is therefore evident that, for such kind of problems, a multiscale approach, consisting of the use of different models for each relevant scale, is needed. In order to link large scales to local microscopic properties, a variety of scale bridging techniques have been developed, which can be broadly classified as either *sequential* or *concurrent*. In the former class of methods, the continuum treatment conveys all microscopic details into appropriate constitutive relations or laws [1–4], and systems are then simulated separately at the different levels of resolution. In the latter class, which is the one relevant to this work, systems are analysed in a single simulation at different levels of space and time resolution. Within this approach, which has been implemented to investigate a vast spectrum of physics problems, including hard matter [5–7], soft matter and molecular fluids [8–11], as well as dilute and dense fluid flow problems [12–17] and fluctuating hydrodynamics [18,19], a small region of the domain is analysed at a finer scale level, whereas the remaining part is treated on a coarser and computationally less demanding level; a hand-shaking region and appropriate exchange of information schemes then enable the coupling and communication across the different parts.

When investigating multiscale dilute gas flow problems, the most common approach is to adopt a solver of the continuum Navier–Stokes equations [12] or Euler equations [13] for the part of the domain where the flow exhibits a small deviation from equilibrium, and rely on kinetic methods, typically the DSMC, wherever substantial non-equilibrium effects are present. While this coupling has led to a number of interesting results, it is plausible to expect that the replacement of the Navier–Stokes equations with a lattice kinetic solver, such as lattice Boltzmann (LB), may prove convenient in several respects. In particular, as well documented in [20], the kinetic nature of LB is reflected in a complete disentangling between non-locality and nonlinearity, i.e. kinetic transport along molecular trajectories (free streaming) is linear, as opposed to nonlinear hydrodynamic advection along fluid material lines. Moreover, diffusion emerges from relaxation to local equilibria, with no requirement for second-order spatial derivatives. Both features may offer significant advantages in practical implementations, especially in connection with solid boundaries. In this work, we depart from the models proposed in the literature, as the flow at the continuum level is simulated with the efficient LB method in the Bhatnagar–Gross–Krook (BGK) single-relaxation-time formulation [21]. Because LB and DSMC are both widely documented in the literature, the reader should refer to [22,23], which represent the reference for the DSMC and LB methods, respectively. Although both DSMC and LB bear strong links to the Boltzmann kinetic equation, they show major differences in purpose and mathematical structure. Indeed, while DSMC is poised to solve the actual Boltzmann equation by retaining microscopic realism, the LB was initially devised to retain only the symmetry/conservation constraints required to reproduce the hydrodynamic equations in the low-Knudsen-number hydrodynamic limit. Thus, LB could be regarded as a Navier–Stokes solver in kinetic vests. It is only recently that higher-order versions of the LB method (LBM) have been developed, which prove capable of capturing (some) non-equilibrium effects beyond the hydrodynamic picture [24–26], but such higher-order versions are not relevant to this work. The main feature that distinguishes the LB from the DSMC is the coarse graining operated at the level of the microscopic velocity space, as the degrees of freedom of the velocity space itself are largely reduced. In fact, while, in LB, probability distributions at each lattice site **x** are allowed to propagate only along a finite set of directions with an assigned speed *ξ*_{a}, in DSMC particles are grid-free.

## 2. Coupling schemes

The most important steps of any hybrid model are the definition of accurate protocols to exchange information between the different solvers at the overlapping region and the decomposition of the computational domain into the particle, continuum and overlap domains. For what concerns the first issue, as LB involves the solution of kinetic equations for the single-particle distribution function evaluated at a set of discrete microscopic velocities, *f*_{a}(**x**,*t*), and not directly the hydrodynamic moments or their fluxes as in classical Navier–Stokes equations solvers, we developed specific mapping schemes that allow information to be passed from the DSMC to the LB domains and vice versa, correctly transferring also non-equilibrium information. In [27], the implementation details are extensively presented. Here, for completeness, we briefly recall the most important features.

The basic idea of the mapping schemes is founded on Grad’s moments method [28] and on Gauss–Hermite quadratures [29]. According to Grad’s formalism, the single-particle distribution function *f*(**x**,** ξ**,

*t*) can be expressed as a series in the dimensionless Hermite orthonormal polynomials, , in the velocity space

**as 2.1where**

*ξ**ω*(

**) is the weight function associated with the Hermite polynomials and**

*ξ***a**

^{(n)}are the rank-

*n*expansion coefficient tensors, which, in turn, can be identified as the hydrodynamic moments (or a combination of those), e.g.

**a**

^{(0)}=

*ρ*,

**a**

^{(1)}=

*ρ*

**u**,

**a**

^{(2)}=

**P**+

*ρ*(

**uu**−

*c*

^{2}

_{s}

**), where**

*δ**ρ*is the fluid density,

**u**is the fluid velocity,

**P**is the momentum flux tensor,

*c*

_{s}is the speed of sound and

**is the Kronecker delta. In equation (2.1), we took advantage of the orthonormality of the Hermite polynomials to truncate the series at the order**

*δ**N*: the complete and the truncated distributions have the same leading

*N*velocity moments [24].

Starting from this common ground, the two mapping schemes can be established (see also figure 1*a* for a schematic representation of the schemes). In particular, the *projection step* (DSMC→LBM) allows one to project the DSMC hydrodynamic moments (fine level of description) onto the LBM discrete distributions (coarse level of description). In correspondence of the overlapping region, where both descriptions are accessible, the coefficients **a**^{(n)} in equation (2.1) are computed using the DSMC moments. In fact, because the distribution *f*^{N}(**x**,** ξ**,

*t*) is completely and uniquely determined by its values at a set of discrete velocities, one has 2.2where

*w*

_{a}and

*ξ*_{a}are the weights and abscissae of a Gauss quadrature of algebraic precision of degree ≥2

*N*and

*q*is the total number of discrete velocities. From inspection of equation (2.2) and taking into account the classical procedure in which moments are computed in LBM, we immediately see that the non-equilibrium discrete distributions to be supplemented to the LBM at the coupling nodes are the scaled values of the continuous distribution function evaluated at the quadrature abscissae, 2.3It is now possible to define the order

*N*at which the distribution in equation (2.1) is truncated, as this is directly related to the algebraic precision of the quadrature. As in this work we present results for the D3Q19 lattice only, we set

*N*=2, because the proposed quadrature is able to correctly compute the distribution moments up to the momentum flux tensor,

**P**. It can be shown [30,31] that the resulting solution generated by the LBM is at an isothermal Navier–Stokes level of description.

The inverse *reconstruction step* (LBM→DSMC) allows one to reconstruct from the LBM discrete distributions (coarse level), the continuous, truncated distribution function (fine level), from which, in turn, the velocities of the DSMC particles can be sampled. At the coupling nodes, the LBM non-equilibrium populations *f*_{a}(**x**,*t*) are used to determine the coefficients of the expansion in equation (2.1),
2.4(cf. equations (2.2) and (2.3)). From the knowledge of the coefficients, the truncated distribution *f*^{N}(**x**,** ξ**,

*t*) is also known. An acceptance/rejection algorithm is then used to generate the velocities of the particles, which, eventually, will evolve according to the DSMC rules. It should be noted that the zeroth- and first-order coefficients in equation (2.4) are routinely available in the LBM implementation, because they are used to evaluate the equilibrium distribution functions. In this work, therefore, as only the second-order coefficient terms have to be computed on purpose, the computational overhead for this coupling step turns out to be very limited.

The domain decomposition is strictly related to the information exchange strategy implemented at the interface between the different methods. In the proposed scheme, to enforce a stronger coupling between the two methods, we extend the LBM to cover the whole simulation domain, whereas the DSMC is confined only within some regions of the domain. In this way, a two-level grid is obtained (figure 1*b*). A few basic subroutines define the implementation of the hybrid simulation adopting the illustrated projection and reconstruction steps. During the time integration of the LBM equations, the discrete distributions *f*_{a}(**x**,*t*) are advanced on a coarse time level increment, Δ*t*_{LBM}, over the whole domain, including the part that lays underneath the DSMC region. Then, the DSMC domain is integrated, and it is advanced to the new time on the LBM level through a sequence of fine level time steps, Δ*t*_{DSMC}. The ratio between the two time steps, depending on the LBM lattice spacing and the DSMC particle characteristic collision time, defines the number of subcycles the DSMC should run for each LBM time step. The LBM solution information in proximity of the fine level grid is passed to the DSMC via the imposition of proper boundary conditions on the particle simulation domain. At the beginning of each DSMC integration step, in fact, particles are generated within the buffer region by sampling Grad’s distribution built according to the reconstruction step. The particles, in both the DSMC and the buffer regions, are allowed to stream for one fine level time step. If a particle from the DSMC region crosses the interface, that particle will contribute to the fluxes of mass and momentum for the LBM discrete distributions corresponding to the node where the crossing occurred, so that the populations will be corrected as
2.5where are the populations before the application of the correction, *R* is the ratio between the LBM and DSMC time-step durations, and the correction contribution, , can be evaluated from equation (2.1) for *n*≥1 (only the fluxes are exchanged). After the moving of the particles, those still in the buffer region, or that moved into it, are discarded and collisions among the remaining particles are performed. Finally, the DSMC solution in the fine level region is transferred to the LBM in the overlap region (red area in figure 1*b*) using the outlined projection step, overwriting the LBM discrete distributions for the lattice nodes in this area with the discrete distributions obtained from DSMC. A new integration of the LBM equations can then be performed.

## 3. Numerical results

To test the proposed hybrid method, we consider the steady isothermal Couette flow problem. The flow is driven by the walls, located at *x*/*H*=0 and *x*/*H*=1, moving with a velocity *u*_{wall} such that *Ma*_{wall}=0.1, where *H* is the channel height. The isothermal condition is guaranteed by the negligible viscous dissipation generated by the moving walls. The DSMC particles are treated as argon-like hard spheres (diameter *d*_{ref}=3.66×10^{−10} m; mass *m*=6.63×10^{−26} kg). The reference temperature *T*_{ref} and the speed of sound, *c*_{s}, are 273 K and 308 m s^{−1}, respectively. The reference density is set so that Knudsen number *Kn*=0.10, based on the channel height equal to 1 mm. At the reference density, 100 particles are initially placed in each DSMC cell, which, in turn, is of length 0.1λ_{0}, λ_{0} being the mean free path. As shown in figure 2, the flow domain is decomposed in such a way that DSMC extends over a region in proximity of the walls, where stronger non-equilibrium effects are expected to occur, equal to *α* and the pure LBM domain extends over the remaining part of the domain. For the results of the simulations of this work, each buffer region occupies a width equal to the length of a DSMC cell. To enforce on the LBM domain the flow *Kn* number, we set the relaxation time *τ* according to the relation [32,33], where for the D3Q19 lattice and *N*_{H} is the number of nodes along the channel height. DSMC and LBM grids are then set, so that the centres of DSMC cells and the LBM nodes, in the overlap regions, are located at the same positions **x**. In order for the LBM to accurately simulate such flow, it is fundamental to implement proper wall boundary conditions, namely *diffuse reflection* of populations impinging into the walls [34], and the *regularization procedure* [35,36], which, projecting the non-equilibrium discrete distributions onto the same Hermite polynomial basis as used to project the equilibrium populations, allows one to get rid of the moments not sufficiently supported by the lattice. Based on the test conditions here investigated, five fine level time integration steps for the DSMC are performed between two consecutive coarse level time integration steps for the LBM.

Taking into account the characteristics of the tested flow, a few additional considerations about the flux exchanges at the DSMC–LBM interface as prescribed by equation (2.5) are due. One disadvantage of imposing the fluxes across the coupling interface, instead of the related densities, resides in the fact that fluxes show larger statistical noise [37]. As a means to increase the signal-to-noise ratio of the fluxes exchanged at the DSMC–LBM interface, besides using the cumulative averages over time, justified by the steady-state nature of the investigated flow, we choose to evaluate them by averaging over the cell volumes, instead of measuring across the interface surface, see [12,38]. Moreover, it has to be noted that, for the tested case, there is no net mass flux across the interface. Nonetheless, we also impose and test the mass flux conservation constraint across the DSMC–LBM interface. We also avoid having to implement fitting of the mean moment profiles to damp out fluctuations faster, as done in [3,4], because such practice cannot be easily extended to complex flows.

The proposed flow problem is at *Kn*=0.10. While the rarefaction effects for the tested conditions can be treated within the bulk of the flow by implementing appropriate boundary conditions for the coarse level solver, non-continuum effects produced by the flow gradients in proximity of the walls cannot be captured by the classical D3Q19 model, so a microscopic description is needed. The flow considered here, therefore, appears appropriate to validate the proposed hybrid model, as a separation of scales is naturally present. In the bulk of the flow, in fact, the hydrodynamic variables vary substantially over the macroscopic scale, in turn dictated by the domain geometry; close to the walls, especially within the Knudsen layer, those changes take place at a length scale much closer to the molecular one. In the literature, several analytical solutions for such uniform shear stress flow have been proposed both for the case of Maxwell molecules, based on the moment-hierarchy method [39,40], and for the BGK collision term formulation [41,42]. In particular, for the latter case, which is of interest here, it was found that, within the bulk of the flow, for distances greater than at least a mean free path from either wall, a *normal solution*, i.e. a solution constructed by the Chapman–Enskog method, exists, at Navier–Stokes order, for a *Kn* based on channel height smaller than about 0.2. From a mathematical point of view, and within Grad’s formalism, this is equivalent to saying that, whereas the bulk flow can be accurately described by a single-particle distribution function truncated at the order *N*=2, the flow at the wall requires that the distribution should include terms above the second-order ones. However, a clear connection between the rarefaction and non-equilibrium conditions and the order at which the series in equation (2.1) should be truncated is still missing [43]. From all these considerations, it appears reasonable that the DSMC domain should cover at least one mean free path from the walls, i.e. *α*/2≥λ_{0}.

Figure 3*a* shows the velocity profiles, , averaged over time and normalized with respect to the wall velocity, as obtained by a full DSMC simulation and by a hybrid simulation when the DSMC region extends over a mean free path away from both walls (*α*=0.2) and the buffer region is 0.1λ_{0} wide. A good agreement between the two profiles is found as can be seen from the inset of figure 3*a*, where the relative error, defined as
3.1is plotted. The largest error, in fact, is below 2%, whereas, in the case of a full LBM simulation (orange dots in the inset of figure 3*a*), a 5% error is found to occur within the Knudsen layer. This result is obtained strictly maintaining mass conservation within the system. In figure 3*b,* in fact, the instantaneous number of particles present in the DSMC region is compared with the averaged one, which, in turn, deviates from the initial one by less than 0.05%. More importantly, we also measured the speed-up provided by the hybrid scheme versus full DSMC (*α*=1) as a function of the parameter *α*.

In figure 4, such speed-up is compared with the speed-up one would obtain by just summing the computational times required for the DSMC to run over the *α* fraction of the domain and for LBM over the whole domain. The computational times are recorded when the solution reaches a global error, for both the velocity and shear stress profiles, , below 3% with respect to a reference full DSMC simulation. The global error for the velocity is defined as
3.2and analogously for the shear stress. The reference DSMC simulation has cumulated enough statistics to confine the difference between the maximum and the minimum values of the shear stress below 0.02%. The recorded times are then averaged over several realizations of the tested flow. It is possible to see that a substantial speed-up is achieved by applying the hybrid method. More interestingly, the required computational time, *t*_{hybrid}, is smaller than what would be just the sum of the times required for DSMC to run over the *α* fraction of the domain and for LBM over the whole domain, *t*_{LBM+DSMC}. In particular, for *α*=0.2, more than a sixfold reduction with respect to a full DSMC simulation is achieved. It has also to be noted that the computational overhead, due to the fact that the LBM extends over the whole domain and not just over (1−*α*), is very limited. For the tested conditions, in fact, a full LBM simulation is about 80 times faster than a full DSMC simulation, performed over the whole domain. The measured extra speed-up, defined as *t*_{LBM+DSMC}/*t*_{hybrid}, was found to be 25% for *α*=0.2 and 16% for *α*=0.36. Such computational extra gain may qualitatively be explained by considering that the LB information processed in the bulk provides a higher quality, less noisy, input to the DSMC in the Knudsen layer than bulk DSMC alone. Thus, LB serves as a sort of DSMC accelerator.

Finally, a clarification about the flow domain decomposition technique and its range of applicability is due. This type of hybrid approach is a suitable choice to analyse multiscale problems if the flow features a bulk region for which a constitutive model is known to be accurate, but it may become not an efficient nor an accurate method if such a condition is not fulfilled. This would be the case, for example, of the Couette flow at a larger *Kn* number than the one presented here. For such a condition, it is known that rarefaction effects spread also across the bulk of the flow and the classical Navier–Stokes equations can no longer accurately reproduce the flow properties found in the bulk. To overcome this limitation, however, the kinetic roots of the LBM may help. In fact, whereas LBM was specifically designed to solve the Navier–Stokes equations [30,31], in the last decade various groups have begun to explore the possibility to extend the LB methodologies towards non-equilibrium flows beyond Navier–Stokes hydrodynamics by adopting higher-order lattices [24–26] or introducing further modelling in order to mimic effects of free particle streaming [44]. This is also equivalent to saying that, if a higher-order lattice is adopted, the size of the domain where the DSMC solver is still needed can be significantly reduced, thus further improving the overall computational efficiency of the simulation. In this view, it should be pointed out that the coupling strategy here proposed can be directly applied to any suitable lattice, provided that the discrete speeds are abscissae of a Gauss–Hermite quadrature. The extension to higher-order lattice models of the hybrid scheme is currently under development.

A remark on the range of application of the proposed hybrid method is in order. While the DSMC can correctly predict thermal effects, the LB model here employed is not able to reproduce such effects, thus limiting the current implementation to isothermal flows. To overcome such limitation, essentially stemming from the lack of sufficient isotropy of the lattice, different strategies may be envisaged. In particular, three main approaches have appeared in the literature to extend the LB method to thermal flows: the mixed approach, where the velocity field is solved using a usual LB model and the energy, or temperature, equation is solved by adopting a different numerical method such as finite difference or finite volume method [45]; the double distribution function approach, where a second lattice is adopted for the energy field [46]; and the multispeed approach, where a single lattice fulfilling higher-order isotropy conditions is used [25]. The most natural choice among the proposed methods is the multispeed scheme, as it shares with the present hybrid model the common Hermite expansion approach. The extension to non-isothermal flows, however, will be pursued in future works.

## 4. Conclusion

In conclusion, a novel multiscale hybrid scheme for dilute gas flow simulations relying on the accuracy of the DSMC for the regions of the flow where a microscopic scale description is needed and on the efficiency of the LBM for the parts of the flow where hydrodynamic properties evolve over macroscopic scales has been proposed. The method, which is based on Grad’s formalism and Gauss–Hermite quadratures, enables the flux exchange across the LBM–DSMC interface, conserving both mass and momentum. We applied the hybrid formulation to the steady Couette flow and we found that, alongside a good agreement with a full DSMC solution, up to a sixfold reduction of the computational time cost can be achieved. The application of such a model could also open the way to study macroscopic flows characterized by the presence of local regions where non-equilibrium effects are strong, and therefore the method is not restricted to classical micro-mechanical systems.

## Data accessibility

The datasets supporting this article have been uploaded as part of the electronic supplementary material.

## Authors' contributions

G.D.S. performed the simulations, data analysis and drafted the manuscript. H.J.H.C. conceived the project and helped draft the manuscript. S.S. conceived the project and helped draft the manuscript. F.T. conceived the project, and helped with the coding and drafting the manuscript. All authors read and approved the manuscript.

## Competing interests

The authors declare that there are no competing interests.

## Funding

No funding has been received for this article.

## Acknowledgements

This research is supported by the Dutch Technology Foundation STW, which is part of the Netherlands Organisation for Scientific Research (NWO), and which is partly supported by the Ministry of Economic Affairs. One of the authors (S.S.) was partially supported by the Integrated Mesoscale Architectures for Sustainable Catalysis (IMASC), an Energy Frontier Research Center supported by the US Department of Energy, Office of Science, Basic Energy Sciences under award no. DE-SC0012573.

## Footnotes

One contribution of 17 to a theme issue ‘Multiscale modelling at the physics–chemistry–biology interface’.

- Accepted June 27, 2016.

- © 2016 The Author(s)

Published by the Royal Society. All rights reserved.