Mesoscopic simulations at the physics-chemistry-biology interface

Mesoscopic simulations at the physics-chemistry-biology interface

Massimo Bernaschi IAC-CNR, via dei Taurini 19, 00185, Rome, Italy    Simone Melchionna ISC-CNR, c/o Dipartimento di Fisica, Università La Sapienza, P.le A. Moro 5, 00185, Rome, Italy IACS-SEAS, Harvard University, Oxford Street 29, 02138, Cambridge, USA    Sauro Succi Center for Life Nano Science@Sapienza, Italian Institute of Technology, Viale Regina Elena, 295, I-00161 Rome, Italy IAC-CNR, via dei Taurini 19, 00185, Rome, Italy IACS-SEAS, Harvard University, Oxford Street 29, 02138, Cambridge, USA
Abstract

We discuss the Lattice Boltzmann-Particle Dynamics (LBPD) multiscale paradigm for the simulation of complex states of flowing matter at the interface between Physics, Chemistry and Biology.

In particular, we describe current large-scale LBPD simulations of biopolymer translocation across cellular membranes, molecular transport in ion channels and amyloid aggregation in cells. We also provide prospects for future LBPD explorations in the direction of cellular organization, the direct simulation of full biological organelles, all the way to up to physiological scales of potential relevance to future precision-medicine applications, such as the accurate description of homeostatic processes.

It is argued that, with the advent of Exascale computing, the mesoscale physics approach advocated in this paper, may come to age in the next decade and open up new exciting perspectives for physics-based computational medicine.

\stackMath
Contents

I Introduction

Thanks to the spectacular advances of computer technology (hardware and software) on the one side, and mathematical modelling on the other, in the last few decades modern science has come to the point of providing a quantitative description of many biological systems, whose complexity would have been regarded as mission-impossible until only recently. In this Review, we shall illustrate the point through several concrete examples.

Notwithstanding such major advances, the challenge of modelling biological and physiological systems remains formidable, as it actually amounts to cover some ten decades in space (from molecules to the human body) and easily twice as many in time. No mathematical/computational model in the foreseeable future can take up such a challenge head-on, i.e. by direct simulation of all the actual mechanisms, scales and levels involved in the process.

Coarse-grained methods come in many flavours and families, depending on the range of scales and problems they are targeted to, but in this review we shall focus on a specific mesoscale technique, known as the Lattice Boltzmann method (LBM), namely a minimal lattice version of the Boltzmann equation Boltzmann (2012); Cercignani and Berman (1976) which has witnessed a burgeoning growth for the description of complex flow phenomena across an impressively broad range of scales.

The Lattice Boltzmann Equation (LBE) was developed as a computational alternative to the discretization of the Navier-Stokes (NS) equations of continuum fluid mechanics Benzi et al. (1992); Succi (2001).

Over the years, however, it has made proof of an amazing and largely unanticipated versatility and ability to describe a broad variety of phenomena involving complex states of flowing matter, beyond the strict realm of continuum hydrodynamics, including non-trivial flows at micro and nanoscale. Thanks to this versatility, and to the coupling with various families of mesoscale particle methods, in the last decade LBE has gained increased status for the simulation of many complex flow problems at the interface between fluid dynamics, chemistry, material science and biology. These include, for instance, multiphase and multicomponent flows with complex interfaces, the motion of suspended bodies under strong geometrical confinement, possibly with chemical reactions Succi (2018); Krüger et al. (2017). After revisiting the main ideas behind the Lattice Boltzmann (LB) theory, in this review we discuss current and future prospects of multiscale/level Lattice Boltzmann-Particle Dynamics (LBPD) simulations at the physics-chemistry-biology interface, in an attempt to identify and portray outstanding problems of potential relevance to clinical applications in a not-so-distant future, i.e. computational medicine.

Computer-assisted medicine is a consolidated branch of modern science, which generally develops around two separate pillars; molecular biology and macroscale physiology. The former is heavily leaning on bio-informatics and data science tools, while the latter usually relies upon the methods of continuum and fluid mechanics. In this Review, we portray a third, alternative approach, based on mesoscale physics, i.e. physics-informed coarse-grained models of microscale biological processes, possibly embedded within their physiological environment. This mesoscale approach is grounded into the intermediate level of the description of matter, namely kinetic theory, in both its versions, Boltzmann’s kinetic theory and Langevin stochastic particle dynamics. As a result, the main tools of the mesoscale approach are fluids, particles and probability distribution functions. The multiscale LBPD concept is thus apparent; Boltzmann naturally connects upwards to continuum hydrodynamics, while Langevin connects downwards to the molecular level. Combining the two in a single computational harness, opens up a direct route from the continuum to the atomistic word, and back.

Several specific examples have already shown the potential of LBPD simulations in areas straddling across Physics and Biology. A selected set of applications will be described in detail in this Review, to provide a taste for the breadth of applications that have been tackled in the recent past. Readers keen on specific details are kindly directed to the original literature.

Clearly, the PCB interface is enormously rich and varied, ranging from nanometric macromolecular phenomena, to peptidic aggregation and biopolymer translocation, to cellular motion and active matter, as sketched in Fig. 1. In a few words, the topic goes way beyond the scope of the present work. Nevertheless, it is important for the reader to appreciate the great flexibility of the method to cope with multiple scales of motion and, more importantly, its ability to incorporate the desired degree of bio-chemical specificity inherent to the problem at hand.

Figure 1: Scales related to biological and medical applications, whose stretch can be covered by the LBPD approach. While the particle description enables the representation of complex macromolecules or cellular organizations by tracking the fate of each individual particle, the LB description is based on the collective motion of solvent molecules. The boundary between LB and PD is blurred and depends on the degree of microscopic detail required by each single application (see also Fig. 8).

The LBPD framework can be enriched in the direction of describing complex flowing states of matter at micro and nanoscales. Although not “rigorous”, such variants prove capable of providing new physical insights into highly complex states of flowing matter. Remarkably, such extensions can be put in place without compromising the outstanding parallel amenability to parallel processing of the method. On the assumption that computing power will keep growing in the next few decades, if only perhaps at sub-Moore’s paces, this state of affairs spawns tremendous opportunities to gain new insights into a series of fundamental problems dealing with complex states of flowing matter in general, and with special focus on those relevant to biology and medicine.

The paper is organised in three main parts. In the first one, we discuss the basic aspects of Boltzmann’s kinetic theory with special emphasis on its lattice version for fluid dynamics and its extensions to soft-matter and biological applications, including the coupling to particle dynamics for the motion of suspended bodies.

In the second part, we describe selected applications to biological and physiological systems, such as biopolymer translocation, ion channels, protein diffusion and amyloid aggregation in cellular environments. For a quick visual appreciation, see figures 2, 3, 4, to be commented in detail later in this manuscript.

Figure 2: Translocation of a biopolymer, as in the case of DNA, through a narrow pore. Different representations can be used to study multiple levels of detail during the translocation process. a,b) the biopolymer is represented as a simple necklace of beads, by neglecting correlations stemming from the local molecular rigidity or backbone charge. c,d) at the next level, the macromolecule is charged and moves in an electrolyte solution, whereby a neutral solvent and counterions and coions migrate due to an externally applied electric field, giving rise to an electroosmotic flow that ultimately causes the molecule to translocate Datar et al. (2017).
Figure 3: Representation of the crowded interior of the cell as obtained from simulations Bernaschi et al. (2013a). The inset illustrates the embedding of a protein in the LB mesh and each protein atom is represented via the DPM particle-fluid exchange scheme.
Figure 4: Left: snapshot of a configuration of amyloid peptides A simulated in a cubic box of size nmand complex flow structure generated by their motion in the surrounding solvent. Right: evolution of the size of the peptidic aggregates as a function of time Nasica-Labouze et al. (2015); Chiricotto et al. (2017).

In the third part, we provide a prospective view of a series of problems at the physics-chemistry-biology interface, which may become accessible once Exascale computers are with us. Special attention is paid on their potential import for clinical applications, such as the direct simulation of biological organelles and the quantitative description of hemostatic processes.

Figure 5: Snapshot of molecules and particles in a multiphase flow , either a) dragging the fluid from one phase to another, b) sitting at the fluid-fluid interface (Sega et al. (2013) – Reproduced by permission of The Royal Society of Chemistry), in absence c) and in presence of colloids d) (Tiribocchi et al. (2019) – Reproduced by permission of The Royal Society of Chemistry) modifying the shape of the interface to a corrugated one to reflect molecular correlations.

Finally, due to the crucial role played by high-performance computing in this story, in the Appendix we provide an extended account of the main issues involved with the implementation of the LBPD scheme on high-end parallel computers in the Exascale range.

The main message we wish to convey in this Review is that a mesoscale physics-based approach to computational medicine may come to age in the next decade.

Ii Boltzmann Kinetic Theory

The Boltzmann equation (BE) is the core of Boltzmann’s kinetic theory, that, in turn, is the cornerstone of non-equilibrium statistical mechanics, a pillar of theoretical physics at large Boltzmann (2012). Besides its paramount conceptual value as a bridge between the microscopic world of atoms and molecules and the macroscopic world of thermo-hydrodynamic fields, the BE also provides a concrete tool for the quantitative investigation of a broad variety of practical non-equilibrium transport problems Cercignani and Berman (1976). However, the BE is all but an easy piece to work with: a non-linear integro-differential equation in (phase-space plus time) dimensions. This motivates a relentless search for new methods to solve the BE either analytically or numerically, the latter option usually covering a broader ground. Graeme Bird’s Direct Simulation Monte Carlo (DSMC) method has played a leading role in this respect and continues to do so to the present days Bird (1994). In principle, DSMC solves the BE directly and in full, i.e. accounting for the specificity of molecular interactions, as well as strong non-equilibrium effects, using a stochastic particle technique, whence the Monte Carlo label. This comes at a major computational cost, which is why various approximations have been developed and considerably refined over the years Ketsdever and Struchtrup (2016); Dimarco et al. (2018). Close to local equilibrium and away from confining elements, however, molecular details become increasingly irrelevant. Universality takes stage and more economical descriptions can be devised. The basic idea is to relinquish the “irrelevant” details while still preserving the basic properties of macroscopic physics, namely the symmetries and conservation laws which secure the emergence of the NS equations from the underlying molecular dynamics. Among others, a description which has gained major interest for the last three decades is the LB method Benzi et al. (1992); Succi (2001, 2018); Krüger et al. (2017). LB was devised with the specific intent of providing an alternative to the discretization of the NS equations for the numerical solution of continuum hydrodynamic problems. This still is its mainstay and, for some authors, also the only place where it belongs. Indeed, the use of LB for flows beyond NS was ruled out Luo (2004); Junk et al. (2005), mostly on account of the lack of a rigorous asymptotic limit. The above no-go has been proven largely over-restrictive and nowadays applications beyond the strict realm of continuum fluid dynamics abound, especially in the direction of soft matter. Since problems in biology and medicine hardly involve fluid mechanics alone, such developments are of direct relevance to computational explorations at the interface between physics, chemistry and biology, the main scope of this Review.

ii.1 The Boltzmann Equation

The Boltzmann equation (BE) of classic kinetic theory is basically a continuity equation in six-dimensional phase-space, namely Boltzmann (2012):

(1)

where is the probability density of finding a molecule at position in ordinary space, with velocity at time . The left-hand side of the BE represents the streaming of the molecules under the effect of a force field and reflects the Newtonian mechanics and of a representative molecule, while the right-hand side describes intermolecular interactions. For the case of a dilute gas, as originally considered by Boltzmann, these interactions are typically taken in the form of two-body local collisions, since higher order encounters are much less frequent, hence negligible in the collision count. Even under such major simplifications, solving the BE presents a very daunting challenge on account of its high-dimensionality, six phase-space dimensions plus time, as well as due to the non-linear (quadratic) integral character of the collision operator Cercignani and Berman (1976).

Regardless of the complexity of the underlying microscopic interactions, the collision operator must comply with the mass-momentum-energy conservation laws, namely:

(2)

In addition, it must also secure compliance with the Second Principle, which amounts to supporting a so-called H-theorem, namely:

(3)

In other words, the dynamics of the distribution function must converge to a universal global attractor, corresponding to the thermodynamic equilibrium.

The macroscopic fluid variables are obtained by a linear and local contraction of the Boltzmann distribution, namely:

(4)
(5)
(6)

where is the mass density, the flow speed, and the flow temperature in spatial dimensions.

Central to the emergence of hydrodynamic behavior is the notion of local equilibrium. This is defined as the specific form attained by the Boltzmann distribution once collisions come in complete balance, i.e.

(7)

Inspection of the Boltzmann collision operator provides the following universal Maxwell-Boltzmann (MB) local equilibrium distribution:

(8)

where is a normalization constant in spatial dimensions, is the thermal speed and

is the peculiar speed, i.e., the molecular velocity relative to the fluid one, in units of the thermal speed. The reader familiar with statistical mechanics will readily recognize the canonical distribution in the co-moving frame of the fluid, with the identification .

A few comments are in order.

First, the MB distribution depends on space and time only through the hydrodynamic fields, , its dependence on the velocity variable being a universal Gaussian distribution. This is a strict consequence of Eq. 2, i.e., the microscopic conservation laws.

Such dependence is largely arbitrary, with the caveat that it should be weak on the molecular scale. More precisely, the macrofields should not show appreciable changes on the scale of the molecular mean-free path, that is

(9)

where designates any macrofield and is the molecular mean free path. The above ratio, known as Knudsen number, serves as the smallness parameter controlling the emergence of the hydrodynamic limit from Boltzmann’s kinetic equation. Ordinary fluids dynamics holds in the range and below.

ii.2 From Boltzmann to Navier-Stokes hydrodynamics

The conceptual path from BE to the NS equations of continuum fluids is based on two fundamental steps:

  1. Projection of the Boltzmann equation upon a suitable basis function in velocity space, typically Hermite polynomials in cartesian coordinates,

  2. Multiscale expansion using the Knudsen number as a smallness parameter, on the assumption of weak departure from local equilibrum.

The projection generates a hierarchy of partial differential equations for the kinetic moments:

(10)
(11)
(12)

where

(13)

and denotes the -th order tensor Hermite polynomial. Note that is a tensor of rank , namely (scalar) is the fluid density, (vector) is the fluid current and (second order tensor) is the momentum-flux, whose trace delivers (twice) the kinetic energy of the fluid and the triple tensor is the flux of momentum flux.

The left hand side clearly shows that the moment hierarchy is open, since the time derivative of is driven by the divergence of . This is the mechanism by which heterogeneities fuel non-equilibrium. Also to be noted that the right-hand-side of the first two equations is zero because collisions conserve mass and momentum. However, they do not conserve momentum-flux, which is why the right-hand-side of the third equation is non-zero, expressing the relaxation of the momentum-flux to its equilibrium expression. Such relaxation takes place on a collisional timescale , which in turn fixes the kinematic viscosity of the fluid. Instantaneous relaxation () denotes the infinitely strong collisional regime whereby collisions do not leave any chance to non-equilibrium to survive, formally corresponding to the idealized case of a perfect (zero dissipation) fluid.

Macroscopically, this corresponds to the inviscid Euler equations.

All real fluids, though, relax in a short but finite time (strictly speaking this is also true for superfluids), and consequently the moment equations present an open hierarchy which needs to be closed, somehow. This is where the assumption of weak departure from local equilibrium takes stage.

By that assumption, one formally expands the distribution function and space-time derivatives in powers of the Knudsen number, replaces the expansion in the moment equations, and collects homologue terms order by order in the Knudsen number. To zero order, the Euler equations are obtained, whereas the first order delivers the NS equations of dissipative fluids, namely

(14)
(15)

where is a second order tensor formed by spatial derivatives of the flow field. Typically:

(16)

where is the symmetrized gradient tensor, is the divergence of the flow and is the unit tensor. The first scalar is the dynamic shear viscosity whereas associates to the bulk viscosity, .

Despite their deceivingly simple physical content, essentially mass and momentum conservation (Newton’s law), as applied to a finite volume of fluid, the NS equations prove exceedingly difficult to solve, as they involve the non-linear evolution of a three-dimensional vector field, often in a complex geometry set up. This sets a formidable challenge to even the most advanced computational methods, whence the ceaseless hunt for more efficient numerical methods. Those methods include a broad array of techniques to discretize the NS equations, using finite differences/elements/volume schemes.

Three decades ago, however, an entirely different route was devised, which consists in attacking fluid dynamics “from the bottom”, i.e., appealing to a microscopic description of the fluid states matter, namely a highly stylized version of molecular dynamics known as Lattice Gas Cellular Automata Chopard and Droz (1998); Rivet and Boon (2005). In a nutshell, the idea is to introduce a boolean lattice fluid, consisting of a set of boolean particles whose dynamics is confined to the lattice sites. Boolean, here, means that the state of the system at a given lattice site and instant in time is uniquely defined by a set of binary digits, coding for the absence/presence of a corresponding particle moving with unit speed along one of the links connecting each lattice site to its neighbours. By a suitable choice of the lattice connectivity and interaction rules, such Boolean system can be shown to reproduce the NS equations of continuum fluid dynamics. An associated LB equation was also derived in the process of taking the boolean automaton to NS, but its computational capabilities went under noticed. Even though the LGCA did not make it into a competitive tool for computational fluid dynamics, it nonetheless paved the way to the idea of computing fluid flows by simulating fictitious particle dynamics instead of discretising the NS equations. LB Higuera et al. (1989) fully inscribes within this line of thought.

ii.3 Boltzmann equation for Biology?

The Boltzmann factor is a household name in biology, as it governs the statistical behaviour of a broad class of equilibrium and non-equilibrium (activated processes) phenomena of utmost relevance to biological systems. But, how about the Boltzmann equation? At first sight, Boltzmann kinetic theory, in its original form at least, offers little scope for biological applications, since it formally applies to dilute states of matter in which molecular collisions are rare, the so-called weakly coupled regime, in which many-body interactions can safely be neglected with respect to binary molecular encounters. On the contrary, most biological phenomena are hosted primarily by condensed and soft matter systems in which many-body effects play a primary role.

Here, however, a change of perspective proves exceedingly fruitful. Rather than the original Boltzmann equation for actual molecules, we shall refer to model Boltzmann’s equations for fluid-particles, the latter denoting the effective degrees of freedom describing the behaviour of representative groups of molecules. Historically, the distinctive feature of model Boltzmann equations is a dramatic simplification of the collision operator, in an attempt to relinquish most mathematical complexities while retaining the essential physics at hand. The most popular Boltzmann model is the celebrated Bhatnagar-Gross-Krook equation, in which the collision operator is replaced by a simple single-time relaxation term Bhatnagar et al. (1954):

(17)

where is the local equilibrium and is a relaxation frequency controlling the relaxation to the local equilibrium on a timescale .

The key advantage of this change is neat: model equations are more flexible by construction, hence they can be modified (extended) not only to simplify the collision operator but also to describe a wide host of physical effects not included in the original Boltzmann equation.

In particular, they can reinstate the effects of, many-body interactions, via effective one-body forces, in the spirit of Density Functional Theory (Vlasov-Boltzmann Equation) Hansen and McDonald (1990) statistical fluctuations, through appropriate stochastic sources (Fluctuating Boltzmann Equation) Ladd (1993) far-from-equilibrium inhomogeneities, via suitably extended collision-relaxation operators in which the relaxation time is promoted to the status of a self-consistent dynamic field Chen et al. (2003); Higuera et al. (1989); d’Humières (1992).

The three extensions above inscribe to the general framework of Reverse Kinetic Theory (RKT), the strategy whereby the kinetic equation is designed top-down, based on prescriptions securing compliance with macroscopic hydrodynamics in the first place, and then adding “molecular” details “on-demand” by the specific application under investigation.

RKT reverses the canonical bottom-up route, whereby kinetic equations are derived from the underlying microscopic models and proves quite effective in bringing Boltzmann-like equations within the realm of condensed and soft matter physics. However, care must be exercised in securing compliance of the top-down approach with the basic principles of statistical physics, namely symmetries/conservation laws as well as evolutionary constraints (the Second Principle and its local form, known as H-theorem).

This all said and done, a practical question still remains: although simplified, the model BE’s still leave in an unwieldy dimensional space. Here, a time-honoured ally of any computational scientist, the lattice, makes its glorious entry.

ii.4 Hydrodynamics for Biology

Hydrodynamics and biology couple across an amazingly broad spectrum of scales, ranging from the macromolecular level in the compartments of living cells, all the way up to industrial bioreactors, rivers, lakes and oceans. From macromolecular motion in crowded cellular spaces, to the deformations of membranes, to the motion of cells and the self-propulsion of bacteria, the common goal is to capture the effect of hydrodynamics under both equilibrium and non-equilibrium conditions. The latter may arise under the effect of an external flow, such as the transport of proteins and cells in the blood stream, or via the consumption of energy in cellular metabolic pathways. Indeed, the response of biological matter to temperature, pH or mechanical forces, plays a key role in most biological processes. For instance, biological macromolecules, cells or tissues, are usually fragile and easily damaged by hydrodynamic or shear forces, as they occur far from equilibrium and under strong confinement. Following in the footsteps of hydrodynamic experimental techniques, customarily used to measure basic molecular properties, such as weight, size and shape, computer simulations can be used in biophysical chemistry to explore and assess the mechanical and dynamical properties of macromolecules far from equilibrium, hence much closer to the in vivo conditions relevant to medical purposes. The benefits of including hydrodynamic forces is even more apparent whenever thermodynamic forces, stemming from solute-solvent interactions, can also be taken into account. A comprehensive computational framework including both hydrodynamic and thermodynamic forces, stands therefore as a most desirable target. Owing to its inherently intermediate nature, between atomistic and continuum descriptions of matter, (lattice) kinetic theory sits at a vantage point to meet this goal.

Iii Lattice Boltzmann for Continuum Hydrodynamics

The basic idea of the LBM is to constrain the velocity degrees-of-freedom to a discrete lattice with sufficient symmetry to protect the conservation laws which secure the emergence of standard hydrodynamic behavior in the macroscopic limit. LB is based on the idea of representing fluid populations on a uniform, cartesian grid.

The standard LB scheme in single-relaxation time (BGK) form reads as follows Higuera et al. (1989); Qian et al. (1992):

(18)

where is the discrete Boltzmann distribution associated with the discrete velocity , running over the discrete lattice, to be detailed shortly. In the above, is a relaxation parameter controlling the fluid viscosity and is the lattice local equilibrium, basically the local Maxwell-Boltzmann distribution truncated to the second order in the Mach number. The truncation is not a luxury, but a requirement dictated by Galilean invariance (GI).

Let us remind that GI refers to the invariance of the NS equations under an arbitrary change of the local fluid velocity ; upon such change, the NS equations stay the same, provided is replaced by . In continuum kinetic theory, Galilean invariance is encoded within the dependence of the local Maxwell-Boltzmann distribution on the relative velocity of the molecules with respect to the fluid one, namely . As a result, an observer in the comoving frame, that is a frame moving at the local fluid velocity, experiences the same local equilibrium as if there were no fluid motion. To be noted that, in the above, is an arbitrary function of space and time, indicating that Galilean invariance is a local and continuum symmetry, i.e., it holds even if different regions of the fluid move at different velocities, which is of course the case in most fluids of practical interest. Galilean invariance is reflected by the specific form taken by the moments of the equilibrium distribution, and specifically by those explicitly relevant to hydrodynamics, namely mass, momentum and the momentum flux-tensor, that is:

where latin subscripts run over spatial coordinates . The above expression follows straight from the property of gaussian integrals in velocity space.

One might naively expect that the same would be true in the lattice, provided the local Maxwell-Boltzmann distribution is retained with the plain replacement . Straightforward algebra shows that the situation is different not for a mere mathematical accident, but as a consequence of the fact that a local and continuum symmetry cannot remain unbroken in a discrete lattice. More precisely, it cannot remain unbroken for any arbitrary velocity field. It turns out, though, that this is possible whenever the fluid velocity is much smaller than the sound speed, i.e., in the low-Mach number limit.

Under such a perturbative approximation, replacing velocity integrals with discrete sums returns exactly the same moments as in the continuum, namely:

(19)

provided that lattice tensors up to fourth order are isotropic. Of course, this by no means implies that GI is fully restored, but simply that the GI-breaking terms are confined to kinetic moments higher than order four.

In other words, the hydrodynamic constraints can be matched perturbatively, by expanding the local Maxwellian to second order in the Mach number, which is sufficient to recover the (isothermal) NS equations, since those equations are quadratic in the fluid velocity. Full Galilean invariance for an arbitrary flow field implies instead an infinite series in the Mach number, corresponding to the full expansion of the local Maxwellian in Hermite polynomials.

The actual expression of the discrete local equilibria reads as follows:

(20)

where and represent the dipole and quadrupole contributions, respectively. Hereafter, repeated indices are summed upon. In the above, is a set of lattice-dependent weights, normalized to unity, which represent the lattice analogue of the global (no-flow) Maxwell-Boltzmann distribution. Finally is the lattice sound speed.

We hasten to note that, at variance with its continuum counterpart, the expression (20), being a polynomial truncation of the Maxwell-Boltzmann equilibria, is non-negative definite only in a finite range of fluid velocities, typically of the order of . This configures the LB method as an appropriate description of quasi-incompressible, low Mach-number flows.

The discrete velocities matching isotropy constraints up to fourth order can be shown to be subsets of the D3Q27 (DQ is a widely used notation to indicate a LB scheme in dimensions using a set of velocities) mother lattice with discrete speeds in three spatial dimensions. D2Q9 and D3Q27 are the direct tensor product of the elementary one-dimensional D1Q3 stencil, , in two and three spatial dimensions, respectively (see Fig. 6).

Figure 6: Examples of standard 2d and 3d LB lattices, with , and discrete speeds, typically denoted as D2Q9, D3Q19 and D3Q27, respectively.

The expansion (20) implies that hydrodynamic LB flows are bound to be quasi-incompressible, i.e., Mach number well below unity. Likewise, third order kinetic moments, describing energy and heat flux, are not correctly reproduced, since these terms require sixth-order isotropic lattice tensors. Such constraints can not be met by lattices confined to the first Brillouin region described by D3Q27: higher order lattices extending beyond the first Brillouin cell are necessary. So much for the low Mach number approximation, which is specific to the lattice. But what about the low Knudsen limit, which, on the contrary, lies at the very roots of the convergence of Boltzmann to NS?

Since the Knudsen number controls the heterogeneity-driven departures from local equilibrium, it is intuitively clear that the low-Knudsen hydrodynamic limit implies further constraints on the non-equilibrium component of the momentum-flux tensor, which amounts to recovering the continuum expression of the stress tensor.

Based on (19) and (20), it can be shown that in order for such constraints to be met, isotropic tensors, again up to order four, need to be exactly reproduced in the lattice. It may come as a surprise that non-equilibrium constraints can be matched at the same order of isotropy of the equilibrium ones. The reason is parity-invariance, namely the fact that each discrete lattice velocity comes with a mirror partner . Indeed, non equilibrium constraints involve lattice corrections driven by fifth-order lattice tensors, which are identically zero due to parity invariance. Readers interested in the straightforward but laborious algebraic details, may look up the vaste literature on the subjectChopard and Droz (1998); Rivet and Boon (2005); Hénon (1987).

Next, let us comment on the source term, , at the right-hand side of the LB equation (18). This term stands for a generic source of mass/momentum/energy, describing the coupling of the fluid to the surrounding environment. Mass sources are typically associated with the presence of chemical reactions, turning one species into another in multi-component versions of the LB for reactive flows.

In the case of inert flows, the source term typically encodes the momentum exchange due to the coupling to external (or internal) fields, such as gravity or more complex interactions, like self-consistent forces reflecting potential energy interactions within the fluid, as well as thermal fluctuations. As we shall see, the two latter cases are crucial for the extension of LB beyond NS hydrodynamics.

From the operational standpoint, the source term acts like a bias promoting the populations which move along the local force field and setting a penalty on those that move against it. It is therefore clear that the strength of such term is subject to stringent stability and positivity constraints.

iii.1 From Lattice Boltzmann to continuum hydrodynamics

The conceptual path taking LB to NS is exactly the same as in the continuum theory, with the crucial caveat of turning around the (many) catches associated with lattice discreteness. To make a long story short, it amounts to securing the proper symmetries of the lattice tensors entering the set of (lattice) moment equations. Like always with lattice physics, the name of the game is to erase the lattice dependence to the desired order. For the case of isothermal, incompressible fluids, the desired order is the fourth one. In equations, and using coordinate notation for tensors Hénon (1987):

(21)
(22)
(23)

where latin indices run over spatial coordinates. In the above, is a set of weights normalized to unity and odd-rank tensors are automatically equal to zero by parity invariance, i.e., every discrete velocity comes with an equal and opposite partner, so that the sum of the two gives a null. Finally is a lattice dependent constant expressing the sound speed in the lattice, as per the expressions above.

With the symmetries secured, everything proceeds like in the continuum theory, with another important caveat though, namely the fact that the effective mean-free path of the lattice fluid is replaced by the lattice spacing whenever the latter is larger than the physical one, the typical case in most macroscopic LB simulations.

This is a very technical and somehow thorny issue, whose details are out-of-the-scope of the present Review. Nevertheless, we wish to caution the reader that the physics taking place at the scale of a few lattice spacings should always be inspected with a big pinch of salt, because it is constantly in danger of breaking the hydrodynamic assumptions.

Once these catches are disposed of, one ends up with a lattice fluid obeying an ideal equation of state and a kinematic viscosity given by:

(24)

where is the natural lattice viscosity. Note that the stability range of the discrete time-marching, , also secures the positivity of the kinematic viscosity.

Both ends of this range must be handled with care. In the low-viscosity regime , typical of turbulence, strong gradients may develop posing a serious threat to the numerical stability of the scheme. A powerful variant of the basic LB includes a self-consistent tuning of the relaxation parameter so as to ensure compliance with local entropy growth (H-theorem) Karlin et al. (1999); Succi et al. (2002). That variant, known as Entropic Lattice Boltzmann (ELB) is normally intended to simulate high-Reynolds macroscopic flows, but lately is proving very effective also to stabilize microscale LB simulations with sharp interfaces.

In the opposite high-viscous regime, , the one most relevant to biological applications, the viscosity may formally diverge, signalling a departure from hydrodynamic behaviour, due to the onset of ballistic motion. As a result, whenever possible, LB simulations should be kept away from both limits, say .

We shall return to this important point in the Section ”Future Challenges”.

iii.2 Boundary conditions

In the early days, boundary conditions were hailed as one of the main assets of LB, and, to a certain extent, they still are. As a matter of fact, since information always travels along straight lines, even complex geometries can be handled by comparatively straightforward computational methods based on elementary mechanical operations. For instance, no-slip on solid walls can be implemented through a simple bounce-back between distributions propagating along opposite directions, i.e., from the fluid to the wall and viceversa (See Fig. 7).

Figure 7: The effect of wall boundaries on the LB populations according to the Bounce Back scheme. Populations hitting a wall node are reflected from the incoming fluid direction, giving rise to no-slip flow conditions.

For the sake of generality, let us consider a fluid flowing in a bounded domain confined by a surrounding boundary . The problem of formulating boundary conditions within the LBE formalism consists in finding an appropriate relation expressing the incoming (unknown) populations as a function of the outgoing (known) ones . Outgoing populations at a boundary site are defined by the condition

where is the outward normal to the boundary cell centered in . Incoming populations are defined by the opposite sign of the inequality. In mathematical terms, this relationship translates into a linear integral equation

(25)

where the kernel of the boundary operator generally extends over a finite range of values inside the fluid domain. This boundary operator reflects the interaction among the fluid molecules and the molecules in the solid wall. Consistently with this molecular picture, boundary conditions can be viewed as a special (sometimes even simplified) type of collisions between fluid and solid molecules. Physical fidelity can make the boundary kernel quite complicated, which is generally not the idea with LBE. Instead, one usually looks for expressions minimizing the mathematical burden without compromising the essential physics.

In particular, one seeks minimal kernels fulfilling the desired constraints on the macroscopic variables (density, speed, temperature and possibly the associated fluxes as well) at the boundary sites . This may lead to a mathematically under-determined problem, more unknowns than constraints, opening up an appealing opportunity to accommodate more interface physics within the formulation of the boundary conditions. However, it also calls for some caution to guard against mathematical ill-posedness. By now, a vast literature which, however, goes beyond the scope of the present Review, is available on this topic Krüger et al. (2017); Zou and He (1997); Chen et al. (1996).

Here, we just wish to point out that the main issue controlling the complexity of the boundary problem is whether the collection of boundary points lies on a surface aligned with the grid (so-called staircase wall boundary) or it is given by a set of off-grid surface elements (sometimes called surfels). The latter case is significantly more complex and requires extra-care, the immediate benefit being that the near-wall physics is second-order accurate as compared to the first-order accuracy of the staircased approximation (for a detailed account see Succi (2018); Krüger et al. (2017)). Advanced applications of the off-lattice boundary method prove capable of dealing with fairly complex geometries, an important but highly technical topic which also goes beyond the scope of the present Review. For full details see the recent books Krüger et al. (2017); Succi (2018)

iii.3 The bright sides of Lattice Boltzmann

Why does LB represent an appealing scheme for simulating complex states of flowing matter ?

Several features merit highlight, but essentially they all track down to the benefits of working with extra-dimensions opened up by the six-dimensional phase-space inhabited by kinetic theory.

More in detail, the upshots are the following.

First, the information always travels along straight lines, regardless of the space-time complexity of the emergent hydrodynamics.

What we mean by this is that the discrete distributions move in sync along the discrete light-cones, defined by:

no matter the complexity of their space-time dependence.

Since the discrete velocities do not depend on space and time, the streaming is exact, no information is lost in hopping from one lattice site to another, literally an error-free operation. This stands in sharp contrast with the self-advection term of macroscopic hydrodynamics, whereby the velocity fields, or any other field for that matter, moves along space-time dependent material lines defined by the fluid velocity itself, typically a complicated function of space and time.

Such independence explains the outstanding LB amenability to parallel computing, a practical asset which can not be underestimated.

Second, space and time always come on the same (first-order) footing. In particular, this means that dissipation is not expressed by second-order spatial derivatives (Laplace operator) but simply as a local relaxation to a local equilibrium, as described earlier. This is a substantial advantage for confined flows, whose dynamics is largely dictated by the presence of solid boundaries, where spatial gradients are usually at their largest. Taking space derivatives near the boundaries is notoriously prone to numerical inaccuracies, especially whenever the geometry of the domain is irregular, as it is often the case for biological flows. LB equipped with suitable formulations for curved boundaries can significantly mitigate such difficulties.

Third, LB carries pressure as just any other macroscopic field, with no need of solving a (usually very expensive) Poisson problem to compute the pressure field consistent with an incompressible flow. This is because the discrete distributions carry the momentum-flux tensor “on their back” and consequently the pressure obeys its own dynamic equation. More specifically, the second order momentum-flux tensor, , obeys the relaxation equation

(26)

where is the third-order energy flux tensor, latin indices running over spatial dimensions. In the limit where the tensor is enslaved to its equilibrium value, the time derivative can be dropped and the energy-flux tensor can be approximated by its equilibrium expression. Under such conditions, the above equation delivers .

Once the due lattice symmetries are fulfilled, the equilibrium component delivers the advection and pressure terms of the NS equation, while the third order non-equilibrium term delivers the dissipative term.

Fourth, coupling of the fluid to a broad variety of other physical phenomena, is readily achieved by formulating suitable expressions of the source term .

Typically this reflects the action of mesoscale forces describing the effective interactions between fluid-fluid and fluid-solid molecules. The positive side effect is that an entire new world of complexity may be supported simply by inserting a few tens of lines of additional code.

This feature is crucial to the extension of LB beyond the traditional realm of dilute gases and particularly to biological flows, as we shall detail in the sequel.

Once again, this rosy picture only refers to the conceptual fresco, actual implementations requiring great care in dribbling lattice artefacts.

Figure 8: The description of Physical-Chemical-Biological (PCB) systems requires a suitable integration of microscopic details (Individuality) within the universal harness dictated by symmetry principles and ensuing conservation laws which govern the macroscopic behaviour (Universality). Kinetic theory is expected to strike an optimal (problem-dependent) balance between Universality and Individuality.

iii.4 The dark side of the LB moon

Following a witty line by Daan Frenkel, it is often more instructive to analyse what can go wrong in a computer simulation, rather than indulging in the glorious description of its success stories. To abide by this wise attitude, in the following we mention things that can still go wrong with LB simulations, the dark side of the LB moon.

LB draws much of its computational simplicity to the features of uniform and regular lattices, that, at times, recalls an ideal Legoland. Realistically complex geometries may sometimes challenge this setup, and call for more flexible and elaborate formulations, such as interpolation procedures for curved boundaries, local grid-refinement, mergers with finite-volumes techniques, to name just a few of some popular options. It is only fair to acknowledge that implementing such extensions maybe laborious and possibly tax the computational simplicity as well.

Another limitation is the finite compressibility. While buying major computational savings through dispensing with the Poisson problem for pressure, a finite amount of compressibility must be tolerated in return. Such effects are usually negligible for low-Reynolds flows, but they need nevertheless to be watched carefully in open flows, especially at outlets, where spurious density waves may eventually back-propagate into the fluid and undermine the accuracy and sometimes even the stability of the simulation.

A similar observation applies to flows with strong thermal transport. For one, it should be appreciated that LB is essentially athermal, in the sense that the discrete Boltzmann distribution is represented by a linear superposition of mono-chromatic beams with no dispersion in velocity space, hence, zero temperature in a strict equilibrium thermodynamic sense. Nevertheless, a kinetic temperature can still be defined as a measure of the dispersion around the mean flow velocity, namely, in spatial dimensions, . The inclusion of strong thermal effects typically requires higher order lattices, accommodating sixth order tensors describing the flux of energy Krüger et al. (2017); Nie et al. (2008). Several alternatives are also possible which go beyond the scope of the present work (for details see chapter 22 in Succi (2018)). Despite major progress, it is fair to say that thermal LB schemes still lag behind their athermal counterparts in terms of numerical robustness. Consequently, LB is often coupled to independent thermal solvers typically based on finite-difference or finite-volume methods.

These limitations off the breast, we believe it is fair to surmise that the appealing features of the LB method largely overweigh its weaknesses.

Iv Lattice Boltzmann for Generalized Hydrodynamics

The ideas and methods presented so far deal with flows of “simple” fluids that can be described in terms of the NS equations (note that simple fluids can give rise to highly complex flows, turbulence being a prime example in point).

Modern high-tech applications, nano-engineering and biology in the first place, set a pressing demand of quantitative understanding of more complex states of flowing matter, in which fluids interact with external or self-consistent fields, undergo chemical reactions, phase-transitions, or interact with a variety of suspended bodies, such as colloids or biological molecules.

This emerging sector of modern science, often referred to as “complex fluids”, or more trendily, “soft matter physics”, portrays a multidisciplinary scenario whereby fluid dynamics makes contact with other disciplines, primarily chemistry, material sciences and biology as well.

There is a growing evidence that LB, and extensions thereof, holds a vantage point as a computational framework for the simulation of complex states of flowing matter. Ideally, LB would fill the gap between fluid dynamics and molecular dynamics, namely the huge and all-important region where fluid dynamics breaks down and molecular dynamics is not yet ready to take over for lack of efficient algorithms and computing power (see Fig. 8) Boon and Yip (1991).

LB is a natural candidate to fill that gap because of its mesoscopic ability to incorporate microscopic details into the kinetic theory formalism via suitable external fields and/or equivalent generalizations of local hydrodynamic equilibria: a microscope for fluid mechanics, a telescope for molecular dynamics Succi (2018).

The extension of LB to generalised flows is based on a number of major upgrades of the basic LB theory. In the sequel, we focus on the following selection:

  • Reactive systems

  • Charged flows

  • Flows far from equilibrium

  • Fluctuating hydrodynamics

  • Non-ideal fluids

  • Flows with suspended bodies

We now proceed to a more detailed discussion of the items above.

iv.1 Advection-Diffusion-Reaction systems

Chemical reactivity is an essential element to deal with when facing the task of simulating systems at the PCB interface Boon et al. (1996); Succi (2018); Coveney et al. (2016a); Alowayyed et al. (2017).

Biochemical reactions lie at the heart of most biological phenomena, controlling species interconversion. They involve the breaking and making of covalent bonds, often catalyzed by enzymes, ubiquitous in all metabolic and synthetic reactions within the cell and in the body. Reactions occur in bulk conditions in single-phase (homogeneous) or multi-phase (heterogeneous) environment, typically at the interface between phases or in a porous-like environment. Biological reactions take place under a wide host of different conditions. The antigen-antibody binding in solution or the binding of small molecules to plasma proteins in blood, clotting reactions on the surface of blood vessels, the hydrolysis of adenosine triphosphate (ATP\nomenclature[ADR]ATPadenosine triphosphate) on the surface of endothelium, ion transport, and the release of nitric oxide are examples of homogeneous reactions. Heterogeneous reactions within tissues include aerobic metabolism, glucose consumption, and receptor-mediated endocytosis.

In all cases, reactions can occur either in no-flow or under flow conditions, as for example the enzyme reactions on the surface of the blood vessels. Diffusion and convection affect the rates of homogeneous and heterogeneous reactions. Because reactants must diffuse or convect to the surface where they react, the mass transfer mechanisms occur in sequence with the reaction process. In many cases, the effects of reaction and diffusion on the reaction rates cannot be easily separated and numerical methods provide the only viable route to study the process.

LB is a powerful framework to capture diffusion and reaction since its mathematical apparatus is by no means confined to fluid equations.

As a matter of fact, by suitable tuning of the local equilibria, relaxation matrix and external forces, it permits to generate a very broad variety of linear and non-linear partial differential equations, including those describing relativistic and non-relativistic quantum phenomena Succi (2002a); Mendoza et al. (2010).

Of relevance are advection-diffusion-reaction (ADR) \nomenclature[ADR]ADRAdvection-Diffusion-Reaction equations of the general form:

(27)

where is a species concentration, say molecules per unit mass or volume. The left-hand side describes the passive transport along the fluid material lines, whereas the right-hand-side describes diffusion plus the effects of chemical reactions occurring in the fluid moving at the barycentric velocity .

Although thermodynamically favorable, many reactions are limited by the energy barrier which needs to be crossed in order to form the activated state, so that, in the absence of the catalyst, reactions simply do not occur. In its presence, the rate of reaction can increase dramatically, although the change in energy between reactants and products is not affected. In fact, the enzymes affect the rate of a reaction, not its equilibrium. For many reactions involving a single substrate, the rate of consumption of the substrate follows the Michaelis-Menten equation

where is the magnitude of the rate of disappearance of the substrate (or reactant) and is the Michaelis constant, the concentration of the substrate at which the reaction attains half of its maximum rate.

Another reaction term particularly relevant to populations biology is the logistic expression

where describes the (Malthusian) growth rate and characterizes the non-linear decay of the species due to competition for, say, space and/or food. The ADR class falls naturally within the formalism, by simply defining a LB distribution for the concentration, such that and specifying the following local equilibria:

where both fluid and lattice speeds are normalised by the sound speed . To be noted that in the above is the fluid speed, which does not match the current , because ADR’s conserve mass but not momentum. The LB is then particularly well suited to solve the ADR equations in complex geometries, such as those occurring in morphogenesis Ponce Dawson et al. (1993); Ayodele et al. (2011), heterogeneous catalysis Falcucci et al. (2016) and related phenomena.

iv.2 Charged Fluids

Electrostatics plays a vital role in biological processes and requires handling electrolytic solutions and charged solutes typically in flow conditions, the realm of the so-called electrokinetics. Solutes range from strongly charged molecules, such as DNA, to weakly charged macromolecules, such as proteins, where partially screened mobile charges and charged surfaces are ubiquitous conditions of biological settings. Prominent examples include i) the stability of proteins as a function of pH and ionic strength, such as salting-in and salting-out effects due to the interplay of protein charges with the aqueous/saline environment; ii) protein-ligand association processes, including enzymatic allostery modulated by salt-bridges, salt-bridges in virus assemblies and thermal stability; iii) membrane proteins and their specific electrostatic properties, whereby the collective charge of the intracellular residues tend to be more positive as compared to the extracellular ones (the so-called “positive inside” rule), and so on. The associated biological flows are driven by electro-osmosis and electrophoresis. Well-known examples of those processes include i) ion channels and the gating of ions across the cell membrane that regulate electrical signalling in secretory and epithelial cells, as much as the cell volume; ii) the electrokinetic flow due to the glycocalyx layer covering cells, a polyelectrolyte exhibiting a surface charge and responsible for an electrochemical gradient that regulates mechanotransduction; iii) the mitochondria and the chloraplast, where proton gradients generate a chemi-osmotic potential, also known as a proton motive force, for the synthesis of ATP.

Although the general principles of electrokinetics are fairly well understood Masliyah and Bhattacharjee (2006), methods that translate these principles into accurate numerical predictions are still in their infancy. In fact, deriving the interactions between charged solutes and the solvent requires computing the interactions of a large number of molecules and the averaging of these over many solvent configurations. Such daunting computational requirements are partly overcome through the continuum mesoscopic approach, whereby solvent and ions are described in the continuum and in a pre-averaged sense Capuani et al. (2004); Reboux et al. (2006); Wang and Kang (2010).

The LB framework solves complex electrokinetic problems through an efficient formulation of the electrolytic solution as a multi-species problem: one species for the neutral solvent, water, and two for the positively and negatively charged ionic components. The Poisson equation provides the solution for electrostatics and the self-consistent forces for the transport of the fluid species, in the so-called Vlasov-Poisson approximation.

Let us briefly survey the LB method for charged multi-component systems. The fluid mixture is composed by three sets of populations labelled by index , two ionic components with charges , being the proton charge, density and velocity . Given the barycentric velocity , the species relative velocity is . The aqueous medium accommodates solute biomolecules, with the i-th particle having position and valence . The electrostatic potential is obtained by solving the Poisson equation,

(28)

where the medium has dielectric permittivity , duly complemented by boundary conditions of the Dirichlet or Neumann kind. For insulating confining walls, where the surface has local charge density , this reads , where is the unit vector normal to the surface.

The dynamics of each species follows the evolution equation:

(29)

where the Maxwellian equilibrium for mixtures is given by Marconi and Melchionna (2011a, b),

(30)

and the force term

The local self-consistent forces are , where is the drag force exerted on species , resulting in a cross-diffusion coefficient , where is a relaxation frequency.

So much for the governing equations. However, a hallmark of charged systems is that local electrostatic forces can be dramatically intense, often exhibiting rapid spatial modulations of the electrolytic densities, as in the presence of double layers.

Such occurrence can lead to severe numerical instabilities in methods based on the direct solution of NS equation, as a consequence of local violations of the Courant-Friedrich-Levy stability condition Courant et al. (1928). Thanks to its inherently small time step (in macroscale units), LB can often handle stiff forces without loosing stability. However, strong polyelectrolytes are problematic to handle and cannot be treated directly, thus some sort of charge rescaling is required Datar et al. (2017). Despite these liabilities, LB simulations of electrolytic systems show stable behaviour over a wide range of solute charges and molarities, in particular in the sub-molar range that covers a large portion of biological conditions.

iv.3 Flows far from equilibrium

Most biological systems operate far from equilibrium, i.e., they draw energy from their environment and dissipate heat back to it, thereby lowering their own entropy content at the expense of the environment. This is the operating principle of the so called “dissipative structures”, a cornerstone of non-equilibrium thermodynamics Prigogine (2017). In the process of dumping entropy to the environment, they manage to migrate across a sequence of different non-equilibrium steady-states (NESS), in which they deliver different functions.

Examples of NESS in biology include molecular machines, cells in motion, metabolic pathways and many others. For instance, the kinesin protein walks along the microtubule, carrying cargos from one part of the cell to another, by absorbing energy from ATP hydrolysis and converting chemical energy into mechanical work, of which is used for motion and the rest is dissipated to the surroundings. Off-equilibrium conditions imply that biological agents exchange mass, momentum, energy and entropy with the environment through transient or steady currents and fluxes. Clearly, such complex space-time dependent network of currents and fluxes, typically operating on a broad spectrum of concurrent space and time scales, needs to be captured by the numerical approach Takahashi et al. (2005). Whence a major need of multiscale methods.

Indeed, these currents and fluxes typically occur through thin (atomic scale) interfaces, which resolve the tension between competing mechanisms, say chemical reactions and molecular diffusion, through a sudden spatial transition between distinct bulk phases, say the liquid and vapour in a multiphase flow.

For the case of a diffusion-reaction system, the width of the transition region can be estimated as

where is the diffusion coefficient and is a typical chemical timescale. By expressing , with being the mean-free path and the collisional timescale, we obtain . This shows that, unless chemical reactions are much slower than inert collisions (slow-chemistry), the interface width is comparable with the molecular mean free path, or shorter (fast-chemistry).

Another example in point are foams and emulsions, i.e. droplets (bubbles) of liquid (vapor) dispersed in a continuum liquid phase, typically water. In these multiphase systems, the transition between the dense and light phases is controlled by the surface tension which is, in turn, dictated by molecular interactions, notably the strength of the potential and its spatial range. That leads to interface widths of the order of the spatial range of such interactions, typically nanometers or below. A typical estimate of the interface width is as follows:

where denotes the surface tension. Typical values, in MKS units are and , deliver , i.e a fraction of nanometer.

As shown in Fig. 9, such nanometric interfaces challenge the low-Knudsen assumption which lays at the foundations of the hydrodynamic description. In fact, given that the interface thickness is comparable with the molecular mean free path, the result are local Knudsen numbers of the order unity. In the last decade, a number of technical extensions of the original LB method have been developed with the aim of gaining insights into these complex non-equilibrium interfacial phenomena (and yet, much remains to be done).

Figure 9: A sketch of a high-density/low-density diffuse interface of width , as typical for many biosystems (dotted line), for different values of the mean free path, of length . Each solid line shows a different profile of the interface, depending on its intrinsic structure. Here we assume that the wavelength of the density profile coincides with . In the limit (solid line, bottom), non-equilibrium effects are negligible and hydrodynamics holds. On the other hand, when becomes comparable to the interface width (solid line, top), non-equilibrium effects can no longer be neglected and higher order kinetic moments must be accounted for.

Here, we limit ourselves to summarize the main upgrades, namely: i) Higher-Order lattices, ii) Kinetic Boundary Conditions, iii) Regularization techniques.

iv.3.1 Higher-Order Lattices

As discussed in the initial part of this paper, hydrodynamics represents the “infrared” limit of kinetic theory, whereby all macroscopic heterogeneities live on much longer scales than the molecular ones (hydrodynamic transport regime). From the formal viewpoint, this means that the Boltzmann distribution is accurately described by its lower-order moments, typically density (order ), current (order ) and momentum-flux (order ).

In the hydrodynamic transport regime, all higher order moments (non-equilibrium excitations) are directly enslaved to the low order ones, hence they have no independent dynamics of their own.

Far from equilibrium, where strong inhomogeneities may persist down to near-molecular scales, such low-order picture breaks down, and more kinetic moments concur to define the motion of matter beyond the hydrodynamic regime.

Without entering details, it is intuitively clear that this far from equilibrium regime requires the use of higher order lattices, securing the proper recovery of the correspondingly higher order moments. A formal theory of LB beyond NS, based upon higher order lattices, was laid down by X. Shan et al Shan et al. (2006). Those authors developed the discrete analogue of Grad’s expansion in Hermite polynomials, with full details on its specific implementation on a series of higher order lattices associated with different numerical quadrature rules. The first observation as compared to Grad’s 13 moment formulation is that the scheme provides a larger set of degrees of freedom Grad (1949). The discrete speeds of the corresponding lattices are typically of order or more, hence many more than the Grad’s moments Meng and Zhang (2011); Zhang et al. (2006).

iv.3.2 Kinetic Boundary Conditions

Higher-order lattices offer room for extra-moments, but this is not sufficient per-se, unless suitable boundary conditions are formulated near solid boundaries. This is where the lattice formulation makes a distinct contribution: while it appears very hard to devise well-posed boundary conditions for the kinetic moments, in fact complicated nonlinear tensors, lattice formulations lend themselves to conceptually transparent and mathematically well-posed formulations. The reason is always the same: the information moves along straight lines and boundary conditions can be formulated in terms of mechanical relations between the outgoing (fluid-to-wall) and the incoming (wall-to-fluid) discrete distributions. Non-equilibrium flows exchange momentum with the solid walls along both tangential and normal directions: describing this exchange is the mandate of Kinetic Boundary Conditions. Empirical forms adapting this constraint to the lattice environment have been formulated in terms of accommodation coefficients Succi (2002b). Following Maxwell, the idea is that molecules impinging on the wall loose track of their incoming speed Maxwell (1878). Consequently, they are re-injected into the fluid along a random direction, with a velocity drawn from a Maxwellian at the local wall speed and temperature. Albeit handy, these accommodation schemes remain empirical in nature. A more satisfactory formulation was developed by Ansumali and Karlin, who basically expressed the accommodation coefficients in terms of the outgoing (fluid-to-wall) distribution functions and their equilibrium version, thus providing a closed and consistent recipe Ansumali and Karlin (2002). The Ansumali-Karlin boundary conditions exhibit a number of appealing properties. First, they preserve the positivity of the distribution at boundary nodes, i.e., if the incident distribution is positive, the reflected one is guaranteed to be positive too. Second, they readily extend to general wall scattering kernels, such as those allowing a blend of slip and reflection. At the time of this writing, these kinetic boundary conditions are the tool of the trade for finite-Knudsen LB simulations.

iv.3.3 Regularization

Regularization is a general and powerful idea across many fields of theoretical physics, to remove various forms of divergences and singularities which arise whenever a given description/theory fails to capture the physics in point.

Kinetic theory is no exception. Indeed, it is known since long that post-hydrodynamic equations beyond the NS level, suffer a number of problems, primarily short-scale linear instabilities~Bobylev (1982). Many regularization schemes have been proposed ever since to tame such instabilities Struchtrup and Torrilhon (2003, 2007). Regularization procedures have been (re)-discovered only recently by the LB community Latt and Chopard (2006); Montessori et al. (2014) and it is not yet clear how they relate to the corresponding counterparts in continuum kinetic theory. LB regularization consists of filtering out the contribution of non-hydrodynamic moments (ghosts) to the hydrodynamic ones: mass, momentum and momentum-flux.

Let us dig a little bit deeper into the subject.

The standard LB scheme in BGK form reads (time-step made unit for simplicity):

(31)

The actual distribution can be split as follows:

where the hydro-component collects terms up to third order Hermite polynomials associated with mass, momentum, momentum-flux and energy flux, while collects all higher order terms, transport plus genuinely kinetic fields with no immediate macroscopic interpretation (ghosts in LB jargon). Formally, one defines a regularisation operator , projecting the actual distribution onto the hydrodynamic subspace, i.e., , that is .

Both hydro and ghost terms further split into equilibrium and non-equilibrium components, the ghost equilibrium being identically null by construction. Thus, what remains to be filtered out is just the ghost non-equilibrium, which is constantly revived at each free-streaming step, the non-equilibrium engine.

By applying the regularisation projector to the right-hand-side of the BGK equation, an operation corresponding to filtering out ghost components after streaming, we obtain the following Regularized-LB:

(32)

The Regularized-LB has recently made proof of providing significant benefits in terms of improving the stability of the LB scheme under finite-Knudsen conditions.

The “post-hydro” LB literature is vast and growing, but not conclusive yet. In particular, it is not clear how the three main ingredients mentioned above should be combined in order to obtain correct finite-Knudsen behavior.

A very valuable step in the comprehension of LB post-hydro capabilities was provided by Ansumali, Karlin and coworkers, who discovered exact solutions to the hierarchy of nonlinear LB kinetic equations for stationary planar Couette flow at non-vanishing Knudsen numbers Ansumali et al. (2007). By using a 16-speed two-dimensional lattice and kinetic boundary conditions, these authors have derived closed-form solutions for all higher-order moments, and solved them analytically.

The results indicate that the LB hierarchy with larger velocity sets does indeed approximate kinetic theory beyond the NS level. If only for a simple set-up, those exact solutions indicate that LB equipped with kinetic boundary conditions is able to carry quantitative non-hydrodynamic information, hence it can be regarded as a kinetic closure in its own right.

It thus appears that the extension of the LB formalism to higher order lattices can result into an effective tool to probe deeper into non-equilibrium regimes beyond the hydrodynamic description Gan et al. (2015).

iv.4 Fluctuating Lattice Boltzmann

For the case of nanoscale flows, reproducing thermal behavior implies that the LB must incorporate the effects of statistical fluctuations.

To achieve this goal, following in the footsteps of Landau-Lifshitz fluctuating hydrodynamics, Ladd added a source of random fluctuations of the momentum-flux tensor Ladd (1993), namely:

(33)

where is the shear tensor and is the amplitude of the fluctuations.

The latter must be tuned so as to comply with the Fluctuation-Dissipation-Theorem. The resulting fluctuating NS equation reads

(34)

where is the fluctation-free momentum-flux tensor and is the body force acting on the fluid.

The LB scheme is modified accordingly, by adding a stochastic source to the right hand side:

(35)

as sketched in Fig. 10. The source term is local in space and time and acts at the level of the stress tensor and non-hydrodynamic modes, with no effect on mass and momentum conservation.

cd

Figure 10: A sketch of the fluctuating distribution, (wiggly line), featuring rapid oscillations with respect to the non-fluctuating one, (smooth line). Fluctuating hydrodynamics emerges as a consequence of the stochastic component of the distribution function.

In actual practice, the source term is constructed via a set of the lattice eigenvectors with orthonormal according to the scalar product . In the D3Q19 scheme, the eigenvectors correspond to the kinetic moments: is relative to density, to current, to the momentum flux tensor, the remaining eigenvectors to non-hydrodynamic modes. The stochastic forcing reads as follows:

(36)

where is a set of random numbers with zero mean and unit variance. Given the fact that the thermal velocity is fixed and equals the underlying lattice speed , the thermal mass is chosen in such a way as to obtain the thermal fluctuations according to .

The stochastic forcing has been subsequently improved so as to produce consistent fluctuations at all spatial scales, in particular at short distances where the effect on the translocating molecule is critical Dünweg et al. (2007); Adhikari et al. (2005). The fluctuating LB passes a number of litmus tests, particularly the compliance of velocity-velocity and force-force autocorrelation functions with the principle of stochastic kinetic theory, in particular the fluctuation-dissipation theorem.

iv.5 Lattice Boltzmann for non-ideal fluids

Boltzmann originally derived his equation under the assumption of diluteness, whereby molecular collisions take place as zero-ranged, instantaneous events, leading to large scattering angles. This rules out long-range, soft interactions giving rise to small-angle deflections. Such interactions, however, are of utmost importance for biological applications, since soft interlayers based on membranes and biopolymers define the spatial boundaries between different phases in biological systems. In aqueous media containing various ions, interactions are governed by a complex interplay of generic and specific interfacial interactions, typically controlled by dispersion and electrostatic forces, and locally shaped as hydrogen bond networks, disulphide bridges, hydrophobic, entropy-induced interactions and so on.

Soft-core interactions can be included in the kinetic equation in the form of effective one-body forces, resulting from the collective interaction of a representative particle with the self-consistent environment (bath) due to all other particles. Formally, such effective one-body force takes the following form

(37)

where is the (unknown) two-body distribution and is the two-body atomistic potential. Following a customary practice, one writes , where , the two-body correlation function, collects the two-body physics. The exact form of the correlation function is known exactly only in a few precious cases; yet one can introduce several useful ansatzs which turn the formal expression 37 into an operational one. A typical ansatz, which has proven very fruitful for LB modelling of non-ideal fluids, looks like follows:

(38)

where is a model Green function and is a local functional of the fluid density . The above expression is quite general and allows a wide degree of latitude in modeling non-ideal fluid interactions.

iv.6 Pseudo-potential models

The most popular lattice transcription of the above, due to Shan and Chen Shan and Chen (1993, 1994), reads as follows:

(39)

The sum runs over the prescribed set of neighbors, typically the first Brillouin region in the original version, and subsequently extended to the second or even the third one. Typically, all discrete speeds in the same Brillouin region share the same , so that in the original Shan-Chen (SC) formulation, there is just one coupling strength (see Fig.11). This is nonetheless sufficient to generate the main ingredients of non-ideal fluids, namely a non-ideal equations of state supporting phase-transitions, as well as surface tension. Moreover, this limitation can be readily lifted by extending the Shan-Chen interaction beyond the first Brillouin region, thereby explicitly accounting for both repulsive and attractive interactions.

The Shan-Chen expression delivers a non-ideal equation of state of the form

(40)

By choosing the generalized density in the form:

it is readily checked that the Shan-Chen fluid becomes critical for at a critical density . Note that even though the interaction is purely attractive (), no pile-up instabilities take place because the force becomes increasingly faint as the density increases. This is the reason for using the generalized density instead of the physical one . This is a very expedient trick to trigger phase-transitions without any repulsive force. The translation from the force to the LB source proceeds through a systematic expansion in lattice Hermite polynomials. To leading order , but higher order terms, at least quadrupole ones, are definitely needed to obtain accurate results.

Figure 11: Sketch of the Shan-Chen piece-wise linear potential, extending over the first Brillouin lattice cell at a distance . The potential mimics the attractive tail of a van der Waals potential (dashed line) while the repulsive one is quenched to zero. Notwithstanding the absence of hard-core repulsion, the Shan-Chen potential does not cause any unstable density pileup because, in the high density limit, , the generalised density flat-tops to a constant , yielding a zero gradient, hence zero force .

More generally, thick interfaces lead to large values of the Cahn number, namely the ratio of the interface width to a typical mesoscopic scale, say the droplet (bubble) diameter:

For most applications is in the range of tens, up to one hundred, microns, leading to very small Cahn numbers, . Replicating this number in LB simulations is totally unviable, for it would amount to placing lattice spacings across the particle diameter. LB simulations are forced to operate at much higher Cahn numbers of the order of , which means that potential inaccuracies due to finite Cahn number effects need to be carefully monitored.

Lack of sufficient symmetry also leads to the appearance of spurious currents, which may seriously affect the physical results whenever the density ratio between the light and dense phases is above 30 – 50, a problem which is further exacerbated in the case of fast-moving interfaces. Finally, the model is not thermodynamically consistent, as it is not derived from a mean-field free-energy functional. That limitation is, however, less serious than it seems, since the method is, actually, equipped with a quasi-free energy Sbragaglia et al. (2009). As a matter of fact, more disturbing, instead, is the fact that the Shan-Chen state equation features a sound speed smaller in the liquid with respect to the vapour phase, with significant consequences on the stability of the light phase across the interface.

Most of these limitations have been significantly mitigated, if perhaps not fully resolved, by subsequent developments. Among others, particularly noteworthy for the simulation of soft-glassy materials, is the so-called multi-range pseudo potential method.

The main idea is to augment the original Shan-Chen short-range attraction in the first Brillouin region (belt, in LB jargon), say the D3Q27 lattice in three dimensions, with a repulsive interaction acting on the second Brillouin region, i.e., including discrete velocities up to . The main advantage of this formulation is that a proper tuning between the first-belt attraction and the second-belt repulsion, makes possible to achieve smaller values of the surface tension, thereby permitting to sustain long-lived, multi-droplet configurations with a highly complex interfacial dynamics. Despite its many limitations, the SC method remains the most popular version of LB for non-ideal fluids, mostly on account of its simplicity and transparency.

However, for applications leading to complex interfacial dynamics, more advanced schemes are required.

iv.6.1 Free-Energy models

Another successful route to LB schemes for non-ideal fluids is (lattice) Density Functional Theory (DFT), whereby the non-ideal physics is encoded within a free-energy functional of the form:

(41)

where is a local functional of the fluid density and its gradients Swift et al. (1996).

The former encodes the non-ideal equation of state, while the latter describes the effect of surface tension.

Variational minimization over density changes delivers the equations of motion of the non-ideal fluids, governed by the Korteweg pressure tensor:

(42)

where

is the non-local pressure, consisting of the bulk contribution , fixing the equation of state, plus an interface contribution associated to surface tension. In the above, stands for the space derivative along direction and is a tunable coefficient controlling the surface tension.

From the practical standpoint, the free-energy formulation leads to non-local equilibria, which involve second order derivatives of the density field, thus taxing the simplicity of the method, and sometimes its stability as well, as compared to the pseudo-potential method. Nevertheless, the free-energy method has found broad use for many applications such as microflows over geometrically or chemically patterned surfaces and various types of droplet motions Liu et al. (2014).

Several variants have been developed over the years, which have considerably improved over the original versions, especially in terms of reaching higher density ratios without compromising numerical stability. Among others, high-order finite-difference schemes Lee and Lin (2005) and mergers with flux-limiting methods Sofonea et al. (2004). These methods are currently used to simulate a variety of multiphase and multi-component fluid applications, such as Rayleigh-Taylor instabilities, droplet collisions, cavitation and free-surface flows. For a review, see the recent books Succi (2018); Montessori and Falcucci (2018) and references therein.

iv.6.2 Chromodynamic models

Another class of LB methods for non-ideal fluids which has been revamped in the recent past is the so-called Color-Gradient model. Here, the two phases/components, call them Red and Blue for convenience, segregate based on a top-down “colouring” rule, which sends the particles along the color gradient, namely Red towards Red and Blue towards Blue, thereby triggering an instability leading to interface formation Gunstensen et al. (1991).

The recoloring amounts to correcting the the Red and Blue populations as follows:

(43)

where is the mass density fraction, is the cosine of the angle between the color gradient and the -th discrete velocity. In the above , where star indicates the population after the application of the non-ideal force stemming from surface tension, and , where superscript denotes the equilibria at zero flow. Finally, is a free parameter controlling the strength of the recoloring procedure, to all effects and purposes an anti-diffusive operator. The crucial term is the second one at the right-hand side of eq. (43), which, by construction, is active only at the interface between the two components. For more details see Leclaire et al. (2017) and Montessori et al. (2018a).

The interface is then stabilised by means of a chromo-dynamic force proportional to the amplitude of the color gradient, but opposite to it, so as to level out the deficit of one species over the other (color gradient). Although essentially rule-driven, modern variants of this scheme have proved capable of accessing parameter regimes which appear to be off-limits for Shan-Chen schemes and extensions thereof, as well as for Free-Energy methods.} For instance, such methods have been recently applied to the design of microfluidic devices for droplet generation, such as flow-focusers and step-emulsifiers \cite{AM2018a,AM2018b}.

The interface is then stabilised by means of a chromo-dynamic force proportional to the amplitude of the color gradient, but opposite to it, so as to level out the deficit of one species over the other (color gradient). Although essentially rule-driven, modern variants of this scheme have proved capable of accessing parameter regimes which appear to be off-limits for Shan-Chen schemes and extensions thereof, as well as for Free-Energy methods.

The recoloring amounts to correcting the Red and Blue populations as follows: where is the mass density fraction, is the cosine of the angle between the color gradient and the -th discrete velocity. In the above , where the star indicates the population after the application of the non-ideal force stemming from surface tension, and , where superscript denotes the equilibria at zero flow. Finally, is a free parameter controlling the strength of the recoloring procedure, to all effects and purposes an anti-diffusive operator. The crucial term is the second one at the right-hand side of eq. (43), which, by construction, is active only at the interface between the two components. For more details see Leclaire et al. (2017); Montessori et al. (2018a).

The interface is then stabilised by means of a chromo-dynamic force proportional to the amplitude of the color gradient, but opposite to it, so as to level out the deficit of one species over the other (color gradient). Although essentially rule-driven, modern variants of this scheme have proved capable of accessing parameter regimes which appear to be off-limits for Shan-Chen schemes and extensions thereof, as well as for Free-Energy methods. For instance, such methods have been recently applied to the design of microfluidic devices for droplet generation, such as flow-focusers and step-emulsifiers Montessori et al. (2018b, c).

iv.6.3 Entropic models for multiphase flows

The lattice Boltzmann method can also be formulated by minimizing a suitable lattice H-function of the form (Kullback-Leibler entropy):

where are suitable lattice weights. The resulting scheme takes the usual form of a standard LB, with the crucial twist that the relaxation time is adaptively adjusted in such a way as to secure compliance with the second principle of thermodyamics, namely:

Leaving a detailed description to the original literature, here we just mention that compliance with the second principle translates into a significant enhancement of numerical stability Karlin et al. (1998). For this reason, ELB has found profitable use for the simulation of low-viscous flows typical of fluid turbulence.

In a recent time, the ELB has been extended to the case of multiphase flows, and shown to provide stability benefits also in the viscous regime of relevance to many biological applications Wöhrwag et al. (2018); Montessori et al. (2017). Although these developments are too recent to permit a solid prononciation, the results appear very encouraging and leave hope that the entropic version of LB may become a major player in the field for the years to come.

V Fluids and Particles: The Lattice Boltzmann - Particle Dynamics Scheme

The multiphase LB schemes discussed in the previous Section have generated a mainstream of applications in soft matter research, since they permit to deal with flows of major dynamic and morphological complexity, such as foams and emulsions. However, they are unsuited to handle rigid bodies suspended in the continuum phase, nor can they describe in detail mechanical properties of deformable objects, say membranes, vesicles, cells and other biological bodies. To this purpose, LB needs to be explicitly coupled to particle methods, tracking the dynamics of the biological bodies immersed in the flow.

Indeed, most flows of biological interest consist of biological bodies of assorted nature: cells, polymers, proteins, floating in a fluid solvent, typically water. Such flows often operate at low, often near-zero, Reynolds number Purcell (1977), but this does not mean that hydrodynamic interactions (HI) can be neglected. To the contrary, HI have repeatedly been shown to accelerate a variety of nanoscale biological transport processes, such as biopolymer translocation across nanopores, amyloid aggregation of proteins in the cell and other related phenomena. This explains why the combination of LB with particle dynamics has made the object of extensive research and the development of multiple simulation schemes.

v.1 The extended particle model (EPM)

Flows with suspended objects have been modelled since the early days of LB research, starting with the trailblazing work of A. Ladd Ladd (1994a, b). Ladd’s original method consists of tracking the motion of rigid spherical bodies under the impact of the surrounding solvent hitting the surface of the body. The scales pertaining to body and solvent are fairly separate, since the former is much larger and heavier than the latter, hence the mass ratio serves as a suitable scale separator.

In the EPM, the exchange of momentum between particles and LB fluid is a boundary-collision method whereby the suspended body interact with each other only through the mediation of local collisions with the solvent molecules (LB). The solvent-body collision is conservative, the momentum lost (gained) by the solvent to the body is gained (lost) by the body on the solvent. The method therefore is based on the local exchange of momentum by computing the amount of momentum that every population hitting the body surface exchanges with the latter. Being the total force and torque the sum of fluid-body momentum exchange (or drag interactions), direct particle-particle mechanical interactions or, in case of microscopic conditions, stochastic forces due to the Brownian motion, the body linear momentum and angular momentum obey the classical equations of motion

(44)

It is worth mentioning that the fluid-particle coupling is hydrokinetic in nature because it is treated at collisional level rather than at hydrodynamic level, so that, hydrodynamic forces such as long-time tails, naturally emerge at times larger than the LB marching one. The approach takes naturally into account the solvent fluctuations, if present, since the latter are transmitted across the body surface analogously to the drag forces (see Fig.12).

Since the solvent-body collisions are conservative, no extra stochastic source is needed on the body side.

Ladd’s method finds its best use in simulating colloidal suspensions of rigid particles and in this sense has found a comparatively limited application in the biological context. However, it provides a first and reliable example of coupling between LB and suspended bodies. At the same time, the method shows some inaccuracy due to the lattice discreteness, especially for near-contact colloid-colloid interactions, occurring below the grid scale where the LB method can not resolve the lubrication forces. Some variants have been developed in the literature, making use of grid refinement and/or dynamic interpolations. Unsurprisingly, they only add to the computational complexity of the method, which is comparatively laborious on its own, as it demands full knowledge of the local fluid-solid connectivity, namely the dynamic list of fluid nodes interacting with the solid ones.

Figure 12: Ladd’s method for a suspended (spherical) body: a) the body-fluid interface is tracked by searching the mesh points inside and outside the particle Ladd (1994a). Here a modified bounce-back scheme is applied to account for the local exchange of momentum. The mid-point of the connecting links identifies the stair-cased surface that corresponds to the no-slip condition; b) the particle center moves in the continuum according to the momentum exchange, and the rotational motion is accounted for by computing the total torque.

v.2 The point-particle model (PPM)

A qualitatively different strategy has been proposed by Ahlrichs and Dünweg Ahlrichs and Dünweg (1998). This is a minimal way to embed point-like particles in the LB fluid which, as opposed to Ladd’s method, is entirely force-driven, hence metric instead of topologic. In essence, particles carry a phenomenological friction coefficient (also known as the bare friction coefficient) and the fluid exerts a drag force based on the instantaneous difference between the particle and fluid velocities, reading

(45)

The addition of the point force into the fluid equations introduces a singularity into the flow field, as represented by Dirac’s delta function. This is automatically smoothed out by interpolation procedures, but requires nonetheless due care in the numerical treatment.

On the other hand, the flow field around a finite-sized particle is generated by a distributed force located around the particle position, as shown in Fig. 13.

In the scheme, the flow field is treated as everywhere finite and the force density acts onto the LB fluid at position as

(46)

where is a function interpolating between the particle position and the surrounding mesh node . In its original formulation, the method replaced and by the values of these fields at the mesh node nearest to the particle position.

A refined version resorts to a simple linear interpolation scheme using the mesh points on the elementary lattice cell containing the particle. Denoting the relative position of the particle in this cell by with the origin being at the lower left front edge, one defines , etc., the interpolation weights. The interpolated velocity then reads , where the sum is over the mesh points on the considered elementary lattice cell. The force exerted back to the fluid can then be chosen with the same weight coefficients or by spreading the force equally on the edges of the cell, the corresponding rule to preserve linear momentum being easily derived.

Due to the regularization of the point-force, one should expect that the mobility of the suspended particle is not simply given by the bare friction coefficient , but is somehow renormalized. Indeed, the mobility is given by the sum of the bare mobility and a Stokes-type contribution due to the lattice discretization. The effective friction coefficient relates to the bare one via , where is a geometrical correction coefficient Ahlrichs and Dünweg (1998).

The PPM finds application in simulating microscopic systems, with a stochastic term added on each particle in addition to the fluctuating LB bath. In this way, the particle does not leak energy on average and the fluctuation-dissipation balance maintains a well-defined temperature of the fluid-particle system. Again, to balance and preserve momentum, such stochastic force is also restituted to the fluid with an opposite sign. The presence of long-time tails, that is, the inherent hydrodynamically-sustained motion of the moving particle generating a long-time decay of velocity, has also been observed Ahlrichs and Dünweg (1998). On the downside, the PPM is permeable to fluid momentum, as the local force is unable to create enough resistance to the incoming flow, and the classical Stokes-like picture of the streamlines does not apply.

Owing to its inherent simplicity, the PPM was first recognized to be useful to simulate topologically connected particles, such as polymers, in the presence of hydrodynamic interations. Whenever the embedded particles represent micro/mesoscopic objects, such as atoms or the beads of a polymer, mass diffusion plays a role comparable to mechanical and hydrodynamic forces. Another interesting extension involves coupling PPM with the Shan–Chen multicomponent LB to simulate efficiently complex fluid–fluid interfaces. The idea is to introduce a solvation free-energy for the particle–fluid interaction proportional to the fluid density gradients Sega et al. (2013) and reads as follows

(47)

Such force drives particles towards maximal or minimal density gradients, depending on the sign of the coupling coefficient . The approach is particularly appealing when used in conjunction with multiphase conditions. As a matter of fact, it makes possible to treat different multiphase fluids in the presence of suspended molecules, such as amphiphilic molecules as surfactants or, by adding another level of detail with electrostatics, polyelectrolytes in bicomponent fluids.

v.3 The diffused-particle model (DPM)

An extension of PPM represents particles with finite-size extension, still relying on a force-based mechanism. Given the finite extension, it is well suited for anisotropic particles, providing a very handy computational flexibility for the description of biological suspensions, cellular compartments or even entire cells (recalling that the interior of the cell is anisotropic due to intracellular organelles). The strategy is to consider the particle roto-translational response as originating from the coupling of the finite-size extension of the particle with the fluid momentum and vorticity. The particle is an effective diffused body, with no need to track its boundary to control its coupling with the environment.

The particle shape is described as an ellipsoid having three major radii along the three principal directions, with . To fit in the discrete nature of the lattice, are taken as three integers, one for each cartesian component (such requirement can be lifted but a few interesting properties described below, would be lost). The roto-translation is governed by rigid body dynamics of eqs. 44. The particle rotational state is encoded by the matrix whose rows are three orthogonal unit vectors aligned along the principle axis of the particle, that is, the basis to transform between the laboratory and the moving frames. The roto-translational state is given by the tensorial product

(48)

where

is a shape function having compact support and the normalization property . Typically corresponds to a spherically symmetric diffused particle with a support extending over 64 mesh points. The translational response of the suspended body is designed according to the fluid-particle exchange kernel

where and the hydrodynamic DPM force is obtained via integration over the particle spatial extension. It reads as follows:

where .

Figure 13: Examples of immersed bodies in the LBPD approach. Panels a)-c) show the different coupling methods between particles and fluid, as detailed in the text: a) PPM with point-like particles exchanging information with the nearest grid node, b) DPM for diffused isotropic particles, c) DPM for diffused anisotropic particles. Panels d) shows the timestepping method between particle dynamics (PD) and fluid dynamics (LB). The timestepping scheme of the LB and PD individual components can be either synchronous or asynchronous. Typically, the PD component ticks more frequently than the LB component, allowing for more efficient simulations. Panel e) illustrates a typical configuration of the fluid around a moving particle. Hydrodynamic response manifests itself as Stokes flow field (upper part) and long-range flow structure around a spherical particle or its far-field flow pattern (lower part).

The coupling between the body motion and the fluid vorticity is represented by a rotational kernel

with a rotational coupling coefficient. The corresponding drag torque is

where