Progress in particle-based multiscale and hybrid methods for flow applications

This work focuses on the review of particle-based multiscale and hybrid methods that have surfaced in the field of fluid mechanics over the last 20 years. We consider five established particle methods: molecular dynamics, direct simulation Monte Carlo, lattice Boltzmann method, dissipative particle dynamics and smoothed-particle hydrodynamics. A general description is given on each particle method in conjunction with multiscale and hybrid applications. An analysis on the length scale separation revealed that current multiscale methods only bridge across scales which are of the order of O(102)-O(103)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\mathcal {O}}(10^2) {-} {\mathcal {O}}(10^3)$$\end{document} and that further work on complex geometries and parallel code optimisation is needed to increase the separation. Similarities between methods are highlighted and combinations discussed. Advantages, disadvantages and applications of each particle method have been tabulated as a reference.


3
68 Page 2 of 38 Knudsen numbers; therefore, the continuum approach breaks down.
Particle-based methods, such as molecular dynamics (MD), are derived from atomistic observations and therefore valid on much smaller physical length scales. It is, however, not feasible to employ these approaches on the continuum level as the computational cost becomes prohibitively expensive. Therefore, particle-based hybrid and multiscale methods were actively developed over the last two decades. MD is capable of simulating the assumed correct behaviour of the slip boundary conditions, and it starts to behave as a continuous medium at around 10-20 atomic diameter taken from the wall (Asproulis et al. 2012). For channels exceeding a couple of hundred atomic diameters in their height (Asproulis et al. 2012;Xu and Li 2007), computational efficiency is expected to be high by using multiscale methods. In addition to this, experimental capabilities such as micron-resolution particle image velocimetry (µ-PIV) (Santiago et al. 1998;Tretheway and Meinhart 2002) can still only operate at microscales. Advances towards smaller scales are hampered by the diffraction limit, noise in the particle image, interaction between fluid and seed particles, and the effects of the Brownian motion. Multiscale approaches are situated between pure particle simulations and experiments and able to obtain results in an efficient and confident manner; however, relatively few studies have been published concerning this overlapping region.
Micro-and nanofluidics became important in recent years as can be seen by the advent of MEMS, µ-TAS and lab-on-a-Chip devices (Gad-el Hak 2001;Ho and Tai 1998;Manz et al. 1993). All of them are operating at length scales in which the continuum approach may not be valid. Squires and Quake (2005) investigated the flow physics at these small scales in terms of dimensionless numbers. They investigated the inertial effects (Reynolds number), advection and diffusion (Péclet number), interfacial tension (Capillary number), elastic effects occurring in deformable objects such as polymers (Deborah, Weissenberg and Elasticity number), density-driven flows (Grashof and Rayleigh number) and the validity of the continuum model (Knudsen number) concluding that the behaviour of the fluid flow at microscales differs due to the increased surface-tovolume ratio. More information can be found in the papers of Bayraktar and Pidugu (2006) focusing on flow physics in microchannels, and Gravesen et al. (1993) investigated micropumps, microvalves and microsensors. Koo and Kleinstreuer (2003) categorised the literature on microchannels into "flow instabilities", "viscous changes" and "no changes compared to macroscale". Understanding the fundamental different behaviours of the fluid flow at microscales is mandatory to appreciate the need for multiscale methods. The combination of two or more fluid dynamic approaches may accurately describe the fluid flow behaviour at the atomistic level with a particle-based approach while simulating the bulk of the fluid domain with an efficient Navier-Stokes solver.
In this paper, we present practical examples that set the scene for multiscale methods and we discuss the main findings in recent years related to several particle-based multiscale and hybrid methods including molecular dynamics (MD), direct simulation Monte Carlo (DSMC), lattice Boltzmann method (LBM), dissipative particle dynamics (DPD) and smoothed-particle hydrodynamics (SPH). Each method has been discussed separately in its own section where a brief overview has been given on the method and its governing equations. Each section presents the current state-of-the-art research followed by an interim conclusion. We summarise our findings on particle-based multiscale and hybrid methods in the last section.

Molecular dynamics method
Molecular dynamics models the movement of atoms at the atomistic scale. We invoke Newton's second law for every atom i directly as The inter-atomic forces are replaced by the derivative of the inter-atomic potential. The most used one is the 12-6 Lennard-Jones (LJ) potential where ǫ is the depth of the potential well, σ is the distance at which the potential between two atoms is zero and r ij is the distance between two atoms. In theory, we have to sum over all atoms for each atom to obtain V ij which is computationally expensive and proportional to O(N 2 ). We define a cut-off distance r cut after which we neglect the inter-atomic interactions. This is possible because the potential [Eq. (2)] asymptotically approaches zero for large values of r ij and therefore, the computational time reduces to O(N log N). From kinetic gas theory, we obtain the equipartition theorem for single atoms as where v i is the velocity of the atom, k B is the Boltzmann constant and T is the temperature. While we can conserve the overall momentum of a system by using Eq. (1), we cannot control the temperature. Therefore, it is common practice to rescale the velocity according to Eq. (3) or to (2) V (r ij ) = 4ǫ σ r ij Page 3 of 38 68 use a thermostat. One widely used thermostat is that of Berendsen et al. (1984) where the equation of motion, Eq. (1), is coupled to a heat bath as where T 0 is the target temperature, T is the current (computed) temperature and γ is the strength at which the equation should relax towards T 0 . The atoms are initialised with a random velocity at a temperature T obeying the Maxwell-Boltzmann distribution as Boundary conditions are not straight forward in MD simulations and require special attention. Where possible, periodic boundary conditions are applied. Particles can leave and enter the domain seamlessly. The problem arises when dealing with open boundary conditions. Particles may be inserted or removed easily; however, the force evaluation near the boundary results in an incorrect interaction potential due to the missing particles beyond the boundaries which may propagate into the domain. One remedy may be to sample the force near the boundary via a second simulation where periodic boundaries are applied, see, for example, Steijl and Barakos (2012). This, however, requires a second set of simulation data for describing correct boundary conditions. Solid boundaries are less problematic to impose. Atoms may be placed rigidly in a lattice structure to model a wall. Asproulis and Drikakis (2011) pointed out the danger of using non-physical, high mass values to effectively freeze atoms into position in that it produces incorrect slip length behaviour. Asproulis and Drikakis (2010) showed that using spring potentials for atoms on solid boundaries produces more realistic values for the slip length. See also Thompson and Troian (1997) for a description on solid boundaries and Delgado-Buscalioni et al. (2015) for an overview on open boundary conditions.

Review on hybrid and multiscale molecular dynamics methods
In 1995, O'Connell and Thompson (1995) conducted the first coupled computation of MD and continuum mechanics. Ever since there has been a broad interest in this field as indicated by the growing number of publications. They investigated the flow in a channel and split it into a continuum (Navier-Stokes) and atomistic (MD) region, separated by a hybrid solution interface (HSI) running parallel to the solid boundaries across the channel. The HSI is a buffer layer in which both descriptions are valid and it provides boundary conditions for each other. In this way, information from one description is passed to the other and an information exchange can take place. The buffer region introduced by the HSI is needed to obtain smooth profiles for the primitive variables across the interface. The MD and continuum region were coupled by using constrained dynamics, and the exchange of variables was done by state, i.e. velocity and density were imposed directly at the HSI. Their test case was the start-up Couette flow, for which good agreement between MD and continuum could be observed. Hadjiconstantinou and Patera (1997) investigated the flow around a square cylinder inside a channel and coupled non-equilibrium molecular dynamics (NEMD) with a spectral element solver in the wake of the cylinder. The solution at the interface was iteratively obtained by the Schwarz alternating method. Since the coupling region was chosen away from the wall, both NEMD and Navier-Stokes solutions were valid and could be compared against a full Navier-Stokes simulation. Good agreement was achieved between the full continuum and hybrid computation although the accuracy was limited due to statistical fluctuations, boundary condition imposition by NEMD and mismatch of the transport coefficient in the two models. Following up on their research, Hadjiconstantinou (1999) investigated the moving contact-line problem (two immiscible but otherwise identical fluids) in a microchannel for low Reynolds numbers. The continuum domain was placed at the channel centre, and the walls were resolved by MD. The initial velocity distribution was obtained via a full continuum solution, and therefore, only a few Schwartz iterations were necessary to converge the multiscale approach to its final solution. Comparison between a full MD solution showed a similar behaviour as in their previous study (Hadjiconstantinou and Patera 1997), where the overall trend was matched with fluctuations noticeably present. Flekkøy et al. (2000) introduced a different approach to couple the atomistic and continuum region by imposing fluxes in the HSI. The Navier-Stokes equations are written in the following form where we define and q as follows (6) ∂ρ ∂t + ∇ · (ρu) = 0, ∂ρe ∂t + ∇(ρeu + · u + q) = 0, = ρuu + P − ∇ · u − µ (∇u) + (∇u) T − ∇ · u , where µ and are the dynamic and bulk viscosities, T is the temperature and k c is the thermal conductivity. We can compare the differences of coupling by states and fluxes in Tables 1 and 2, where we have defined ρ, u and e as the macroscopic density, velocity and total energy, A as the cross-sectional area across which the flux is evaluated, V c as the volume of the cell inside the HSI, n the normal vector on A, m as the mass of an atom, s as the number of added/ removed particles, v ′ i and ε ′ as the velocity and energy of the added/removed particles, F ij as the force between two atoms i and j (and F i as the force acting on a singular atom i ), J Q as the molecular energy flux, E p,i as the potential energy of atom i and 〈 〉 as a time average. and q have been defined in Eqs. (9) and (10), respectively. We can see from Table 1 that density, momentum and energy are directly imposed on the atoms from the continuum when coupling by states is used. The macroscopic quantities for density, momentum and energy are obtained for each computational cell individually in the HSI by averaging over the atoms inside the cell c. Fluxes have to be calculated first before they can be imposed on the continuum or atomistic side, as shown in Table 2. To conserve the overall (10) q = −k c ∇T , mass, momentum and energy, it is necessary to remove or add atoms in the HSI. Since particles are free to move and interact inside computational cells, their energy level will change locally inside the cell. The energy level is associated with the temperature of that cell and so by randomly inserting and deleting atoms, the temperature may change and the overall conservation of energy is not guaranteed. Delgado-Buscalioni and Coveney (2003a, b) introduced the USHER algorithm to remove this shortcoming. In their algorithm, particles are introduced into a cell at a predefined energy level which has been determined by a steepest-descent approach so as to conserve the total energy. The USHER algorithm is equally applicable to coupling by states and fluxes. Both approaches were successfully applied, see O'Connell and Thompson (1995), Hadjiconstantinou and Patera (1997), Hadjiconstantinou (1999), Nie et al. (2004), Wang and He (2007) for state coupling and Flekkøy et al. (2000), Delgado-Buscalioni andCoveney (2003a, 2004) and Delgado-Buscalioni et al. (2005a, b) for flux coupling. Delgado-Buscalioni et al. (2005b) demonstrated the need to incorporate the fluctuating component of the fluxes as they can impact the overall solution of the coupling scheme. They modified the boundary conditions for the continuum which did not influenced the flux conservation but allowed for fluctuations during the exchange. Applied to the oscillatory shear flow, the influence of the fluctuations was clearly seen and qualitatively good results achieved. Markesteijn et al. (2014) coupled the Landau-Lifhitz fluctuating hydrodynamic (LL-FH) equations with MD motivated by the fact that a peptide immersed in water showed a strong correlation between its deformation (dihedral angle) and the density fluctuation. Tested in two simulations (rectangular domain with separated LL-FH and MD domain), one with zero mean and one with a drift velocity on the continuum side, while in both cases the MD domain was initialised with no net momentum, the MD solution was approaching the continuum solution while fluctuations were preserved. O'Connell and Thompson used constrained dynamics to exchange data at the HSI. They introduced a coupling strength parameter ξ which in their study was set to ξ = 0.01. They observed that higher values negatively influenced the results although increasing ξ would cause lower computational times. Nie et al. (2004) were able to use a value of unity for the coupling strength in conjunction with an external force applied to the particles in the HSI to prevent them from escaping. Wang and He (2007) derived an equation for ξ allowing them to dynamically update the value of ξ, while Kamali and Kharazmi (2013) used yet another approach and imposed an arithmetic relation on ξ driven by the molecular time step. The simulation was initialised with ξ(t MD = 0) = 0.1 and increased ξ via a thirdorder polynomial until ξ = 1 was reached at a user-defined time step. They argued that fluctuations at the beginning of the exchange were responsible for the divergence. Using a low coupling strength parameter at the beginning suppressed possible divergence and the results converged as the coupling parameter was increased.
One of the challenges in coupling the continuum with an atomistic description is partially due to the imposition of boundary conditions from the continuum onto the atomistic level. While molecular data are easily averaged and imposed on the continuum, the reverse is not as evident due to the disparate degree of freedoms. Praprotnik et al. (2005) developed the adaptive resolution scheme (AdResS) in which MD is coupled with a coarse-grained version of MD which is able to blend between the two descriptions, changing the degrees of freedom on the fly. The inter-molecular force is blended between the two descriptions as where w(x) is a blending function taking values from zero to unity, and α and β are the centre of masses of the two interacting molecules. The superscript ex denotes the explicit treatment (MD) and cg the coarse-grained version. While removing degrees of freedom and effectively increasing the molecular size, the approach is blended from a microscopic to a mesoscopic level at which boundary conditions may be imposed easier than on a pure MD boundary. Ren (2007) investigated the stability of the state and flux coupling scheme as well as a combination of the two. Specifically, the VV, FV, VF and FF coupling schemes were investigated, where we have V = velocity (state), F = flux and the first letter indicates the direction from continuum to atomistic and the second from atomistic to continuum. (11) VV and FV were found to be stable for both small and large sampling intervals. VF performed well for small sampling intervals, while the purely flux-based FF scheme was claimed to be weakly stable. A stability analysis was carried out which confirmed the above findings and showed an amplification factor of unity for the FF scheme which, in combination with statistical noise, caused the weakly stable nature.
When dealing with multiscale methods, there are several ways to link the atomistic level with the continuum. The one used and discussed so far is known as the domain decomposition (DD), where the computational domain is split into sub-domains for which either an atomistic or continuum description is used while the coupling happens at the HSI. An alternative description is the heterogeneous multiscale method (HMM) (Weinan et al. 2003) where the entire domain is covered by a continuum solver and the microscopic part enters the computation locally at nodes where the continuum description is invalid. Asproulis et al. (2012) developed the point-wise coupling (PWC) based on the HMM and showed that this approach delivered good results for the velocity profiles in a start-up Couette flow. In their investigation, they placed the local atomistic regions at the wall nodes and analysed a range of parameter in the Lennard-Jones potential, as well as different channel heights and wall geometries. As evident, the approach lends itself to create bespoke multiscale flow domains and use the atomistic description only where it is necessary. Borg et al. (2013) extended the idea of PWC to fieldwise coupling (FWC). It operates in a similar fashion to the PWC approach and combines its strength with the domain decomposition. Sub-domains (fields) are placed continuously in the continuum domain but do not need to coincide with the continuum nodes as in the PWC. In a series of 6 steps, the continuum solver projects its solution on the atomistic level for which a new solution is obtained and imposed back onto the continuum. The stresses are corrected afterwards. As with the PWC, the atomistic domain can be tailored to the flow field and freely placed inside the macroscopic domain. Applied to a 1D Poiseuille flow with Newtonian and non-Newtonian fluids, they showed that by using only one microelement (field) near the wall, the resolution was enhanced while further elements needed to be placed inside the fluid domain to reduce the error compared to a full MD solution.
The Couette flow is usually used to test hybrid schemes as it is easy to set up and has an analytical solution. It does not require a complex computational domain, and hence, unstructured meshes, as are widely used in continuum CFD solvers, are not usually used. Borg and Reese (2008) set out to develop a framework that incorporates unstructured meshes and described their implementation in the opensource CFD solver openFOAM (Open Field Operation and Manipulation, http://www.openfoam.com/) towards a general purpose coupling approach. Although the implementation has been described in detail, no actual simulations were presented. Macpherson and Reese (2008) introduced the arbitrary interacting cell algorithm (AICA) which was designed to obtain the particles inside the cut-off radius r cut , i.e. for building neighbour lists, within an unstructured framework. Geometrical constraints as well as parallel implementation issues have been addressed and discussed in detail. Using an unstructured mesh for a domain decomposition-based multiscale approach would require the HSI to handle complex interfaces. Borg et al. (2010) used a state controller to impose macroscopic conditions on the atomistic level. It has an actuator which is controlling the atoms, e.g. changing their velocity, and a sensor which is measuring the atomistic properties which are then used to drive the actuator in an iterative procedure. They successfully applied this technique to handle the exchange of boundary data inside the HSI for a Couette flow with an obstacle attached to the lower wall (inside the atomistic region). This method is, however, not only applicable for (complex) multiscale boundary data exchange but could also be used to impose non-periodic boundary conditions in pure MD simulations. Steijl and Barakos (2012) reviewed and improved the treatment of non-periodic boundaries. Here, the cut-off distance of particles close to the boundary may exceed the boundary itself and hence will exhibit a lower intermolecular force. They stated that neglecting this force could contribute 30-40 % in density fluctuation. To remove this shortcoming, a similar MD simulation is usually conducted for which the normal force of the particle component can be evaluated which is sampled and then imposed on particles in the simulation of interest. Furthermore, the authors also treated the tangential force component which has its equivalent in the continuum stress. They applied their methodology to a channel flow with the continuum domain in the centre and MD domains at the walls and obtained good agreement with the analytical solution for the Poiseuille flow. Holland et al. (2014) used MD to simulate a small portion of a channel and obtained the boundary conditions at the wall from it. They extended the channel to a high aspect ratio and imposed the boundary conditions from their MD data. The results showed that without the MD data, the continuum solver did not produce the correct physical behaviour, while with the MD data, a cheap pre-simulation technique was found to reproduce the expected results for geometrical similar, high aspect ratio channels. Even when placing an obstacle at the channel centre, the profile for velocity, density and pressure were accurately captured compared to a full MD simulation. Alexiadis et al. (2013) developed a novel Laplacianbased HMM in which the shear stresses are directly obtained from MD computations (an alternative is to use the Irving-Kirkwood equation which is commonly employed). The idea behind their approach is to solve for the Laplacian term, i.e. the momentum equation of the Navier-Stokes equations becomes We have simply solved for the Laplacian term and collected the rest of the equation in Ω(R), where R is the macroscopic position vector. We can do the same on the atomistic level with r being the atomistic position vector as and thus have found a way to approximate the right-hand side of the Navier-Stokes equation, i.e. Ω(R) and, therefore, have solved for the stresses implicitly. Alexiadis et al. (2014) introduced the hybrid taxonomy to classify HMMs based on their order of approximation going from the atomistic to the continuum level. The order of approximation is equal to the order of the gradient that is being approximated. The Laplacian-based scheme is second-order in their taxonomy as a second-order gradient is approximated. HMMs approximating first-order derivatives are therefore first order, while zeroth-order schemes only approximated scalar transport coefficients. With increasing order, the noise may become dominant and they have stated that second-order schemes are currently the limit. To reduce noise, a filter may be applied as done in their study using the Savitzky-Golay filter which could allow higher-order approximation. Strictly speaking, hybrid-hybrid schemes are possible where, for example, velocity is obtained with the second-order Laplacian scheme, while the temperature could be supplied by a full continuum simulation without any coupling. This hybrid-hybrid scheme was then tested in a gravity driven microchannel, and good results were achieved after fine tuning the MD cell dimensions. It was found, similar to Borg et al. (2013), that the MD cells configuration played a crucial part in the accuracy of the overall solution and could deteriorate if inadequate dimensions were chosen.
So far we have limited our discussion to monoatomic flows, but extension to multispecies flow may be simply achieved by accounting for each atomic mass as done by Kim et al. (2012). Instead of treating each atom equally, summing over its mass and force due to the Lennard-Jones potential allows for several species to be considered. They validated their approach for an incompressible flow with argon in a Couette flow and then used the same set-up but different weight ratios for the two types of atoms which otherwise satisfied the properties of argon. Incompressible flows have not been widely used for multiscale simulations, mainly due to high statistical scatter at low-speed flows. Ko et al. (2014) proposed an array of solutions to overcome this shortcoming. The first approach mentioned simulated the same geometry with different initial conditions for a number of times and then averaged the results. Clearly, for a large number of simulations it is expected that the noise will substantially reduce, but this gain comes with a high computationally cost. The second approach made use of a spatial regression which took data from neighbouring cells to improve the results in the cell at which information was exchanged. The third and final approach made use of temporal regression where data from the same cell, at which information was exchanged, were collected over the sampling time. The solution with the lowest statistical noise was obtained using several simulations and averaging the results; however, this approach can be repeated arbitrarily often to arrive at a defined noise threshold. The spatial and temporal regression techniques are more difficult to converge as sampling data are not abundantly available. However, encouraging results were obtained for the Couette flow and could pave the way for incompressible, hybrid schemes. Cosden and Lukes (2013)  In the introduction to this section, we have mentioned the work of Hadjiconstantinou (1999) in which the moving contact-line problem was studied using a multiscale approach. The contact line was created between two immiscible fluid phases which has been achieved, on the MD side, by removing the attractive force in Eq. (2), i.e. by flipping the minus to a positive sign. This was only done for the potential between the two phases, while inter-phase potentials followed the standard form of Eq. (2). The topic of immiscible fluids is little studied in the context of multiscale simulations but has recently gained some interest in the field of polymer blends. Fermeglia and Pricl (2007) investigated the miscibility of such compounds found in industrial applications where MD was used to obtain coefficients which were mapped onto a mesoscale model. A widely found polymer blend in the automotive industry is formed of polymethylmethacrylate and polycarbonate, also referred to as PC-PMMA. They investigated a 70/30 mix of PC-PMMA for different shear rates and found that the miscibility was invariant with respect to shear. The blend remained immiscible and only changed its morphology, where spherical-like structures were found for low shear rates which elongated under higher shear. Only using a compatibiliser may the two phases be mixed. Similar results were found for polycarbonate-acrylonitrile-butadiene-styrene (PC-ABS) blends. Furthermore, the inhomogeneous density field distribution was mapped onto a FEM solver which then produced homogeneously distributed variables on the macroscale. Posocco et al. (2012), who investigated self-assembled monolayers composed of hydrocarbon and perfluorocarbon as a surfactant on gold nanoparticles, mentioned the importance of multiscale modelling in the context of immiscible flow domains. In their multiscale study, MD was utilised to obtain parameters for a mesoscale model (in this case dissipative particle dynamics, see also Sect. 5). These parameters had direct influence on the morphology on the mesoscale. The authors concluded that by influencing the parameters on the microscale, which can be easily done, 3D patterns can be controlled and created on the mesoscale. Furthermore, they stressed that tailored nanoparticles of complex nature may be studied and created for medical applications. Hence, a multiscale approach may offer new insight into the morphology of polymer blends, but its application may be far reaching and enabling the creating of new materials.

Interim conclusion on hybrid and multiscale molecular dynamics methods
Multiscale MD computations have allowed to accurately capture the flow physics at the microscale while extending the computable domain into the mesoscale (Xu and Li 2007), where experimental measurements are feasible (Santiago et al. 1998;Tretheway and Meinhart 2002) and could provide validation data. Experiments may be performed at the microscale for low-speed flows, and hence, incompressible multiscale methods need to be further developed to cope with statistical scatter. Most of the coupling schemes presented were either tested for the Couette flow (steady and unsteady), the oscillatory shear flow or the Poiseuille flow. Efforts towards more complex geometries have been presented, but simulations making use of complex HSIs and unstructured meshes are not known to us. This may be circumvented by using PWC or FWC as the underlying mesh can be structured or unstructured; however, the domain decomposition approach would benefit from a more general description for flows where the domain can be split into different domains, as, for example, in the case of a wallbounded flow. The imposition of boundary data between atomistic and continuum domains remains challenging. We will discuss in Sects. 5 and 6 how the inclusion of a description between MD and the Navier-Stokes equations can remove some of the noise while simplifying the overall coupling procedure.

Direct simulation Monte Carlo method
The direct simulation Monte Carlo (DSMC) method is consistent with the Boltzmann equation but not derived directly from it (Palharini 2014;Shen 2005). It is a method developed for rarefied gas conditions, i.e. where the mean molecular diameter is much smaller than the mean molecular path. In DSMC, real atoms are grouped together and represented by a single particle which is then used for the simulation. The procedure of a DSMC simulation is as follows: 1. Populate mesh with particles, initialise simulation 2. Advect particles 3. Index particles 4. Perform particle collisions 5. Sample flow properties 6. If t = t max proceed, else go to step 2 7. Output solution In step 1, particles are randomly seeded in the mesh obeying computational constraints. Each cell should have about 20-30 particles for low statistical noise (Bird 1994;Palharini 2014), and the time step t needs to be small in comparison with the local mean collision time (Shen 2005). In step 2, particles are advected according to the temporal scheme used and could be as simple as r n+1 p = r n p + v p t with r p being the particle coordinates at time n or n + 1 and v p its corresponding velocity.
Step 3 determines the cell in which each particle currently resides. The collisions between particles are carried out in step 4, and the macroscopic flow properties are obtained in step 5. This process is repeated until the maximum time step is achieved or, if a steady solution is sought, after a time asymptotic solution is achieved.
While the particles are advected deterministically, the collisions occur statistically. The most common collision technique used nowadays is the no time counter (NTC) method introduced by Bird (1989) which removed unrealistic collisions rates found in highly non-equilibrium flows in its predecessor, the time counter (TC) method, previously introduced by the same author. First, the total number of collisions per cell is calculated as where N is the number of simulated particles inside the cell, N is the average number of particles during the sampling time, F N is the number of real atoms represented by each simulated particle, σ T is the total collision cross section and obtained from an appropriate molecular model, c r is the relative speed and V c is the cell volume. The simplest of molecular models to obtain σ T is the Hard Sphere (HS) model where we have σ T = πd 2 12 and d 12 = (d 1 + d 2 )/2 , d 1 and d 2 being the diameters of the colliding particles. Its simplicity comes at the cost of lacking physical soundness, and various models have been introduced in the past to overcome these shortcomings, a detailed description of which can be found in Shen (2005) and Bird (1994).
For each cell, we calculate N coll and loop over it. At each instance, we select a particle pair and test whether collision occurs via Equation (15) gives the collision probability. It is compared against a random number R n ∈ [0; 1] and if P coll > R n is true, the selected particles are accepted for collision. We should note that (σ T c r ) max is difficult to estimate a priori and so is subject to be updated if (σ T c r ) > (σ T c r ) max . This does not tamper with the validity of the solution as it appears in the nominator of Eq. (14) and denominator of Eq. (15).
The collision occurs elastically, and energy is conserved. The velocities of the particles prior to collision are where we have the centre of mass velocity as and the relative velocity is The velocities after the collision for each particle can be expressed as Since an elastic collision is assumed and energy is conserved, the magnitude of relative and centre of mass velocity has to be conserved, i.e. we have |c r | = c * r and |c m | = c * m where post-collision values denote with an asterisk. In the HS model, the scattering occurs isotropically, and therefore, each particle's outbound direction is equally likely. We can therefore randomly select collision directions from an appropriate distribution which yields c * r and c * m and Eqs. (20) and (21) can thus be solved. The macroscopic flow properties are then obtained with the following relations where n represents the number density. Boundary conditions are easily imposed in the DSMC method. The simplest form is periodic boundaries in which particles reenter the computational domain via a periodic interface. If a particle travels past a solid boundary during the advection step, its normal velocity component to the boundary is simply reversed and the position updated accordingly. At open boundaries, particles are removed when they exit the domain. At inflow boundaries, one needs to determine the velocity and number flux of particles. A popular choice to obtain the velocity of the particles is the acceptance-rejection method in which a random velocity is drawn from a suitable distribution function f(v). This velocity will be accepted if f (v)/f max > R n , where f max is the maximum value of f(v) and R n a uniformly distributed random number. The number flux of particles across an open boundary can be calculated as where β is related to the most probable thermal speed as β = (2RT ) −1/2 with R being the specific gas constant and T the temperature and K t is a normalisation constant and given by where S is given by S = βV, V being the normal component of the velocity on the boundary and we have introduced erf() as the error function. For further information on boundary conditions, the reader may wish to consult Bird (1994), Lilley and Macrossan (2003) and Xiaohai (2005).
At this point, we wish to emphasise that despite a computational mesh is used, no information is exchanged with neighbouring cells, and hence, the choice of structured or unstructured meshes does not pose any higher implementation effort. When using unstructured grids, however, the neighbour cells need to be determined in order to change the cell ID once a particle passes from one cell to another. Fast neighbour search algorithms and implementation guidelines can be found in Löhner (2008). A further indepth review of the DSMC method can be found in Oran et al. (1998).

Review on hybrid and multiscale direct simulation Monte Carlo methods
The DSMC method thrives for flow condition in the rarefied gas regime. This is encountered at near space altitudes of 20-100 km (Xu et al. 2009) and hence lends itself for space applications. Modelling spacecrafts re-entering the atmosphere requires a more complex treatment as the full Knudsen number range is encountered, i.e from free molecular flow (Kn > 10) in the outer atmosphere to the full continuum (Kn < 0.001), while descending towards earth. La Torre et al. (2011) showed that at around Kn = 0.1 the error from a Navier-Stokes solution was twice that of a full DSMC and confirmed that the continuum description is invalid in this region. Although the DSMC method could be extended into the low Knudsen number regime, it would be too computational expensive in comparison with a continuum-based Navier-Stokes solver. It is here that a multiscale DSMC solver benefits from the advantages of both worlds and computational efficiency is gained. Boyd et al. (1995) introduced the gradient length local Knudsen number as where φ could be any flow variable of interest. They referred to Bird (1970) who showed that Kn GLL > 0.05 implies a departure from the continuum behaviour. Lofthouse et al. (2007) investigated a Knudsen number range of 0.002 < Kn ≤ 0.25 for a Mach 10 flow around a cylinder using two independent solver: one based on the Navier-Stokes equations and the other on the DSMC method. They calculated the pressure, shear stresses and heat transfer on the cylinder surface and correlated the results with the local Kn GLL to judge whether the region can be treated as a continuum or not. The shear stresses and heat transfer obtained from the Navier-Stokes solution showed considerable differences compared to the DSMC results in the Kn GLL > 0.05 region. They explained this behaviour to be caused by the inaccurate boundary conditions imposed at the cylinder surface. The pressure was not as adversely affected and showed a better agreement. For Kn = 0.002, the total drag and peak heating was within 1 % of the predicted DSMC result, while for Kn = 0.25 the error increased to 26.2 and 32.1 %, respectively. The study showed the need for hybrid solvers where the fast Navier-Stokes solution should be applied in regions of Kn GLL < 0.05 and the DSMC method elsewhere. This study crucially pointed out once more that an inaccurate boundary condition could cause a departure from the expected results despite the fact that the used method may be applicable away from the boundaries. We have shown results of Holland et al. (2014) in Sect. 2.1 that boundary conditions obtained and sampled from an appropriate molecular or particle approach can extend the Navier-Stokes equation further into the smaller length scales. Thus, multiscale method present a cure where the exact boundary conditions are uncertain on the macroscale, but not on the microscale. Wu et al. (2006) developed a multiscale DSMC solver based on unstructured meshes. The breakdown parameter was based on Eq. (27), where φ was chosen to be the density, velocity and temperature so that Kn GLL was obtained as Kn max = max Kn ρ , Kn u , Kn T , and a second breakdown parameter was used as where the indices Tr and R refer to the translational and rotational temperature. The initial solution was obtained via the Navier-Stokes equations, and if one of the two continuum breakdown parameters exceeded a threshold value, of which several values were tested, the domain was flagged for the DSMC solver. They tested their solver against a 2D flow over a wedge at 25 • of inclination and a 3D near-continuum parallel orifice jet that expanded into near-vacuum. Good agreement between the full and multiscale DSMC results were attested and further experimental data for the 3D flow case confirmed this. Expanding on this work, Lian et al. (2011) used Eq. (27) with φ being the pressure. This has the advantage that the high-velocity gradient in the boundary layer, which may still be in thermal near equilibrium, is not excluded from the continuum domain, and hence, a costly DSMC calculation is avoided. They applied this to the same 2D, 25 • inclined wedge flow at Mach 4 and to a Mach 12 square cylinder flow. The L2 norm of the density and temperature showed that the error was decreasing faster and achieved a slightly lower overall error compared to the breakdown parameter used before without diminishing the overall accuracy of the solver. Pantazis and Rusche (2014) coupled a DSMC and Navier-Stokes solver in open-FOAM for complex, three-dimensional and unsteady flows. The authors emphasised that for unsteady flow coupling, numerical constraints become more restrictive and may deteriorate the solution and the efficiency of the multiscale scheme. Numerical noise, for example, may not as easily be reduced if the sampling interval is low and cannot be increased due to the time step requirement. The efficiency may suffer from a poorly designed parallel implementation. Any problems of this sort will build up over time and could cause the hybrid solver to perform worse than a monoscale approach (longer computational times) and diverge due to self-induced numerical instabilities. For geometric flexibility, unstructured meshes were supported. They tested their hybrid solver on the shock tube problem and the unsteady flow through an orifice. Their results closely matched those obtained with a pure DSMC solver. A comparison of the ideal and real speed-up further revealed good strong scaling capabilities. So far in our discussion, we have only considered domain decomposition and heterogeneous multiscale methods and variations of these. These may not always be the best choices. If we consider a channel of height H that is of the order of the mean free path , as is the case in a Knudsen compressor, then either method (DD and HMM) would be difficult to implement. Using a continuum mesh and then refine it using the HMM approach, the mesh spacing inside the channel would be smaller than and, therefore, microelements placed on the nodes would overlap. One possible solution would be to make use of the fieldwise coupling approach introduced in Sect. 2.1, which is what Docherty et al. (2014) have done. They used a hybrid Navier-Stokes DSMC code and considered the 1D micro-Fourier flow. As has been mentioned by Borg et al. (2013), the size of the microelement, which consists of a sampling and relaxation zone, needs to be appropriately defined in order to not adversely affect the results. Hence, a parametric study is needed a priori which has been done by the authors and the size of microelements had been adjusted to match a full DSMC simulation. Their approach worked well for the tested case but is currently limited to 1D flows. Therefore, the reported speed-up was only marginal as a wide scale separation could not be achieved. Another possibility is to join all microelements along transversal lines together, to form one single, cross-sectional microelement. This has been done by Patronis et al. (2013). They refer to this method as the internal multiscale method (IMM). They applied their new technique to a converging-diverging channel driven by an external acceleration, a curved, high aspect ratio channel driven by a pressure difference and the developing flow confined by two eccentric cylinders at different rotational velocities. Their approach agreed well against full DSMC simulations. The computational speedup was reported to be 6, 50 and 300 times faster, respectively. The very high speed-up was due to the steady-state continuum solver for the multiscale method, whereas the full DSMC solver was fully unsteady. Patronis and Lockerby (2014) extended this method to the low-variance deviational simulation Monte Carlo (LVDSMC) as introduced by Homolle and Hadjiconstantinou (2007). It is based on the Boltzmann equation but retains the structure of the DSMC algorithm. The key idea is that the velocity distribution can be decomposed into an equilibrium part which can be solved analytically and a deviational part which is solved for by the simulated particles. Since the equilibrium is obtained in this way, small deviations from equilibrium can be simulated with higher signal-to-noise ratio and the overall noise decreases. For the case of the convergingdiverging channel, good agreement was obtained compared to the DSMC results.
Statistical scatter is omnipresent and will remain a challenging problem due to the inherent stochastic treatment of the collisions. Averaged quantities based on the particles contain a high variance and the restriction imposed by current and foreseeable future generation of high-performance computing facilities means that researcher has to opt for noise reducing method to control the accuracy instead of increasing the particle count. Burt and Boyd (2008) introduced a low diffusion variant of the DSMC which replaces the collision by a set of explicit deterministically steps. Particles are seeded in cells and not allowed to freely travel from one cell to another, which is enforced by specular walls on the cell. This removes the necessity to index particles at every time step but requires to transfer and exchange momentum in-between neighbour cells. Furthermore, the cells are not rigid but move with a cell velocity obtained from the mass averaged particle velocities. Motivated by the successful application to the 1D shock tube, they coupled the LD variant with a conventional DSMC solver and simulated the flow around a cylinder at Mach 6 and a nozzle/plume expansion (Burt and Boyd 2009) as well as the flow around a cylinder at Mach 20 (Burt and Boyd 2010). The use of the LD-DSMC method to increase computational efficiency was justified by the fact that particle trajectories were correctly integrated while moving over cellbased length scales much larger than the mean free path. The results compared well to the reference solution and the speed-up for the Mach 20 cylinder flow was given to be about 2.6 higher compared to a full DSMC simulation. Jun et al. (2013) calculated the flow around a cylinder as well, this time at Mach 10. Compared to a full DSMC simulation, their speed-up was only about 1.2 using the same method as described above. They argued that a considerably large part of the calculation is spent finding the initial solution with the DSMC method which, for small Knudsen numbers, i.e. continuum regions, is rather time-consuming. Using a continuum solver instead to obtain the initial solution, improved their speed-up to a factor of 2.6. Sun and Boyd (2005) introduced the sub-relaxation technique in which the history of a signal is removed during averaging once it starts to deteriorate the averaged result. A relaxation parameter θ can be freely set to influence the averaging procedure but may negatively affect the results and so a parametric study is advisable. Schwartzentruber and Boyd (2006) used the sub-relaxation technique in a hybrid Navier-Stokes DSMC solver and used it for boundary condition imposition at the interface between DSMC and continuum. They tested their approach on a onedimensional shock tube for liquid argon and nitrogen gas. Their technique was able to provide smooth boundary conditions for the hybrid solver which compared well to full DSMC and experimental data. A pure Navier-Stokes solution was also tested which gave inaccurate results with the reciprocal shock thickness overestimated by almost a factor of 2. This is consistent with the observation of Carlson et al. (2004), who computed the flow inside a shock tube for Mach numbers of 1.55, 5 and 10. In their multiscale solver, the Navier-Stokes equations were coupled with the DSMC method based on the IP method (information preserving, Fan and Shen 2001). In this method, each particle is assigned a second, information preserving, velocity. This velocity is solely used to obtain the macroscopic quantities through particle averages but does not influence the movement of the particles and hence is less prone to statistical scatter. The IP velocity follows its own set of rules to update itself, a detailed list of which can be found in Shen (2005).
The results for the Mach 1.55 case yielded overall good agreement amongst the Navier-Stokes, DSMC and multiscale solver. Increasing the Mach number decreased the accuracy of the continuum solution, whereas DSMC and the hybrid version matched each other well. The reciprocal shock thickness confirmed a deviation by a factor of about 2 for the Navier-Stokes equations from the experiment, while DSMC and the multiscale solver matched the experiment well. The speed-up was also mentioned, and despite the multiscale solver being about 1.6 times faster than the pure DSMC version, the Navier-Stokes equations still yielded results 9 times faster than the multiscale solver. Thus, for those cases where the Navier-Stokes equations can be modified to incorporate appropriate boundary conditions, it is still cheaper while retaining a good degree of accuracy. Boyd (2008) investigated the flow past a cylinder-flare and blunt planetary probe using a hybrid DSMC Navier-Stokes solver. Different mesh densities and time steps for the respective solver were used in order to speed up the calculation. The initial solution on which the domain has been separated into rarefied and continuum regime was obtained by the Navier-Stokes equations as they have previously proven to provide the initial solution in a shorter time (Jun et al. 2013). It has been mentioned that the mean free path varied over one order of magnitude for the cylinder-flare case (speed-up 1.4) and over two orders of magnitude for the blunt planetary probe (speed-up 12.5), where the speedup in parenthesis are comparing to full DSMC simulations. This shows that multiscale solvers are becoming more efficient for a wide range of Knudsen number separation and it is under these conditions that they thrive.  Schwartzentruber and Boyd (2015) pointed out that experimental studies are vital and needed for the validation of current DSMC codes. Currently, most multiscale and hybrid DSMC codes are validated against pure DSMC solvers, see Scanlon et al. (2010) for a range of 1D, 2D and 3D test cases. They pointed out that MD could be used as a validation tool but is computational rather expensive for rarefied gases. Another issue that arises when coupling DSMC with a second solver is the noise penalty that arises due to the stochastic treatment of the collision. A range of approaches have been presented herein (LD-DSMC, LVDSMC, IP, sub-relaxation) and have all been able to reduce the noise level. As seen in Ko et al. (2014) (see Sect. 2.1), using spatial or temporal regression reduced the noise with little extra computational effort while being independent of the method and hence could prove to be a viable alternative or even be used in conjunction with one of the abovementioned methods.

Interim conclusion on hybrid and multiscale direct simulation Monte Carlo methods
The reported computational gains were highest for the cases with the highest Knudsen number separations, and hence, applications may be limited to these sorts of flow regimes. However, space applications have shown to exhibit the free molecular to continuum regime and thus rendered themselves as prime candidates for multiscale DSMC schemes.

Lattice Boltzmann method
The lattice Boltzmann method (LBM) is strictly speaking not a direct particle method. It is derived from the lattice gas cellular automata (LGCA) which in turn is based on physical particles. In the LGCA method, particles reside on so-called sites which are connected via links. It is essentially a computational mesh where sites refer to nodes and links to edges. The governing equation simply states that the number of outgoing and residing particles at a site is equivalent to the collision occurring at the same site, i.e. Several collision models have been introduced, most notably the one of Hardy, Pomeau and de Pazzis (HPP model, Hardy et al. 1973Hardy et al. , 1976 and Frisch, Hasslacher and Pomeau (FHP model, Frisch et al. 1986). Due to its binary nature (particles either do or do not exist at sites), the method is free of any numerical round-off errors. However, statistical noise is high and averaging over larger regions or several simulations is necessary to obtain macroscopic quantities. The model further lacks Galilean and rotational invariance. The interest in the method declined once the problems surfaced. In the LBM, the single particles are replaced by a density distribution function which removes some issues encountered in the LGCA. That, however, also means that the collision operator cannot be treated explicitly and has to be approximated numerically. The easiest and still one of the most widely used model is the one of Bhatnagar, Gross and Krook (BGK model, Bhatnagar et al. 1954). It is sometimes referred to in short as the LBGK (Lattice-BGK) model and given by where τ is the collision time, f k is the density distribution at link k and f eq k is the corresponding equilibrium distribution. The collision time is given by where ν is the lattice viscosity and ω is the collision frequency. The full governing equation with the BGK operator is where Eq. (32) can be seen as a simple advection equation with a source term. If we ignore the unsteady term, we see the similarities to the LGCA method, i.e. Eq. (29), where we have streaming on the LHS and collision on the RHS. To solve the above equation, we need to define c and f eq k . The lattice speed c is determined by the speed model. The speed model determines how the distribution function is streamed to its neighbours. Typically on a 2D structured mesh, there are 9 sites (nodes) to consider, located at r neighbour = r(r x ± �x, r y ± �y) which includes the site at x = y = 0. Therefore, the speed model would be termed D2Q9 (2D, 9 sites) and each site has a corresponding speed c k and weight w k . The speeds stream from one site to another and, therefore, take integer values of −1, 0, 1 . They are obtained as where The corresponding weights are The weights are to some degree arbitrarily set but need to fulfil certain constraints, one of which is that w k = 1 . More information on the weight constraints and how to develop higher-order speed models can be found in Succi (2001).
The density distribution function at equilibrium is given by where we have the macroscopic density ρ(r, t) and velocity u, while the lattice speed of sound is obtained from c s = 1/ √ 3. The macroscopic quantities are then obtained from Eq. (32) is now solved in the following way: Calculate Ω from Eq. (30) 4. Stream the updated density distribution function to their neighbour sites 5. Obtain macroscopic quantities for ρ and u via Eq. (36) and Eq. (37) 6. Repeat from step 2 until t = t max 7. Output results Due to the nature of the LBM in which macroscopic properties are determined via the density distribution function, we need to find boundary conditions for each contribution of the density distribution function f i on the boundaries. For the D2Q9 model that we have discussed above, there will be 3 components of f i outside the computational domain, 3 components on the boundary itself, i.e in-plane and 3 components pointing into the computational domain for which we need to formulate boundary conditions as they cannot stream from outside the boundary into the domain. There are several possibilities to find appropriate formulations for the inward pointing components of f i and we will briefly describe the method of Zou and He (1997) which is one popular choice. Let us denote the inward facing components of f i as f − , f c and f + where f c is the component normal to the boundary and f − and f + the component to its left and right, respectively. We first need to evaluate the density, which we can write for a Cartesian domain in compact form as where ±u ⊥ is the perpendicular velocity component on the boundary and the sign is positive if the normal direction (pointing outside the domain) of the boundary is positively aligned with the coordinate system. We differentiate further between the components of f i which are inplane, i.  (40) and (41) will be positive if the normal (pointing outside the domain) on the boundary is negative, and negative if the normal is positive. 1/2 ± ρu � in Eq. (40) and Eq. (41) on the other hand will be positive for f + and f − , if their direction with respect to the boundary is in positive coordinate direction and negative otherwise. The velocity component of u ⊥ is either prescribed via Dirichlet-or Neumann-type boundary conditions and u can be set for moving boundary problems. A comprehensive overview of various boundary conditions can be found in Latt et al. (2008) and Chen et al. (1996). Classical and recent applications of the LBM can be found in Chen and Doolen (1998) and Aidun and Clausen (2010).

Review on hybrid and multiscale lattice Boltzmann methods
The LBM has found wide acceptance and applicability among various research discipline, mainly due to its ease of implementation, handling of complex geometries and low computational cost along with its local behaviour which makes it a prime candidate for parallelisation. However, the LBM has also been found to work less well for flows involving heat transfer, compressibility and high Reynolds numbers. The main reason for the low-to-moderate Reynolds numbers that can be achieved with the classical LBGK model is due to the instability that arises with low lattice viscosities. One possible solution is to decrease the lattice spacing at the cost of increased computational cost. One could also use a multiple relaxation time (MRT) scheme (D'Humières et al. 2002) instead of the BGK operator for the collision. The main idea of the MRT scheme is that the collisions occur in moment space rather than velocity space and each moment is relaxed towards equilibrium at its own rate, while in the BGK model, all moments are relaxed at the same rate 1/τ. That allows for lower viscosities and hence higher Reynolds numbers. Karlin et al. (1998) introduced the notion of the entropic lattice Boltzmann method (ELBM), where the H-theorem is incorporated into finding the equilibrium distribution function f eq k in an attempt to adhere to the second law of thermodynamic, i.e. positive entropy production. In its limit, the ELBM will result in the classical LBM and subsequently allows lower viscosities as well. Compressible effects are difficult to model since the small number of discrete velocities in the speed model allows only small temperature variations. Successful applications have been presented, for example, by Guangwu et al. (1999), who simulated Sod's and Lax's shock tube problem. Despite following the trend of the reference solution, spurious oscillation, as well as numerical dispersion, was present. Compared to other published data obtained with the Navier-Stokes equations using different Riemann solvers and numerical schemes, its L1 norm in density, pressure, velocity and energy compared less favourable. Joshi et al. (2010) constructed a hybrid LBM in conjunction with a finite volumebased Euler solver where the primitive variables at the cell centres were obtained from the Euler equations, while the inter-cell fluxes were approximated by the LBM through a modified equilibrium distribution function to account for thermal effects. They applied their scheme to the 1D and 2D Sod shock tube problem as done by Guangwu et al. (1999). Despite the absence of tabulated data in the work of Joshi et al. (2010), their scheme qualitatively exhibited the same behaviour as the pure LBM, see Guangwu et al. (1999). Compared to two solutions obtained from the pure Euler equation making use of the Godunov and Flux Vector Splitting (FVS) scheme, their solution compared well to the Euler results with spurious oscillation and numerical dispersion still present. To incorporate heat transfer and to compute the temperature field, as stated previously, the equilibrium distribution function needs to be modified and higher-order speed models should be used, i.e those that take their immediate neighbours and their neighbours into account. Various models have been introduced to account for the energy and thus for the temperature (Guangwu et al. 1999;Kataoka and Tsutahara 2004a, b;Shi et al. 2001) and are reviewed in Kun (2008). Mezrhab et al. (2004) discussed the thermal models available in the Lattice Boltzmann framework and concluded that these are still premature and need further development. They opted for a hybrid approach where the LBM was solved alongside an advection-diffusion heat equation to obtain the temperature field. They applied their coupled solver to a range of two-and three-dimensional cases and obtained results that compared well to reference data. The same approach was chosen by Wu et al. (2012) to simulated forced convection around a cylinder at various Reynolds numbers. The cylinder was resolved by the immersed boundary method (IBM) and its influence introduced into the advection-diffusion heat equation by a source term. Since the boundary itself did not coincide with the computational nodes as in body fitted meshes, evaluating gradients at the boundary was not straight forward. To calculate the Nusselt number, which contains the normal temperature gradient at the cylinder surface, it was recast using Fourier's law −Q/k = ∂T /∂n| Ω to eliminate the gradient and then finding an appropriate relation for the heat flux Q. Their results were compared against two other numerical simulations, and the discrepancies were at most 5 %. Chen et al. (2013) extended the advection-diffusion equation approach to handle species transport. Their test cases were species convection and diffusion with bulk reactions, species diffusion in a channel with surface reactions and natural convection inside a square cavity. The first two test cases were compared against analytic solutions and good agreement was attested while the same could be said about the third test case compared to the solution of a commercial Navier-Stokes solver. With their methodology tested and validated, they set out to simulate a more complex example of a wall-coated microreactor to simulate heat transfer, mass transport and chemical reactions. The velocity profiles were compared to a full-scale lattice Boltzmann simulation for which good agreement was obtained. To further accelerate the procedure, the grid spacing for the continuum domain was increased while keeping the spacing of the lattice Boltzmann domain the same. It has been found that at a grid spacing ratio of 10:1, results were still good while maximising the computational turnaround time.
Yet another study has been done by Filippova and Hähnel (2000), using the same approach as chosen by the authors in the previous studies, and extended their approach to variable density flows in order to simulate chemical reacting flows. The LBM was further extended with adaptive mesh refinement capabilities and boundary fitting to account for the misalignment of the boundaries. They have applied their method to the flow around a porous burner, where the burner itself (cylindrical shaped) has been resolved with a refined grid. The ratio tested was 3 and 6 times finer than the far field grid. Since the time step is directly proportional to the lattice spacing, i.e. �x k /�t = c k = const., a refined grid means that the time step decreases as well. For reactive flows, this becomes important as the timescales for reactions are rather small. Studies were conducted for a range of Reynolds and Damköhler numbers, matching the results of full Navier-Stokes solutions. Li et al. (2014) simulated the melting and solidification process in a 2D rectangular domain where the west wall was above the melting temperature, the east wall below the melting temperature and the north and south wall adiabatic, i.e. ∂T /∂n| Ω = 0. The velocities were obtained by the LBM, while the temperature was calculated from the energy equation which was solved by a finite volumebased SIMPLE (semi-implicit method for pressure linked equations, Patankar and Spalding 1972) scheme, see also Patankar (1980), Versteeg and Malalasekera (2007) for a detailed description of the SIMPLE scheme and its derivatives. The variables for the finite volume solver were stored at the cell centres while the LBM stored its variables at the vertices. In order to provide face-centred velocity components for the SIMPLE algorithm, the velocity needed to be interpolated to its required location. The interface was obtained via the interfacial tracking method. The method identifies the cells in which the interface currently resides and assumes that the cell temperature is at the melting condition, T m . This is, however, only the case if the interface is exactly at the cell centre and, therefore, a correction was introduced to allow the interface to be off-centre. Simulations were then carried out for various Stefan numbers, Ste = c p �T /L where L is the latent heat, and compared to experimental data. The bulk of the interface was agreeing well with the experiment. At the wall, discrepancies were observed and attested to the difficulties to maintain the adiabatic condition in experiments, as well as the assumption made in the numerical simulation, that the volume change during the phase transition can be neglected. Brent et al. (1988) developed an enthalpy-porosity technique for melting and solidification processes in which the Navier-Stokes equations remain valid throughout the flow for the liquid, transition and solid phase. The energy equation contains an additional source term in which the latent heat evolution is modelled for T > T m ; otherwise, the source term is zero. The momentum equation is equipped with a source term that mimics the Kozeny-Carman equation. It is zero for the liquid phase and increases in the transition region. For a fully solidified solution, its value is great enough to force all the other terms to vanish, enforcing a zero velocity field. The advantage of this approach is that the interface comes out as a product of the solution and does not need to be explicitly tracked. Chatterjee and Chakraborty (2006) reported results of a dendrite growth in an undercooled melt making use of the enthalpy-porosity technique. They used a modified thermal LBM from which hydrodynamic variables were computed and coupled it with an enthalpy density distribution function from which the thermodynamic variables were derived. Their results were reported to compare well against reference data although only qualitative figures were presented. However, their approach coupled the fast, modified (hybrid) LBM with the elegance of the enthalpy-porosity technique and hence lent itself to fast computational times for complex interfaces. Ladd (1994a, b) developed a general framework and solver (SUSP3D) to simulate suspended particles. Forces are imposed on the particles via fluctuating stresses rather than applying random forces directly to the particles. They are advected using classical Newtonian dynamics including angular momentum of the particles. The solver has been developed with an emphasis on good scalability for a large number of particles and hence has been subsequently used by various other researchers. Beetstra et al. (2006) made use of Ladd's solver and investigated the drag behaviour of various particle formations. They observed that particles in a cluster formation exhibited a lower average drag force per particle than a single particle. This was explained by the increased surface-to-volume ratio as each individual particle, on average, has less surface area to resist the oncoming flow. For a cluster of particles whose formation represented that of a sphere, at the closest packing possible (r/d = 1 ), the average drag force per particle was only about 20 % of that of a single particle and only reached the same level of drag per particle for an inter-particle distance of r/d = 7 . Hence, the effect of clustering felt by each particle sustained several inter-particle distances away. Similar results were observed for other shaped clusters. They further investigated how the drag behaviour changed by increasing the Reynolds number and compared that against a single-particle drag law presented in Clift et al. (1978). They showed that for increasing Reynolds numbers, larger deviations from the drag law were observed. Following their work, Beetstra et al. (2007) investigated mono-and bidispersed particle systems, the latter with diameter ratios of 1:1.5-1:4. They derived a drag force relation for both cases and showed that failing to include the bidisperse nature in the drag law yielded significant differences. Their results were able to not only match the trend but also to accurately predict the particle drag for various Reynolds numbers and inter-particle distances. Furthermore, they showed that their model and simulation data could differ by as much as 3 times from current drag laws which, despite their shortcomings, are still widely used. Shah et al. (2013) investigated the drag behaviour of particle clusters using SUSP3D. Their cluster was surrounded by randomly seeded particles and, therefore, two voidages, that of the cluster and that of the surrounding cell, was defined. They explored an extensive range of various parameters, confirming that compared to the Ergun and Wen-Yu drag law, particle clusters exhibit a different drag behaviour.
Furthermore for a constant cluster voidage of 0.7, they found that a minimum for the drag existed at an overall (cell) voidage of 0.96. Increasing the overall voidage further resulted in the drag to approach that of randomly seeded particles while decreasing the voidage resulted in a similar behaviour at a voidage of 0.85, which again demonstrated the importance of cluster formation in particle driven flows, such as the fluidised bed reactor. Usta et al. (2006) studied the migration of polymers in microchannels using SUSP3D, where polymers were modelled via particles connected by stiff linear Fraenkel springs. Their work was motivated by the discrepancies in the literature where polymers were reported to either migrate towards the channel centreline (including hydrodynamic interactions) or wall (without hydrodynamic interactions). From a thermodynamic perspective, no migration preference in a uniform shear flow is expected. Kinetic theory, on the other hand, predicts migration towards the centreline for uniform shear and Poiseuille flows. Their findings showed that for uniform shear flow the polymer migrated towards the centre of the channel as long as the channel height was large enough. Decreasing the size increased the influence of the walls, and for small enough channels, the polymer could migrate towards the walls as observed by previous studies. It is worthwhile to point out that this study presented an opportunity to study micro-and nanoscale phenomena using a mesoscale approach by simply incorporating a suitable polymer model into the solver making it a hybrid method. Thus, in comparison with a full molecular or multiscale approach, even more computational savings are to be expected while extending the capability of the solver to handle polymeric effects. Ollila et al. (2011) mentioned that for thermal fluctuations to be present, the density needs to fluctuate as well. Using a squared speed of sound (as, for example, in Ladd 1994a) will zero out the viscosity which is relevant for fluctuating density flows. They incorporated their fluctuations via a forcing (source) term on the RHS of Eq. (32) and showed that the variance of the fluctuation is related to the fluctuation-dissipation in the stress tensor. They applied the new thermal fluctuating LBM (TFLBM) in a multiscale approach in Ollila et al. (2013c). Here they showed results, among other particle drag related studies, for MD particle in a cubic box where the solvent was modelled by their TFLBM. Tracking the mean square displacement over time of each particle, they were able to calculate the macroscopic diffusion coefficient without Langevin noise added to the system. Although the diffusion coefficient was observed to be smaller than anticipated, it has been attributed to the finite size of the simulation domain. Ollila et al. (2013b) further explored the effect of confined polymers using either the TFLBM or Langevin Dynamics (LD) in 2D and 3D. First, they studied the radius of gyration (a measure of how far the polymer is stretching) and the static structure factor (a measure of the equilibrium size of the polymer) in both planar and perpendicular direction with respect to the confinement and showed their different behaviours for various confinement levels. They showed that the TFLBM was able to capture the decrease in planar diffusion coefficient for increasing confinement levels, while LD was not able to accurately describe the polymers migration due to the absence of hydrodynamic interactions. In yet another study, Ollila et al. (2013a) investigated a microfluidic T-junction with one inlet at the south and two outlets at the east and west side. They showed that the separating streamline, i.e. the streamline which separates the streamlines going to the outlet to the east or west, can be controlled by adjusting the outlet pressure at both boundaries, which thus influenced the influx-to-outflux ratio (mass per unit time). Being able to control the separating streamline, they placed two particles in the channel downstream of the inlet with different geometrical parameters (initial separation of the two particles, distance of their centre of mass to the separating streamline, radius and channel geometry) and studied which outlet they would choose to exit. This has been validated against experimental data. In fact, they showed their findings in a phase diagram which made it possible to predict which outlet the particle preferred, based on their initial state. Their findings can thus be used to construct a microfluidic NOT gate (one particle) or NAND gate (two particles).  used the TFLBM to couple it via local forces, conserving momentum and energy, to MD particles. Here, the mass at the lattice nodes was used and coupled with the MD particles directly. MD particles were used to construct spheres so fluid could be captured inside. By changing the mass of the MD particles, either the mass of the fluid inside the spheres was much larger or smaller than the overall mass of the sphere (combined MD particles). Studying the velocity autocorrelation function (VACF), they showed that the long-time asymptotic behaviour was independent of the particles mass and agreed with theoretical findings. The particle mass only influenced the short-term behaviour of the VACF.  implemented the TFLBM into the open-source solver LAMMPS and validated their solver against drag on a single spherical particle, particle motion near a solid wall, hydrodynamic interactions between four particles without collision, two particles with collisions (hard sphere interactions) and a confined colloid undergoing a Couette flow. They supplemented their validated cases with a scalability analysis based on two high-performance computing clusters. Cui et al. (2012) simulated the cavity formation in a particle covered domain, resembling a pipe leakage that initially was covered by soil particles. They used the discrete element method (DEM) for the particle motion, with and without surface energy for particle cohesion, and the LBM for the fluid surrounding the particles. The particle velocity and force were obtained from the fluid motion which in turn was influenced by the particle motion. They used a rectangular computational domain with an orifice at the bottom wall, located in the middle from which fluid was injected. The domain above was covered by particles. For relatively low exit velocities, the pore pressure was found to rise and stabilise over time. Increasing the exit velocity showed the formation of a cavity which was accompanied by a sudden drop in pore pressure. For even higher velocities, a blow out failure was observed where the fluid was flowing in a straight path to the particle surface, rupturing the layers of particles above. For moderate velocities, the cavity reached a stable size and stopped growing. The inter-particle cohesion was found to have no influence on the rate of orifice pressure build up. A greater pressure was, however, required to initiate the cavity formation as the inter-particle attraction forces needed to be overcome.
So far in our discussion, the particles suspended in the fluid have been modelled as rigid bodies. Dupin et al. (2007) coupled the LBM with deformable particle (DP) dynamics. The advantages of DP are that their membrane deforms as a direct response to the fluid motion. This approach is favourable in cases where the particle properties need to be resolved explicitly and influence the flow, however, that comes at an increased computational cost. Coarse meshes for a single DP have been reported to have approximately 500 nodes. Three types of particles were considered; capsules, spherical in shape, having a finite thickness and unbreakable membrane; vesicles, spherical in shape, having an infinite thickness and unbreakable membrane and red blood cells (RBCs), biconcave in shape. The motivation for this study was driven by the fact that studies up to date, using the computationally more expensive membrane potential function, only allowed small domains to be simulated. They validated their approach against several test cases and used their solver to simulate 200 RBC at low Reynolds numbers in a microchannel. RBC were initially placed 1 micron apart in a regular formation and then allowed to reach an equilibrium state. They found that the pressure field varies significantly in time and space. This fact is usually ignored in biological models and cells subjected to this pressure filed showed great deformation which could have a significant impact on biological processes such as gene expression, a process in which gene information are used to synthesis a functional gene such as a protein. Despite the focus on high scalability, with current limitation imposed by HPC facilities, simulations over the order of 100 cell diameters were not possible. However, it was argued that for such length scales the effect of the particles on the fluid can be modelled and used as a mesoscopic model to formulate boundary conditions for macroscopic simulations, as seen, for example, in Holland et al. (2014), Sect. 2.1.
Driven by the same motivation, to construct a highly efficient solver for general fluid dynamics applications, Zou et al. (2014) combined the efficient velocity, pressure and density computation of the LBM with the Navier-Stokes equations for the constitutive relations. They validated their solver against a Poiseuille flow, Taylor-Green vortex and sudden contraction with a ratio of 4:1 and obtained good agreement. Their focus then shifted to the computational efficiency of the hybrid solver. Comparisons for all three test cases showed that the solver coupling lattice Boltzmann and Navier-Stokes needed 20-25 % of the time compared to a full Navier-Stokes solution. The reduction was attributed to the removal of the pressure correction loop inside the PISO (pressure implicit with splitting of operator, Issa et al. 1986) algorithm that has been employed in the full Navier-Stokes solution. This is in accordance with the findings of Premnath et al. (2005) who commented that the pressure calculation in the Navier-Stokes equations usually takes about 80 % of the whole computational time. Salimi et al. (2015) presented a hybrid lattice Boltzmann/Navier-Stokes solver incorporating thermal effects directly into the LBM. Here the thermal LBM was coupled with the Navier-Stokes equations and the temperature, and the other hydrodynamic variables were exchanged at an interface. One of the goals was to elucidate the efficiency of the hybrid LBM solver, for which they first validated it against reference data for different Reynolds numbers, porosity levels, solid to fluid thermal conductivity and blockage ratios. They then simulated a fixed size domain with a pure LBM, pure Navier-Stokes and hybrid LBM for one lift cycle. When the time step in the hybrid and pure Navier-Stokes solver was set to that of the pure LBM, they needed 6.14 and 8.48 more computational time compared to the pure LBM solver. This is to be expected as the LBM itself performs less floating point operations per time step. At a ratio of approximately �t NS /�t LBM = 50, the hybrid and continuum Navier-Stokes solver performed as fast as the pure LBM solver, taking only 0.92 and 0.84 of its time, respectively. Pushing the time step ratio further to �t NS /�t LBM = 200, the computational time further decreased and took 0.55 for the hybrid and 0.33 for Navier-Stokes solver compared to the full LBM solver. This is not surprising as both methods had to perform fewer time steps to advance to the next time level. It is interesting to note that the pure continuum solver performed better for larger time steps than the hybrid solver. This indicates that for true macroscopic simulations the Navier-Stokes equations are still computationally efficient. In contrast, the hybrid solver was able to incorporate multiphysics phenomena on smaller scales at only slightly higher computational costs which may not be easily achieved with a conventional Navier-Stokes approach.
Thus far, we have seen applications in which the LBM was successfully coupled with particle and continuum domains. Due to its mesoscopic nature, however, it is a prime candidate for coupling with a microscale approach. Dupuis et al. (2007) coupled a LBM solver with MD and applied it for the flow past and through a carbon nano tube (CNT). The domain decomposition technique was used in conjunction with the Schwartz alternating method. State coupling was employed and the influence of imposing only velocities from MD to LBM and velocity gradients from MD to LBM investigated. Compared to a full MD solution, the error decreased for about 10 Schwartz iterations after which no considerable gain in accuracy could be achieved. The velocity gradient imposition was 3 times closer to the full MD solution than the solution in which only the velocities were exchanged. Using this approach, for the flow past and through a CNT, the average error (based on velocity differences) was 1.3 and 2.1 %, respectively. The speed-up for both cases was reported to be 7 and 2.4 times faster using the multiscale approach. This is due to the reduced computational demand away from the CNT where inter-molecular motions are adequately resolved by the LBM. Lobaskin and Dünweg (2004) used MD to simulate a charged colloidal system. Here the MD particles used the Lennard-Jones potential for the inter-particle forces and additionally finitely extendable nonlinear elastic (FENE) springs. Two test cases were presented, the first in which the colloid was given an initial horizontal velocity and the second with an initial angular velocity. For the first case, the colloid velocity autocorrelation function exhibited the expected longterm asymptotic behaviour in which momentum is transported away by diffusion. The colloid with initial rotation showed that its decay of angular velocity agreed well with Debye's law and again showed the expected long-term asymptotic behaviour. Results were then presented for a charged colloid during a self-diffusion process. Counterions were placed around the colloid to give a neutral overall charge. The velocity autocorrelation function of the centre of mass velocity showed good agreement with reference data for both neutral and charged colloid. Farahpour et al. (2013) investigate the translocation process of a singlestranded DNA (ssDNA) in which the ssDNA was passing through a pore, driven by an electrical field. The hydrodynamic interactions were captured with the LBM and the ssDNA was represented by a monomer through a MD simulation. The Laplace equation has been solved for the electric potential and counterions were added to overall balance the charge of the simulated system. By studying the velocity gradient tensor by its eigenvalues and eigenvectors, the principle axes of extension and compression could be obtained. It was found that the monomer stretched before approaching the pore, aligning itself with the electric field lines. It then compressed once passed through the pore.
Further investigations revealed that replacing the counterions by a charged monomer alone produced incorrect approaching probabilities of the monomer to the pore. If the hydrodynamic interactions were further neglected, then field gradient effects were overestimated. Thus, the inclusion of both counterions and hydrodynamic interactions via a multiscale approach was deemed necessary to accurately predict the correct translocation behaviour of the ssDNA. Mackay et al. (2014) investigated the effect of shear flow on colloidal particles in a microchannel. Particles were initially placed in ordered layers into the channel, and each layer had the same volume fraction. An order parameter was defined as ψ 6 (r jk ) = 1/n n k=1 exp[i6θ(r jk )]. Here, θ(r jk ) is the angle between an arbitrary fixed vector and the vector from particle j to k. For �|ψ 6 |� = 1.0, perfect order is achieved in the particle layer. A threshold was defined for �|ψ 6 |� = 0.6 above which an ordered state was still encountered. Systems with a lower value of �|ψ 6 |� were consequently classified as disordered states. Their results showed that �|ψ 6 |� dropped rapidly at the start, reaching a steady state shortly after. The linear shear imposed by the opposite moving walls allowed particles close to the wall to move faster on average, allowing particles in their adjacent layers to jump into the layer due to the inter-particle forces. This behaviour caused a change of volume fraction in each layer, and it was identified as a main cause for a drop in the order parameter. For a greater amount of particle exchanges, the order parameter eventually dropped below the critical value of 0.6 and disorder started to form. Particles close to the wall generally stayed in an ordered formation even though disorder was possible at the channel centre. Their findings have been summarised in a phase diagram showing the disordered, ordered and the coexistence layer for different shear rates and volume fractions.

Interim conclusion on hybrid and multiscale lattice Boltzmann methods
In our discussion, we have outlined that the thermal LBM is still immature and a great amount of hybrid solvers have been developed to use continuum mechanics to solve the temperature field on the mesoscopic scale. This is a valid approach which, however, is more concerned with overcoming the inherent shortcomings of the LBM rather than increasing the computational efficiency. Therefore, most studies discussed herein were concerned with hybrid solvers, while true multiscale approaches with bespoke flow domains in a domain decomposition fashion were less common. In a similar way, high Reynolds number flows are challenging due to the inherent numerical instabilities of the LBM. We have mentioned the MRT collision operator and the entropic LBM which currently extend the range of Reynolds number but do not remove the instabilities.
For applications, however, in which solver efficiency is of importance and the temperature does not need to be computed, LBM enjoys a reputation of being easy to implement for complex geometries and handling multiphysics effects, while it outperforms comparable methods at the same resolution in terms of computational time.

Dissipative particle dynamics method
Dissipative particle dynamics (DPD) is a particle-based method and has high resemblance to MD which can be regarded as a coarse-grained version of it. By grouping atoms together into representative DPD particles, length and timescales are extended into the mesoscopic regime. DPD has evolved as an application driven need to simulate molecular systems for which the current limitations in space and time for MD are too restrictive. We can use DPD if we are only interested in the mesoscopic manifestation of the underlying molecular details. The governing equation is similar to that of MD as where we have assumed DPD particles of same mass m i and unity in magnitude. We see the difference to MD in the splitting of the force into internal and external forces. External forces are those imposed to particles by the system, such as gravitational forces. Internal forces arise due to the inter-particle interactions. We can further decompose the internal forces into where the superscripts C, D and R stand for conservative, dissipative and random, respectively. The conservative force models the particle interactions and is given by where a ij is the maximum repulsion force and the weight function w C (r) is defined as Here, r ij is the distance between two particle centre of masses and r ij the vector between particle i and j of unit length.
The dissipative force mimics the effect of viscosity at the atomistic level and is given by where v ij is the particle velocity, γ a coefficient and again w D (r ij ) a weight function. The random force is given by where σ is a coefficient, ξ ij is a random number satisfying a Gaussian distribution with a variance of unity, and w R (r ij ) is a weight function. The coefficient γ and σ of Eqs. (46) and (47)  Much of what has been said about the boundary treatment for MD is applicable to DPD as well. Periodic boundary conditions are easily imposed, while open or solid boundaries require an appropriate treatment of the conservative, dissipative and random force. We refer the reader to Revenga et al. (1999), Mehboudi and Saidi (2014) and references therein for a detailed treatment of boundary conditions. We have outlined the similarities between DPD and MD and the coarse-grained nature of it. In addition to the conservative force, which essentially replaces the inter-particle potential as seen in MD, we have a dissipative and random force term. The random force adds heat to the system, while the dissipative force is reducing the particle velocities. Their effect together causes the temperature to be approximately constant accompanied by small fluctuations. Hence, DPD is an implicitly thermostated coarsegrained MD version.
An in-depth review of the DPD method and recent developments has been given by Liu et al. (2014).

Review on hybrid and multiscale dissipative particle dynamics methods
Although DPD and the LBM are both mesoscopic in nature, they have been used for different multiscale applications due to their inherent different descriptions. While the LBM is Eulerian based, DPD is a true Lagrangian method and due to its MD background, but upscaled mesoscopic regime, it has primarily found applications in medical and biological flows. Symeonidis et al. (2005) developed a new time integration technique and reported that DPD simulations can be 10 6 times faster than conventional (1 − r ij /r c ) s r ij < r c 0 r ij ≥ r c MD simulations. In its limit, it adheres to both microscopic and macroscopic properties, making it a prime candidate for coupling in both directions. Natural occurring molecules such as wax, fat and vitamins are referred to as lipids, which consist of a head and a tail. In water, lipids form bilayered membrane where the hydrophilic lipid heads form the surface of the membrane with their hydrophilic tails pointing inwards. Biological membranes are made up of different lipids, proteins and other molecules, contributing to the membrane properties. Mercker et al. (2012) studied the curvature modulated sorting of lipids in such a membrane. In the absence of chemical interactions, no sorting is to be expected, but they showed that sorting of lipids occur if one of the macroscopically defined elastic moduli for either lipid or proteins were different. Surface curvature and properties influenced the membrane's structure. Increasing the difference of lipid and protein elastic moduli increased the stability of the formed membrane patterns. The membranes themselves were simulated by a finite element-based membrane dynamics approach. The macroscopic parameters were directly derived from a DPD simulation that was coupled with the FE solver. They concluded that in the absence of a multiscale approach, results could deteriorate and produce incorrect physical systems. Thus, an intermediate DPD solver was needed to feed the macroscopic solver with micro/mesoscale information in order to capture the correct membrane dynamics. A similar approach was adopted by Ghoufi and Malfreyt (2012) who studied the effect of salt concentration in water on its surface tension. They used a method called many-body dissipative particle dynamics (MDPD) which at its core is essentially the DPD method, but the inter-particle force, i.e. Eq. (44), includes not only the inter-particle interactions but also a local particle density. This poses the difficulties that model parameters are not readily available. Therefore, atomistic Monte Carlo simulations were performed for a fixed parameter space for which a linear relationship could be established. This relation was then used in their MDPD simulation to derive the required parameters for which they simulated various salt concentration for two salts (NaCl and NaF) and showed that with an increase in salt concentration, the surface tension increased correspondingly. Their trends were matched by experimental data. The approach is similar to that of Holland et al. (2014), discussed in Sect. 2.1. Fedosov and Karniadakis (2009) introduced their socalled triple-decker approach where they coupled MD, DPD and a Navier-Stokes solver. This allows a continuous description from microscopic through the mesoscopic regime up to the macroscopic layer. Not only is the length scale separation reduced via an intermediate, mesoscopic layer as seen in Sect. 2.1 where the microscale was directly coupled with the macroscale, but also this approach benefited from the fact that the mesoscopic solver, here using DPD, was particle based which made the coupling to the microscopic level simpler. On the other hand, the increase in time and length scale of the DPD method also meant that statistical scatter, in comparison with its MD counterpart, was reduced which made it easier to couple to a macroscopic solver. In their approach, MD particles were removed and inserted using the USHER algorithm (Delgado-Buscalioni and Coveney 2003a, b) and DPD particles randomly inserted near the boundary. MD boundary conditions were enforced by specular walls, pressure forces to reduce density fluctuations and shear forces for the tangential velocity component. They applied their approach to low Reynolds number for Couette, Poiseuille and lid-driven cavity flows using a domain decomposition approach. Comparison against a full Navier-Stokes solution showed the applicability of the method. From a qualitative point of view, statistical scatter was low and profiles at the HSI were smooth for all cases. Reducing the HSI to zero, i.e. exchanging information directly at an interface without overlap, was also investigated for the Couette flow. They showed that their multiscale approach started to deviate from the analytical solution and showed a maximum discrepancy of about 15 % at the interface which subsequently reduced towards the channel walls.
Reducing the length scales even further, Kacar et al. (2010) also used a three-layer interaction model, here coupling quantum mechanics (QM) with MD and DPD. QM was employed to find the atomistic interaction energies which in turn were used on the molecular level to obtain the parameters for the DPD solver, specifically the solubility and cohesive energy densities. The results obtained on the mesoscopic level were then mapped back onto the microscopic level which completed one coupling cycle. They investigated styrene-co-fluorinated acrylate oligomers and found that an increase in oligomer concentration caused the system to self-assemble from spherical micelles to hexagonal cylinder-shaped structures. The use of a multiscale approach made it possible to reveal structures on the mesoscale that were influenced by the molecular motion on much smaller scales. Furthermore, the study showed that in certain cases a MD description might not be enough and information from even smaller scales are needed as quantum effects may have an impact on the solution at much larger scales, as seen here on the mesoscale. Dzwinel et al. (2002) also investigated the self-assembly of micelles into hexagonal structures and further the Rayleigh-Taylor instabilities using a three-way coupling approach. Here, MD and DPD were used as well as the fluid particle model (FPM) which was introduced by Español (1998). The FPM is similar to DPD but further accounts for a non-central (shear) force. This is necessary as an orbiting particle around a central particle would not feel Page 21 of 38 68 any drag, which is eliminated with the non-central force in the FPM and so conserves angular momentum. As the FPM is slightly more complex, it is best suited for length scales above those of DPD so that FPM particles are large enough to only influence their immediate neighbours. Similar to Kacar et al. (2010), Dzwinel et al. (2002) found that the self-assembly of colloidal particles into hexagonal, worm-like structures depends on its interaction parameters (here, the Lennard-Jones interaction potentials). Furthermore, they investigated the cluster growth rate and found a power law behaviour where the exponent initially took values of 1/2, whereas for longer times a transition towards an exponent of unity was observed which was consistent with theoretical predictions. In the case of the Rayleigh-Taylor instability, they showed that cluster fragmentation occurred over time where larger clusters separated into smaller ones while the average cluster size increased over time, which is characteristic for a mixing problem. Li et al. (2013) used DPD and constructed a multiscale approach containing MD and a Monte Carlo-like polymerising model to study the effect of crosslink reactions and diffusion on the formation of carbon fibre and epoxy resin. They constructed a simple rectangular domain with a carbon surface at the bottom and periodic boundary conditions on the side. Initially, the hardener and epoxy resin were grouped together, away from the carbon surface and separated by a sizing agent. The simulation was then carried out at a constant temperature of 423 K and results presented at 60 min of simulation time. All three components, hardener, resin and sizing agent diffused into their neighbouring domain and combined into epoxy. The hardener was found to be chemically inadequate near the surface to form epoxy which was accompanied by a drop in crosslink density, compared to the bulk matrix density. The interface region, separated by the matrix and carbon fibres, was found to agree with experimental data. DPD presents a valid approach for studying larger atomistic groups but has shown to be inadequate for epoxy as greater molecular details were needed, and hence, MD was used in conjunction with a polymerising model to adequately represent such components.
A vesicle is a fluid-filled cell constructed from bilayered lipid membranes as described previously. They are involved in metabolism, transport processes or chemical reactions, which take place inside the vesicle. Sevink et al. (2013) used a DPD description to study vesicles with a radius of approximately 0.45 nm containing O(10 4 ) lipids for a duration of O(10 6 ) time steps. DPD was used to resolve the membrane structure and lipid interactions, conserving the necessary molecular details. The solvent surrounding the membrane was resolved by a combination of Brownian dynamics (BD) and dynamic density functional theory (DDFT) which allowed the solvent to be purely expressed by a continuous concentration variable. As a proof of concept, they tested their approach by initially placing lipids in a water solution with two lipid volume fractions. Lipids of approximately 5 v% transitioned into small aggregates which, due to the low percentage of lipids, remained at their initial small size. Increasing the volume fraction to approximately 23 v%, the results were substantially different. The lipids formed a discontinuous network, i.e. a membrane-like structure, but were prevented to form a fully closed membrane. Finally, they investigated a fully pre-assembled vesicle and removed a cylindrical portion of lipids at the membrane wall. After equilibration, the simulation showed that vesicle healing occurred by means of lipid reordering in the membrane to achieve a lower-energy state associated with a closed membrane. They concluded that the multiscale approach was needed to remove the unnecessary solvent details, which could be represented by a continuous field, in order to simulate realistic dimensions. Driven by the same motivation, Berezkin and Kudryavtsev (2013) studied reaction-diffusion flows. Their approach consisted of DPD which accounted for all the molecular reactions taking place, as a result of diffusion processes of initially separated particles or polymers. DPD was used and covered the entirety of the first particle/ polymer domain and extended into the second particle/polymer domain, covering the interface. From there, a simple finite element-based diffusion equation was solved for the concentration in the form of c t = Dc xx where D is the diffusion coefficient and c is the concentration field variable. The coupling, including an overlap region, was achieved by matching the continuum concentration field with the number density of the DPD particles. They validated their approach for a single reaction of different particles for which an analytical solution for the concentration/number density can be found. The FE solution showed a smooth profile throughout the domain and overlap region, while the DPD number density was subject to small fluctuations, however, overall matching the analytical solution. They then tested their approach on reactive coupling of immiscible polymer melts and interfacial polymerisation of two different monomers. In the first case, the number of copolymers per unit area n c at the interface was observed to increase initially with time as n c ∼ t and then slowed down to scale with the square root of time as n c ∼ t 1/2 . The fast built up of n c is explained in the reaction process taking place, while once the reaction had finished, diffusion became dominant and reduced the scaling behaviour. The same trends were observed for the second test case, in which monomer conversion into polymers exhibit the same time scaling law, as did the number-averaged and weight-averaged degree of polymerisation. Their simple coupling procedure allowed for greater computational efficiency by removing unnecessary molecular details in favour of a macroscopic concentration field description which is solely governed by diffusion.
A network of polymers can exhibit a loaded stage when external shear or normal forces act on the system. This behaviour can be observed in tissue engineering, paper manufacturing or biological systems. Masoud and Alexeev (2010) investigated such loaded polymer networks to study the permeation and hindered diffusion using DPD where polymer dynamics were obtained by a bond-bending lattice spring model. They validated their approach against experimental data for a non-loaded network for which they matched the data for permeability and diffusivity over a range of porosity values. For loaded networks, they found that the permeability and diffusion along the parallel direction to the normal forces experienced an increased rate compared to its non-loaded stage, while the components perpendicular to the force direction decreased. Under shear, greater permeability was observed for the components in the shear plane, i.e. for a shear force γ zx , the components in x and z were increased, while it had no influence on the permeability normal to the shear plane, here the y direction. The diffusivity only increased for the x component, while both y and z direction exhibit a decrease in diffusion. To characterise the degree of network alignment, they introduced a second-order tensor for which an isotropic network (no loading) showed each main diagonal component to equal 1/3 of the tensor's trace, while the off-diagonal components were zero (and nonzero for an anisotropic network). Normal forces only altered the main diagonal components while shearing introduced nonzero, off-diagonal components. Aligning the main axis of deformation with the coordinate system by transformation revealed a quasilinear master curve onto which all main diagonal components fell. Fedosov et al. (2012) considered only loaded networks under shear and investigated structural and rheological properties of star polymers for various Weissenberg numbers. These consist of f identical polymers which all connect to a centre node. The star polymers were modelled using a simple spring potential, and the excluded volume interactions were approximated by the shifted Lennard-Jones potential. The coupling occurred through forces in the DPD equations. They considered two hybrid approaches: one of which coupled DPD with the polymer chains and another using multiparticle collision dynamics (MCD). We will focus on the DPD results here although both MCD and DPD produced very similar results. The relative deformation of the star polymer showed to scale with Wi 2 up until Wi ≃ 1, after which it reached a saturated state and reduced. The alignment with the flow scaled initially with Wi −1 and then reduced to Wi −0.43 for Wi greater than unity since these polymers, for high shear rates, were able to align with the flow. Interestingly, the shear rate only scaled as Wi 1 until Wi ≃ 1 after which a weaker scaling was observed, hinting that some shear energy at higher rates was transformed. All of the above findings were tested for five different concentration levels and two different number of nodes making up the polymer; however, they all fell onto one single master curve.
We saw in Eq. (44) that a maximum repulsion force needs to be provided in order to calculate the conservative force. This parameter is a microscopic dependent value and for complex systems best derived directly from MD simulations. This was done, for example, by Maly et al. (2008) who studied the morphology of self-assembling nanoparticles in copolymer mixtures, which are made up only of two different monomers. When mixing polymers, the energy associated with the mixing process is measured by the dimensionless Flory-Huggins interaction parameter χ AB . For lamellar polymers, that is, layers of pure polymers stacked on top of each other, and the parameter a ij of Eq. (44) can be assumed to be a AA = a BB , where we have two polymers of types A and B. The solubility parameters δ i for both type A and B were obtained from MD simulations. With that knowledge, the Flory-Huggins parameter could be obtained as where V DPD is the volume of one segment of a DPD bead. This led to an empirical equation for a AB as a AB = a AA + 3.27(1 + (3.9/N 0.51 DPD ))χ AB . We refer the reader to Maly et al. (2008) for a detailed derivation on how to obtain the Flory-Huggins parameter and briefly discuss their findings (see also Fermeglia and Pricl 2007) for a general description on how to obtain the Flory-Huggins parameter from MD simulations). They constructed nanoparticles with a central, neutral DPD particle covered by 12 DPD particles of different combinations of A and B. They showed that the introduction of these nanoparticles initiated their self-assembly into distinct cylindrical shapes for various combinations of type A and B, covering the central DPD particle. Scocchi et al. (2007) investigated polymer-clay nanocomposites (PCNs) and used a similar multiscale approach by coupling MD to their DPD solver. PCN is a material that exhibit increased mechanical strength, heat resistance and reduced permeability for gases in comparison with a pure polymer matrix while having only a fraction of the weight. Using conventional macroscopic FE solvers, PCNs are of interest to material engineers but depend on their microscopic interactions, which are not fully understood. Critical parameters remain uncertain, making mesoscopic solvers a prime candidate for which the critical parameters are derived directly from MD simulations. They studied and compared the structure of a Nylon-6-MMT PCN to published data and found close agreement with expected results, paving the way for Toth et al. (2009) who studied different PCNs using MD and DPD and further added a macroscale FE solver. Microscopic details were preserved and passed from the lowest scale through the mesoscale up to the macroscale, and thus, no model constants needed to be defined. Careful validation of each stage showed the applicability of MD and DPD at the micro-and mesoscale. Results on the macroscale were then used to derive parameters of macroscopic interest, here the coefficients of thermal expansion and electrical conductivity. Four different polymer sizes in the PCNs were tested and showed that the thermal expansion coefficient decreased for all PCNs over the clay loading, while the rate of decrease was the steepest for the smallest polymer. The electrical conductivity, however, showed an increase over the clay loading, independent of the polymer size. The importance of this study was the elimination of model constants at the cost of a triple-decker approach. A trade-off in accuracy and computational speed has always to be found and assessed for each case individually. However, where accuracy is the sole motivation, this approach presents an excellent coupling procedure that involves all physical scales.
Biological and medical flow applications are often dominated by mesoscopic cell transport, such as red blood cells and vesicles. Filipovic et al. (2012) simulated different particle flow scenarios using a multiscale solver consisting of a macroscopic FE, and a DPD or LBM solver at the mesoscale. Their influence was introduced at the macroscale through forces at the mesh nodes which were then incorporated into the FE equations. They considered four different test cases: (1) one cylinder and (2) two cylinders depositing in a channel, (3) four elliptical particle subjected to shear flow and (4) particles in an arterial geometry of different shapes. For the first case, the drag force exerted on the cylinder agreed well with Oseen's approximation for both FE-DPD and FE-LBM. Comparison with experimental data of the cylinder displacement showed an initial, undisturbed settlement which was gradually disrupted by harmonic oscillation due to the vortex street forming at higher Reynolds numbers. Results for the second case were only presented qualitatively. For the shear flow, the FE-DPD results were compared to a pure FE solution and the trends were matched, and finally, the arterial geometry, essentially a Y-channel, was qualitatively simulated for cylindrical and elliptical particles. In their proof of concept study, they presented simplified flow geometries (often at the core of biomedical flows) and showed that mesoscopic quantities could successfully be projected onto the macroscopic domain. In fact, not in the spirit of multiscale methods, yet important to our discussion, Feng et al. (2012) investigated arterial (channel) flows with different values of stenosis (channel contraction) using independently a Navier-Stokes-based macroscopic solver and DPD. The transversal velocity profiles at different axial locations agreed well and both Navier-Stokes and DPD matched experimental data. The authors also investigated the recirculation zone just behind the stenosis, caused by the quasi sudden expansion. Smooth results were obtained from the Navier-Stokes solution while DPD suffered from larger statistical scatter, although still reproducing a distinct recirculation zone. The velocity magnitude was integrated in a plane covering the whole recirculation area for different stenosis values. For a stenosis of up to 40 %, the vorticity magnitude differed by a factor of two and only at higher values approached each other. The agreement at 75 % stenosis was within a few per cent. The standard deviation, due to the statistical scatter, was also higher for the DPD solution than for the Navier-Stokes part. The integrated backflow for 64 and 75 % stenosis showed a difference by a factor of 4 and 1.5, respectively, which was explained by the fact that the density immediately behind the stenosis varied while the incompressible Navier-Stokes solution assumed a constant density flow. Their study showed that DPD was able to capture the trend of the Navier-Stokes equation while being prone to statistical scatter. Therefore, it can be concluded that DPD possesses the ability to reproduce continuum effects while containing a larger level of noise compared to the Navier-Stokes equations. For a smooth solution, a hybrid or multiscale method is expected to work best, while details on the mesoscale are retained.
Platelets are small parts which are transported in the blood. They are smaller than red or white blood cells and are responsible to stop bleeding by coagulating (clotting) around wounds which on a macroscopic scale can be interpreted as changing from a liquid to a solid state. This process is referred to as haemostasis, and platelets are said to form a haemostatic plug. Zhang et al. (2014) investigated the platelet dynamics using DPD to solve for the blood and coarse-grained MD to study the platelets. The coupling was done by forces. Platelets needed to get activated before they attached to an opening in the endothelium (the inner most layer in the blood vessel) which was done via the platelet stresses. Therefore, the authors investigated rigid and deformable platelets and compared their shear stress distribution on the surface. The particles were subjected to a shear flow with moving walls applied to both upper and lower wall, equal in magnitude but opposite in direction. Initially, they compared the angular velocity and angular acceleration for which differences as high as a factor of 2 were observed. The average induced stresses on the membrane surface were 2.5 times higher for the rigid body than for the deformable one. Since the rigid body would not yield, it resisted the flow and hence higher values were possible. The stresses were also faster accumulated on the rigid platelet than on the deformable one. Their study showed the importance of accurately describing the platelets as the reduced shear stresses on the surfaces not only altered their flow behaviour subjected to shear flows but also suggested that rigid models, overestimating the shear stresses, yielded incorrect activation potentials which can have a significant impact on the correct modelling of the haemostasis. Tosenberger et al. (2013) investigated the mechanism of clot growth (coagulation) as the mechanism is still not fully understood and current theories were rendered inaccurate. In a first simulation, they only modelled the platelets in a blood flow where inter-platelet connections strengthened over time, while platelets already attached for a sufficient amount of time were not able to attach new ones. In this model, the essential process of clot growth was reproduced where a substance known as fibrin is implicitly included as the connections between platelets were growing over time. Fibrin is a protein that forms the haemostatic plug with the platelets. It can be regarded as the glue holding the platelets together. They showed that the obstacle formed by the clot disturbed the velocity distribution, while for large aggregates, the weak forces were not able to resist the drag caused by the platelet clot and ruptured, leaving a smaller platelet clot of only strong connections. In a second investigation, DPD was coupled with a reaction-diffusion-advection equation which was used to calculate the fibrin concentration explicitly. The second simulation was similar to the first one, in which platelets were not able to attach to other platelets in the clot if their fibrin concentration exceeded a threshold value, irrespective of the time they were in the clot. The mechanism of coagulation was essentially the same, while the added level of information available through the fibrin concentration field revealed a constant growing platelet clot. Ruptures of aggregates were still possible, however, a much finer distinction could be made between weak and strong connections. Furthermore, it was possible to show whether the platelets were covered by fibrin or not. Kojic et al. (2008) coupled DPD with a finite element representation of the Navier-Stokes equations. Here, DPD cells were locally introduced in the FE mesh and coincided in a way that the FE mesh edges were the boundary of the DPD region, which could span over several FE cells. This approach resembled the heterogeneous multiscale approach (HMM) as described in Sect. 2.1, where mesoscopic details are only locally computed. The FE solver covered the whole domain and the locally introduced DPD region was coupled to it by forces. This allowed for further computational savings by only considering mesoscopic details where necessary. They applied their approach for two simple benchmark cases: the channel and lid-driven cavity flow. For the channel, the DPD region was placed inside the centre of the channel, without contact to the wall. Compared to an unsteady channel flow, the multiscale approach matched the results of the analytical solution. For the lid-driven cavity, the top right corner (in this case the lid was moving from the right to the left) was presented by the DPD region and velocity profiles compared to the pure FE solution showed very good agreement. In fact, DPD provided very smooth velocity profiles for both cases, while no extra smoothing apart from averaging was done. It is interesting to note that the Reynolds numbers were of the order of 10 1 −10 2 for which high statistical noise would be expected. Even spatial and temporal regression as discussed in Sect. 2.1 were not able to achieve similar level of noise reduction. This approach was recommended for use to study platelet dynamics as the clot region can effectively be modelled by DPD (and further coupled to a microscopic description if needed), while the vast majority of the blood can be modelled with an inexpensive continuum solver. In this approach, one would simplify the computations in favour of computational resources. The opposite approach would be to retain all molecular details in favour of high-fidelity results. For systems of realistic real-world sizes, however, simulations and high-performance computing (HPC) must go hand in hand. For that matter, Grinberg et al. (2011) investigated blood flow in a brain with an aneurysm (a sack-shaped expanded blood vessel). The geometry came from an MRI scan, and a coupled DPD/Navier-Stokes solver was developed from open-source tools. LAMMPS was used for the DPD simulation, calculating the blood plasma, red blood cells and platelets, while NekTar (http://www.nektar.info/, see also Cantwell et al. (2015) for NekTar++, its successor), a spectral element solver, was used on the continuum side for the blood flow. The mesoscopic side contained 10 11 particles and the continuum solver 3 × 10 9 unknowns. In their study, they used up to 2 × 10 6 CPU cores on two different clusters, namely a Blue Gene/P and Cray XT5 cluster. For both systems, the strong scaling followed a linear trend up to about 10 6 CPUs. Only at 186,624 cores, run on the Cray XT5, the strong scaling dropped to 68 % efficiency. Their work was important towards efficient code design and coupling. Due to its general nature, it is equally applicable for a broader range of applications and is not limited to medical flows.
Mukhopadhyay and Abraham (2009) coupled a DPD and MD solver to study the flow over a square obstacle in a two-dimensional channel. The geometry was the same as studied by Nie et al. (2004). The domain was split in two parts: the bottom one containing the square obstacle was modelled using MD and the bulk flow above using DPD. While their validation showed good agreement to an analytical solution (Poiseuille flow without obstacle), they revealed further insight into the efficiency of multiscale methods (in this case for the Poiseuille flow). There were essentially two mechanisms that can increase the computational time: coarse graining, i.e. using different length scales and thus fewer particles (as is the case for DPD particles compared to MD particles) and considering timescale separation, for which fewer integrations need to be performed in order to advance to a given time level. They showed that the speed-up to be expected from only considering length scale separation scales with O{[((N m } where the ratio of particle masses of each respective simulation is given by N m = m DPD /m MD , and κ is the ratio of the volume of the DPD and MD region. Further introducing timescale separation yields O{[((N m √ τ + κ)/(1 + κ)N m √ τ )] 2 }, where τ is the relative timescale between DPD and MD. For different time steps of MD and DPD, i.e. t DPD > t MD , the numerator is growing faster than the denominator and further computational savings can be achieved.

Interim conclusion on hybrid and multiscale dissipative particle dynamics methods
We have shown some examples of hybrid and multiscale DPD applications. Due to its mesoscale nature, it can provide coarse-grained atomistic details for a macroscopic description or speed up a microscopic solver by extending the length scale. DPD is a versatile particle approach that can handle molecular details but provides less statistical noise as its MD counterpart. We have seen that coupling between DPD to either the microscopic or macroscopic level can be achieved by exchanging forces between the solvers. Furthermore, while coupling to the macroscale produced less fluctuations, coupling to the microscale was easily achieved due to the similarities of DPD and MD. Hence, the triple-decker approach used here by various authors coupled successfully the microscale with the macroscale, without experimental data input (i.e. the interaction parameter a ij for the DPD solution) and provided a way to couple complex systems.

Smoothed-particle hydrodynamics method
Smoothed-particle hydrodynamics (SPH), originally introduced and used for astrophysical applications (Gingold and Monaghan 1977;Lucy 1977), is a particle-based, meshless approach for the solution of the Navier-Stokes equations. The underlying approach is the integral representation method that in steps is approximated to yield the final form of the governing equations of the SPH framework. We will briefly explain the key equations and concepts, while full derivations can be found in Liu and Liu (2003). We start with a trivial mathematical expression which is the underlying core of SPH. We can recast any function into an integral form as where the Dirac delta function is defined as As long as f (x) is continuous in Ω, Eq. (50) is an exact representation. In the SPH framework, the Dirac delta function is replaced by a smoothing function, also referred to as the kernel, which yields is the kernel with a support domain h, which is the region where the kernel is active and can be freely chosen. It needs to fulfil the following requirements: 1. Its integrated value needs to be unity, i.e.
2. It needs to be positively defined, i.e.

It needs to have compact support
where κ is a scaling factor 4. It should approach the Dirac delta function in its limit, i.e. lim h→0 W (x − x ′ , h) = δ(x − x ′ ) 5. It needs to be sufficiently smooth, be symmetric about its origin and decay away with increasing distance.
A Gaussian kernel may be used but only approaches zero at infinity. Research has been dedicated to find appropriate kernel functions; for example, the B-spline function is defined as where R = |x − x ′ |/h and α depend on the computational dimension with values of 1/h, 15/7π h 2 and 3/2πh 3 for one, two and three dimensions, respectively.
Note that in Eq. (52) we follow the SPH convention and used the kernel approximation operator 〈 〉 to indicate that we approximate Eq. (50) with a kernel. We make further use of a particle approximation in which we replace the integral by a summation over all particles within Ω as where N is the number of particles within h and a shorthand notation can be introduced as W ij = W (x i − x j , h). By particle approximation, we mean that we fill the domain Ω with a sufficient amount of particles and solve Eq. (54) for each one of them. The integral Ω dx ′ reduces to the volume of each particle V j and is replaced by its mass and density in Eq. (54). The first derivative is found by differentiating both sides. Applying the Gauss theorem and dropping terms which are zero, one arrives at We see from Eq. (55) that in order to find the derivative of f (x i ), we need to find the derivative of W ij which is analytically known and therefore readily available.
To this point, we have only approximated Eq. (50) by particles which is applicable to a general set of differential equations. We apply this methodology now to the Navier-Stokes equations for which we first write them in Lagrangian coordinates. Denoting the spatial coordinates by α and β, we have the continuity equation as the momentum equation (without external forces) as and finally the energy equation as The total stress tensor is decomposed into its pressure and viscous contribution as We now combine Eqs. (56-58) with Eqs. (54) and (55) to arrive at the final equations for the SPH method. The full mathematical derivation is omitted but can be found in Liu and Liu (2003).
The continuity equation becomes for the momentum equation we have and for the energy equation we obtain Boundary conditions are not as easily implemented in the SPH framework as in other particle-based methods. Similar (59) σ αβ = −pδ αβ + τ αβ . (60) as with MD in which particles interact with their neighbours, SPH particles interact with those in their support domain. There are no particles outside of the boundaries, and hence, the kernel will be truncated. The idea is similar to MD to place other particles outside of the boundaries which are also called ghost particles in the SPH framework. However, these ghost particles do not possess a well-defined physical property and need to be modelled via any reasonable physical assumption about their nature. Liu and Liu (2003) reviewed and presented a novel technique by combining existing treatments for solid boundaries. They introduced two new sets of ghost particles, the first of which is residing on the solid boundary itself. These particles repel particles within their vicinity back into the computational domain via a repulsive force. The second set of particles is residing outside of the boundaries and ensures that the support domain will not be truncated for particles near the solid boundary. They are mirrored particles of the real particles near the boundary and have the same physical properties except for the velocity component which is opposite in the boundary normal direction. Keat et al. (2015) reviewed various open boundary conditions and showed that there is still active research in this field due to the difficult nature of imposing correct open boundary conditions. Federico et al. (2012) provided an easy and intuitive approach in which two new sets of ghost particles were introduced, resulting in a total of four particle sets, i.e. fluid, ghost, inflow and outflow particles. Each of these inflow and outflow particles are placed at the open boundaries and have several layers, just as the ghost particles, to prevent a truncated support domain. The inflow particles are assigned with the desired pressure and velocity, while the physical properties for the outflow particles are frozen. The inflow particle are advected in the same manner as the fluid particles and change their set from inflow to fluid once they cross the open boundary. New particles are inserted into the inflow particle set to conserve mass. Similarly, the outflow particles are advected with their frozen velocity until they reach the end of the outflow domain and are then deleted. Each particle evolves with the equations given above for which standard CFD time stepping procedures can be used. Note that due to its Lagrangian particle nature, SPH represents a meshfree method for solving the Navier-Stokes equations and thus particles can easily be refined in regions of interest. This adaptive behaviour is easily done, in contrast to classical CFD techniques, where, for example, either structured meshes with hanging nodes or tetrahedral elements on unstructured grids need to be employed. This requires further solver modification which is not the case in the SPH method. Multi species transport is easily achieved by defining particles of different properties. However, it has been found that higher-order kernels may not conserve physical quantities as W ij may locally have negative values. For the approximation of non-negative quantities such as the density, this could yield nonphysical results. Stability, accuracy and convergence have been mainly tested for simple applications with equispaced particles in lower-dimensional space so that it is difficult to fully assess the strength and weaknesses of the method. We refer to Monaghan (2005) for a detailed discussion about the advantages and disadvantages as well as applications.

Review on hybrid and multiscale smoothed-particle hydrodynamics methods
The breakdown of the Navier-Stokes equations and flow phenomena at either very small length scales or rarefied gas conditions has fuelled research on particle methods. Inevitably, some applications can be described using several methods at comparable accuracy and computational cost. For a range of industrial applications, channel, pipe and flow through co-axial cylinders are often encountered. Filipovic et al. (2009) investigated those applications and compared DPD and SPH results independently to elucidate their strength and shortcomings in comparison. They investigated two laminar, low Reynolds number flow cases. The first was a simple Poiseuille flow and comparison with an analytical solution showed that there was little difference between the two methods, both showing accurate results. The second test case was the flow through a co-axial cylinder for three different conditions. The first set-up saw both cylinders rotate at the same angular velocity and the second only the outer cylinder, and for the third case, both cylinders were rotating with the same angular velocity but in opposite directions. The agreement with the analytical velocity profile was good for both methods and the tabulated error revealed that DPD, at best, was marginally more accurate than SPH. The density profile, on the other hand, showed variations between the cylinder walls for the DPD results which were absent in the SPH solution. They showed, however, that the computational time was about two orders of magnitude less for the SPH method which was due to the larger time step permissible as thermal fluctuations were not resolved as in the DPD method. Therefore, while maintaining a comparable level of accuracy, SPH is a fast method and may possess advantages features to be exploited by multiscale applications. Unlike the continuous description of the Navier-Stokes equations, SPH models the same physical process as a particle which makes it easier to couple it to other particle methods. Liu et al. (2002) coupled SPH directly to MD without the need for an HSI. This coupling, compared to direct Navier-Stokes/MD coupling, has several advantages. For instance, particles do not need to be inserted or deleted, and the initial SPH particle density can be matched at the interface and then smoothly decrease away from the interface. By further making SPH particles at the interface interact with MD particles through an interaction potential, MD particle does not feel a sudden loss of interaction potential at the interface, as on open boundaries, and hence, no correction needs to be done. Since SPH particles have an in-build smoothing ability by averaging over neighbouring particles, thermodynamic properties are constant across the interface. They tested a simple Poiseuille flow for which the wall regions were modelled by 1600 MD particles, while the centre part of the domain was covered with 600 SPH particles. The results matched the analytical solution, and indeed, a smooth profile across the contact interface could be observed. Jorn and Voth (2012) argued that MD is a very versatile approach but leads to disparate length scales for proton exchange membranes (PEM) which are found in fuel cells. Therefore, they developed a hybrid DPD-SPH approach in which DPD was modelling the water and polymers, while SPH accounted for the proton transport through the membrane. They validated their approach for a simple lamellar test case for which the proton concentration was fixed at both ends, i.e. high proton concentration at the anode (proton production) and low concentration at the cathode (proton consumption) and showed good agreement with an analytical solution for both concentration and proton flux profiles. They applied their hybrid model to a PEM geometry to study the conductivity behaviour with various water contents, which is the number of water molecules per sulphonate group. Their results revealed that explicitly accounting for the proton produced a linear increasing relationship between conductivity and water content, while excluding electrostatic effects yielded a nonlinear relationship. Furthermore, the conductivity was approximately two times greater for the highest water content when including the electrostatic effects, while for the lowest water content, no difference was noticeable which showed that electrostatic effects have an important influence on the conductivity. Experimental results for the conductivity were approximately a third higher and also showed a linear increase over the water content suggesting that proton exchange needs to be included for a correct physical representation. Murashima and Taniguchi (2010) used the modified smoothed-particle hydrodynamics (MSPH) method, introduced by Zhang and Batra (2004), and coupled it to a dumbbell model to obtain the stresses in a dilute polymeric fluid. MSPH is an extension of the conventional SPH method in which f (x i ) and its derivative is obtained via SPH and a set of equations, derived from a Taylor series expansion, are satisfied leading to more accurate results in the vicinity of boundaries compared to the classical SPH approach. The dumbbell model is a simple evolution equation of the end-to-end vector of a polymer chain, i.e. the distance between the end points of the polymer. The hybrid procedure was as follows: velocity, pressure and the new location of particles were calculated via MSPH and the divergence of the velocity then calculated. The dumbbell model was used to obtain the stresses and MSPH again to find its derivative. This procedure was repeated until t max was reached. For a simple channel geometry, they showed that if the relaxation parameter is chosen too small, thermal fluctuations become dominant. Increasing the relaxation parameter reduced this behaviour and qualitative agreement with a pure macroscopic solution showed the expected velocity profile. Yuan et al. (2009) investigated several approaches to simulate the deposition behaviour of metal particles in low temperature, high-velocity air fuel spraying. The surface and impinging particles were modelled either continuously with finite elements or particle based with SPH, and thus for surface-particle interaction, the following three methods were applied: FE-FE, FE-SPH and SPH-SPH. They found for increasing particle temperature and velocity that the interface temperature and plastic strain were locally affected by the impinging particles. Due to the high particle energy, upon impact, particles deformed and thus heated/strained the surface over time. The effective plastic strain was further found to be a function of the impact velocity as particles below a critical velocity remained mostly intact, while for higher velocities, scattering was observed. Comparison between the different simulation approaches, i.e. SPH and FE for particles and surface, showed similar behaviour. Tartakovsky et al. (2008) used the core idea of SPH, i.e. Eqs. (54) and (55), and applied that to a reaction-diffusion equation (RDE), which solves for a concentration field, to simulate flow through porous media. Two mixing solutes A and B were used which formed a precipitate product C. It was assumed that C only accumulates at pores, and thus, a domain decomposition strategy could be adopted in which one RDE was solved at the pore level, while an averaged version of the RDE was solved on the macroscopic level. Since precipitation growth only occurred at the pores, reactions only needed to be captured with the RDE at the pores, while reactions were absent in the macroscopic counterpart. Different levels of resolutions were therefore used to reduce the computational time. The coupling was easily achieved in a non-iterative manner as both descriptions solved the same underlying equation. Ghost particles were introduced that could accumulate mass over time. They were transformed into solid particles if their mass exceeded a threshold value, and thus, precipitation growth could occur. They validated the approach against analytical solutions for planar and circular precipitation growth and tested the multiresolution scheme against a single resolution simulation of a passive scalar in porous media. Finally, they tested their approach in a porous medium using the hybrid scheme which showed that the porosity decreased over time due to precipitation accumulation in the pores. They tested it against a full-scale simulation which matched their hybrid approach. The advantage of their study was that the interface did not need explicit tracking but rather was a product of the solution.
Apart from molecular and particle-based applications, SPH has also been used to investigate geophysical problems. In rock removal for mines or tunnel excavation, rock blasting has become an important alternative to drilling. Fakhimi and Lanari (2014) used a hybrid solver comprising of SPH for gas dynamics and discrete element method (DEM) to model rock specimens. Rock particles were represented by circular discs and attached to each other. The domain was circular and housed DEM rock particles in its entire domain except for the centre, at which SPH particle resided. After a simulated explosion occurred, SPH particles expanded into the rock material. Collision via SPH and DEM particles were accounted for by conserving momentum through ideal plastic collisions. Once the contact normal and shear stresses of DEM particles were exceeded, bond break ups and microcracks were formed. By further introducing a plane in a second simulation, located off-centre and slicing the circular domain into two parts at which the rock parameters were 1 % of that of the intact rock, simulations showed that the shockwave was not strong enough to cause further damage beyond the plane. Microcrack growth drastically reduced while the waves that had been reflected at the plane caused substantial damage to the rock material. Another geophysical example is that of a breaking dam. Hilton and Cleary (2012) constructed a multiscale method from SPH and the shallow water (SW) equations. Their initial study, to investigate each method on its own and compare it to experimental data and a multiscale method, was a simple dam break scenario. Here, water was initially held at a height h for a small portion of the domain, whereas the rest of the domain was left empty. Downstream of the water box was a pillar that interacted with the water. Velocity measurements were taking just upstream of the pillar and showed that all methods followed the experimental data. The pure SPH implementation was more computational expensive than the SW method so that the multiscale method was able to reduce the computational cost. Information exchanged between SPH and SW was the velocity and water height which were imposed at a single interface. With their multiscale approach validated, they applied it to a realistic scenario based on the Gordon dam in Australia spanning a region of approximately 71 × 62 km. They reported results for various intervals and showed that after a period of 80 min, a secondary lake formed due to flooding that spanned over 15 km. For such applications involving large length and timescales, a multiscale approach is necessary to capture the small and necessary details to accurately predict the outcome. The use of a multiscale method meant that, for a particle spacing of 15 m, the number of particles could be reduced by two orders of magnitudes (from 2 × 10 7 for a full SPH simulation to 1.5−2 × 10 5 for the multiscale simulation).
SPH has proven itself to be a very robust and fast method for visual applications and has been used in major movie productions (Monaghan 2005). Raveendran et al. (2011) noticed that due to the kernel approximation, sudden changes in pressure are not transferred instantaneously, as only the region within the kernel is updated. Hence, pressure changes manifest themselves as artificial acoustic waves which can be suppressed by increasing the stiffness or reducing the time step. They introduced a new hybrid approach where a coarse grid was superimposed onto the SPH domain and the pressure was solved via a Poison equation. In this way, the pressure could be obtained directly by interpolation, while a density correction procedure ensured that a constant density was kept throughout the domain. Their approach required 10 and 18 % more computational time for an incompressible and weakly compressible SPH solver; however, since numerical artefacts such as the spurious acoustic waves were reduced, much greater time steps could be employed, overall increasing the efficiency of the proposed algorithm. They applied their solver to several different water splashing cases using 10 5 to 5 × 10 5 particles and showed that their speed-up was 3-4 times compared to a standard SPH algorithm. Despite the visual effects orientated approach, the inclusion of a Poison solver on a coarser grid represented an elegant way to include the elliptic nature of pressure in the SPH method. Wang et al. (2013) followed a different approach for visual applications in which the fluid phase was solved by a fast LBM-based solver. The emphasis has been on the free surface dynamics, and hence, the use of LBM removed the problems encountered in Raveendran et al. (2011). The free surface front was tracked with massless particles where a coupling band existed that separated the fluid from the gas phase. Its centreline coincided with the free surface and expanded a distance d into both phases. Particles were inserted into the coupling band and were simulated via SPH. The velocity was imposed onto the particles via a weighted interpolation, gradually blending from a pure LBM velocity to a pure SPH velocity. The purpose of the added SPH particles was for visual effects only. Particles penetrating into the LBM domain from the coupling band were rendered as bubbles, while particles inside the coupling band were rendered as foam. Particles escaping the coupling band, outside the fluid domain, were subsequently rendered as spray. Their hybrid solver produced similar results compared to other visual approaches for various scenarios, while the computational speed, owing to the fast LBM solver, was among the fastest.
We now turn our attention to a new, hybrid method that has emerged in the last decade and has successfully been applied to various cases. In the introduction to this section, we have mentioned the work of Filipovic et al. (2009) and showed the close relation of SPH and DPD. Español and Revenga (2003) introduced a new method based on both SPH and DPD to harvest their respective advantages and termed it the smoothed dissipative particle dynamics (SDPD). This hybrid method reproduces the behaviour of the Navier-Stokes equations (SPH part), while thermal fluctuations on the mesoscopic level are retained (DPD part). Hence, the new method is able to be used for a broader range of scale separation, whereas no explicit coupling is needed in-between methods. The introduction of SPH into the DPD framework also means that SDPD particles have a physical dimension, unlike the arbitrary size of DPD particles. Vázquez-Quesada et al. (2009) investigated the scaling behaviour of the SDPD method for particles and polymers undergoing Brownian motion in suspension. In the case of a single colloidal particle (constant volume) surrounded by solvent particles that have been tested for different solvent particle volumes (by refinement), the thermal fluctuation of the colloidal particle has been shown to be scale invariant (over the solvent particle volume), while the thermal fluctuation of solvent particles showed a square root dependency over its volume. This behaviour is to be expected as the thermal fluctuation should only depend on the size of the particle and not on that of the surrounding ones. This mechanism is not observed in DPD as no scale information is available. Alternatively, either MD has to be employed to derive accurate microscopic parameter that can be imposed on the mesoscopic level or a trial and error approach with fine tuning parameters a posteriori. These are two common choices to circumvent this shortcoming of the DPD method.
Due to the DPD nature of the SDPD method, it can also be applied to microscale flows. Litvinov et al. (2008) investigated a suspended polymer chain in solution where the solvent and polymer molecules were modelled using SDPD, while polymer beads were bonded by a FENE spring model. First, they tested a two-dimensional polymer with zero imposed velocity and studied the stochastic motion of the polymer. The radius of gyration, which is the trace of the gyration tensor, increased linearly with the number of beads as expected from the theory. Here, the gyration tensor is a measure of how far the molecules of the polymer are apart, i.e. G = 1/2N 2 i,j r ij r ij where N is the number of beads. In a second simulation, they applied their method to a polymer in a microchannel and altered the channel height to study the behaviour of confinement. Decomposing the radius of gyration to take only the x and y components of the polymer chain into account, i.e. the radius of gyration parallel and normal to the channel wall, showed anisotropic behaviour for channel heights H < 10 where H was normalised by the unconfined radius of gyration. For H > 10, both parallel and normal component collapsed onto one single curve, showing the diminishing effect of channel confinement and return to isotropic behaviour. In a similar study, Bian et al. (2012) investigated a single particle in suspension and applied it to several test cases ranging from typical examples found in microto macroscale flows which included flow through porous media, a particle with initially imposed velocity (kick started), particles under shear, interaction of two approaching spheres and a neutrally buoyant particle under Brownian motion close to and sufficiently far away from a solid boundary. Their test cases agreed well with reference data. In their last test case of a particle suspended near the solid surface, they also investigated the anisotropic behaviour of a single particle based on its diffusion coefficient. They decomposed it into diffusion parallel and normal to the wall, but unlike in Litvinov et al. (2008), they changed the position of the particle from the wall, while the second wall on the other side was placed sufficiently far away so that only the influence of the closest wall was felt by the particle. They observed that the diffusion parallel to the wall was always greater than normal to it, independent of the distance from the wall. When placed directly on the wall, the normal component was zero, while diffusion parallel to the wall was still permissible. When placed further away, the diffusion coefficient approached each other. Comparison with theoretical results showed that close to the wall (up to a couple of particle radii), SDPD and theory showed discrepancies of up to 30 % which vanished with increasing distance away from the wall. Hence despite the loss in accuracy close to the wall, the general behaviour of the microscale (anisotropic behaviour) was captured in both approaches. Accuracy close to the wall may be increased by using MD for near-wall regions and SDPD for the rest of the domain. Litvinov et al. (2009) derived an analytic expression for the diffusion coefficient to be used in conjunction with SDPD liquids. An accurate knowledge of the diffusion coefficient a priori is important for cases where the Schmidt number Sc influences the flow, as, for example, in the study of non-equilibrium polymer properties. They applied it to a periodic box and showed that flow properties were accurately captured. Furthermore, they showed that no secondary structures were formed during the simulation (crystallisation) by showing the mean square displacement to monotonically increase without plateauing (indicating solidification, decrease in diffusion). The radial distribution function followed the expected behaviour without secondary peaks, which would have indicated particle clustering.
SDPD in itself can be regarded as a hybrid approach combining SPH and DPD. However, it can be further coupled with a microscopic or continuum solver to further enhance details or gain computational efficiency. Moreno et al. (2013) presented a preliminary multiscale approach in which finite elements for the macroscopic Navier-Stokes equations were coupled with SDPD for biomedical applications. Fine scales were locally resolved with SDPD, while the majority of the domain could be represented by the FE method. Although a two-way coupling has been implemented, the coarse scale has been assumed to be steady as the time to reach equilibrium on the smallest scale was less than the macroscopic time step, i.e. t fine ≪ t eq ≪ t coarse . Their test case was the lid-driven cavity, and they showed that the horizontal velocity profile of coarse and coupled velocity matched well. The extension to non-Newtonian fluids is necessary to study true biomedical applications such as blood flow which has been left for future work. Kulkarni et al. (2013) further investigated the scaling behaviour of the SDPD method and introduced a novel multiresolution scheme. First, they reported on the equilibrium properties of SDPD particles in a three-dimensional box. They showed that for different particle resolutions, the density stayed constant throughout the domain and the temperature stayed within 2 % of its imposed value. The velocity PDFs of the different particle systems showed that with increasing the particle resolution, the Maxwell-Boltzmann distribution showed narrower tails and higher peaks, approaching the Dirac delta function in its limit. This, again, showed the correct scaling behaviour of SDPD. They devised a multiresolution scheme in which particles were able to change their resolution from one domain to another. For example, if a particle was refined by a factor of 2, once it entered the refined domain, it would be deleted and two new particles inserted. To conserve mass and momentum, both would need to contain half the mass of the original particle while sharing the same, original velocity. One has to note that boundary conditions are sensitive to scale changes and, for such an approach, need to be considered carefully. They applied the multiresolution approach to an equilibrium test case showing that while the particle resolution was changed smoothly, the density was not adversely affected in the refinement region and showed a constant profile across the domain. The second test case was the Couette flow for which a fine resolution was chosen near the stationary, bottom wall, while a coarser representation was chosen at the moving, upper wall. The resolution was changed at the channel centreline in a region of finite thickness. They were able to match the analytical solution in this way. The ability to smoothly change between resolutions, while SDPD is applicable over a range of length scales, makes it a powerful and computational efficient method. Petsev et al. (2015) extended the idea of Kulkarni et al. to a multiscale, multiresolution SDPD-MD algorithm. First, they derived a multiscale approach based on the AdResS scheme of Praprotnik et al. (2005), see Sect. 2.1 for a detailed discussion, in which the pairwise force evaluation was blended between both description for a smooth Page 31 of 38 68 change. Applied to an equilibrium test case, they showed that the temperature and density profile across the domain differed by 1.8 and 3.2 % at most, where the domain was divided into MD and SDPD with a buffer in-between. They applied it to the Couette flow where the moving wall was set along the MD domain in the first set-up and along all three domains (MD, SDPD and buffer layer) in the second set-up. For both cases, the correct start-up behaviour was observed with minor statistical scatter. In another simulation, they divided the SDPD region into a fine and coarse region where the fine region was located near the MD interface to match its resolution. The moving wall was only applied to the MD domain, and again a good comparison with the analytical solution was observed with minor noise. Although the resolution was successfully changed within the SDPD domain, decreasing the computational demand, the same time step had to be used across all domains which removed some of the potential computational savings. A remedy could be the use of multiple timescale integrators which were not investigated. However, the ease of blending between MD and SDPD via the AdResS scheme and the multiresolution character of the SDPD method showed great potential of the SDPD method to be accurately used on the smallest scale, while upscaling showed the results to converge to the Navier-Stokes equations.

Interim conclusion on hybrid and multiscale smoothed-particle hydrodynamics methods
We have briefly discussed the progress in hybrid and multiscale SPH methods. Its drawback in accuracy has favoured other comparable methods; however, its ease of implementation and computational efficiency has led to applications such as visual effects. However, SPH remains a pure particle-based method that adheres to the Navier-Stokes equations ] and hence has shown superior coupling properties to other particle-based method, most prominently MD and DPD. The recent introduction of SDPD, a combination of the favour-able properties of SPH and DPD, provided new research interest into the SPH method, especially for the mesoscopic domain. It is, however, not limited to this regime as microscopic details are captured as well, although the inclusion of MD on the smallest scale has shown to be beneficial. In the study of Petsev et al. (2015), using the AdResS scheme, SDPD (and therefore SPH) showed that each simulated particle can be seen as a system carrying a set of equations which govern its dynamic behaviour. Changing from the microto a meso-or macroscale description and vice versa does not mean that the particles need to be destroyed or new particles inserted as in a Navier-Stokes-based coupling approach with MD particles. Instead, they may smoothly change from one description to another (changing the equations of motion via blending). This behaviour is truly advantages compared to traditional, grid-based Navier-Stokes algorithms, and further research into particle-based Navier-Stokes methods, such as the SDPD approach, could prove to simplify current challenges arising at the coupling interface.

Summary and conclusions
In the previous five sections, we have looked at different hybrid and multiscale methods that were developed and used in recent years. Two features are commonly the driving motivation for opting to a multiscale or hybrid description over a monoscale approach, to increase the computational efficiency and to enhance details at smaller scales.
We have focused our attention in this review on the most recent and most cited work so as to give an overview of the current state of research. In this spirit, we have not limited ourselves to just one single hybrid/multiscale approach but rather have tried to give a comprehensive overview in various particle-based methods. This also implies that it is not feasible to consider all published literature but rather to give an overview of what each method has achieved to date. It is hoped that the presented studies show an array of different approaches to tackle similar problems and that each method has advantages and disadvantages which may dictate the most suited approach for a particular problem. We have included the lattice Boltzmann and dissipative particle dynamics method in our discussion as two mesoscale representatives. We saw that coupling to other length scales, especially DPD to other particle-based methods, was easily achieved. Molecular dynamics, on the other hand, coupled directly to the Navier-Stokes equations presented difficulties due to disparate degrees of freedom arising from the length and timescale separation. Furthermore, due to the continuum-particle nature in multiscale MD simulations, particles need to be restricted or added/ removed in the hybrid solution interface (HSI) in order to properly model the particle behaviour in the overlap region. This is not a trivial task, and various approaches have been discussed herein. We have discussed possibilities where more than two methods have been combined as seen in the triple-decker approach. MD to Navier-Stokes coupling may benefit from a mesoscopic middleman to reduce the effort to link both descriptions. For example, one could use MD at the atomistic level and couple it with the SDPD method. The SDPD method can then be blended from small to large scales (as done by Kulkarni et al. 2013) and then be coupled with the Navier-Stokes equations on the macroscale. Boundary conditions between MD and SDPD as well as SDPD and Navier-Stokes are easier to impose compared to a direct MD and Navier-Stokes coupling.
At the macroscale, SDPD particles could easily exchange macroscopic quantities such as velocity, pressure and temperature, while on the smaller scales, they could make use of their particle nature to interact with MD particles. This approach presents a simpler coupling between descriptions, and indeed, a trend towards multiscale schemes that make use of more than two methods can be seen. The combination of micro, meso and macroscale together into one single approach may prove to yield more accurate and stable coupling methods, while physical phenomena are adequately modelled across all scales. A different approach was proposed by Asproulis and Drikakis (2013) in which a microscale model was trained by using an artificial neural network (ANN). These are networks aiming to mimic the human's neural network. When presented with training input and knowing the corresponding output, an ANN can be trained to reproduce results for similar data (inside its confidence interval). ANNs are usually used for pattern recognition tasks, for example transforming handwriting into digital text or face recognition. However, Asproulis and Drikakis recognised that in MD to Navier-Stokes coupling, most of the MD simulations on the small scales are presented with similar input from the macroscale (Navier-Stokes solution) and thus produce similar outputs. They trained an ANN with MD data for a Couette flow and showed that it could speed up the computational time by disregarding computationally similar calculations and using the ANN to predict the microscopic properties. The benefit of an ANN is that it can take an arbitrary amount of input data and transform it into an arbitrary amount of output data. Latino et al. (2007) investigated the accuracy of ANNs for various parameters to predict the potential energy surface (PES) by training two different network architectures with MD simulations. They found that for a sufficient amount of training data, the error was less than 1 %. Studies in solid mechanics have also been presented that made successful use of ANNs, see, for example, Hambli (2011) for human bone geometry adaptation over time where long-time intervals prohibit microscale simulations and Unger and Könke (2009) for crack propagation in a reinforced concrete beam.
Returning to multiscale simulations, coupling two descriptions may prove to be a cumbersome task, especially deriving an efficient communication algorithm that works across processors. Tang et al. (2015) recently introduced the multiscale universal interface (MUI), which is a headeronly library that takes most of the programming-related issues off the user's hands. A data sampler, derived from texture sampling concepts, is implemented which allows for spatial data interpolation between descriptions. This is needed in regions where there is a geometric mismatch between the two coupled schemes. The library operates in a multiple program, multiple data (MPMD) approach and uses point-to-point, non-blocking, asynchronous communication protocols for greater flexibility and efficiency. Furthermore, the concept of frames has been introduced which allows the storage of several time steps. This is needed to preserve time information, especially in cases where two solvers iterate at a different pace. Time information can be exchanged at different stages during the execution and the frames can also be used to calculate average quantities. These are sometimes needed where average flow quantities are imposed in the HSI. They applied it to a simple Couette flow and coupled two LAMMPS solver (SPH). Less than 100 lines of code were needed to establish the coupling. The computational efficiency showed good strong scaling over a range of up to 512 cores. In times where opensource fluid dynamic tools are readily available, their MUI library may provide a handy tool to couple different solvers at high efficiency.
In Fig. 1, we have summarised the length scale separations that have been achieved by the different particle-based multiscale methods discussed in the previous sections. Here we plot the smallest characteristic length scales against the size of simulation domain. The smallest characteristic length scales are: molecular diameter (for MD), mean free path (for DSMC), mesh spacing x min (for LBM), particle diameter (for DPD) and average distance between particles (for SPH). These characteristic length scales are chosen arbitrarily but have been found to be a good indicator for the minimum resolution required for each respective method. Since simulations can be done using either dimensional or non-dimensional units, only those studies stating the physical length scales explicitly were considered. In some cases, the smallest characteristic length scales were approximated based on other geometrical information. A few observations can be made by studying Fig. 1. All the molecular dynamics-based multiscale simulations have the same smallest length scale (0.34 × 10 −9 m) since all studies used liquid argon as their working fluid. This is a common choice and fixes the smallest scale, however, the domain size varied in orders of magnitude. A linear trend (in the log-log plot) is evident which shows that all presented studies operated at the same order of length scale separations of approximately O(10 2 )−O(10 3 ). This is also largely true for the non-dimensional simulations which were not considered in Fig. 1. This is explained by the fact that the channel flow is primarily used which in turn indicates that multiscale methods to date are still under development and primarily validated but less frequently applied to complex geometries. The channel flow offers a well understood example for which an analytic formula is available for both the Poiseuille and Couette flow. Hence, the expected scale separation does not change drastically for different particle-based method as the underlying geometry is not changing. However, to make the most use of multiscale simulations that range across scales, new methodologies need to be developed to couple the microscopic with the macroscopic world. It should be mentioned at this point that only those multiscale methods combining two approaches have been considered. We have addressed various studies combining more than two methods, see, for example, Karniadakis (2009), Kacar et al. (2010), Petsev et al. (2015), and it is expected that by combing more than two approaches, further length scale separation can be achieved. In terms of Fig. 1, multiscale approaches ranging over several length scales would be situated in the upper left part of the plot. There is only one example which can be considered a true multiscale approach which bridges the microscale O(10 −9 ) with the macroscale O(10 −1 ). This is the study of Grinberg et al. (2011) who elaborated on the difficulties of parallel computations that occur at these length scale separation. See Sect. 5.1 for a more detailed description. Another observation can be made from Fig. 1: if we define the micro-and macroscale to be approximately of the order of O(10 −9 ) and O(10 0 ), respectively, we can define the mesoscale to sit in-between these two scales. As seen from the plot, all methods, except for molecular dynamics, are located in this area. The DSMC and SPH methods are able to simulated flows on the macroscale as well but also seem to be successful on much smaller scales. In fact, the SPH method ranges from O(10 −9 ) to O(10 2 ) (here we include the SDPD method which makes it possible to be used on such small  Hybrid MD Hybrid DSMC Hybrid LBM Hybrid DPD Hybrid SPH scales), the DSMC method from O(10 −8 ) to O(10 −1 ) and the DPD method from O(10 −9 ) to O(10 −4 ). These methods are valid over a wide range, but their capabilities are not leveraged due to HPC limitations and limited knowledge available on coupling approaches for complex geometries. We have summarised the advantages, disadvantages and applications of each presented particle method in Table 3. It may be consulted together with Fig. 1 to choose a suitable particle method to construct a multiscale scheme.
Finally, we would like to mention the work of Martin Karplus, Michael Levitt and Arieh Warshel who in 2013 received the Noble prize in chemistry for having devised a physical sound manner to communicate information between quantum and classical mechanics in complex chemical systems. Multiscale schemes haven risen over the past decades and have demonstrated that current problems may be successfully simulated with such an approach. They have established themselves, and it is expected that their research interest will continue to grow, helping to demystify current and future challenges in fluid mechanics and beyond.