GENERIC-compliant simulations of Brownian multi-particle systems: modeling stochastic lubrication

A stochastic Lagrangian model for simulating the dynamics and rheology of a Brownian multi-particle system interacting with a simple liquid medium is presented. The discrete particle model is formulated within the GENERIC framework for Non-Equilibrium Thermodynamics and therefore it satisfies discretely the First/Second Laws of Thermodynamics and the Fluctuation Dissipation Theorem (FDT). Long-range fluctuating hydrodynamics interactions between suspended particles are described by an explicit solvent model. To this purpose, the Smoothed Dissipative Particle Dynamics method is adopted, which is a GENERIC-compliant Lagrangian meshless discretization of the fluctuating Navier–Stokes equations. In dense multi-particle systems, the average inter-particle distance is typically small compared to the particle size and short-range hydrodynamics interactions play a major role. In order to bypass an explicit—computationally costly—solution for these forces, a lubrication correction is introduced based on semi-analytical expressions for spheres under Stokes flow conditions. We generalize here the lubrication formalism to Brownian conditions, where an additional thermal-lubrication contribution needs to be taken into account in a way that discretely satisfies FDT. The coupled lubrication dynamics is integrated in time using a generalized semi-implicit splitting scheme for stochastic differential equations. The model is finally validated for a single particle diffusion as well as for a Brownian multi-particle system under homogeneous shear flow. Results for the diffusional properties as well as the rheological behavior of the whole suspension are presented and discussed.


Introduction
Particles suspended in liquids are ubiquitous in nature, with slurries, gels or pharmaceuticals as some common examples [21]. When the particle size starts to enter the sub-micrometer regime, thermal fluctuations of the surrounding environment start to play a relevant role, leading to an increased diffusional properties, i.e Brownian motion [8]. An accurate and efficient description of Brownian multi-particle systems is critical in many applications and the understanding of their dynamics and mechanical response is both of high scientific importance and industrial relevance [20]. In the context of microrheology [35], for example, it is possible to infer the non-Newtonian properties of suspending fluids and/or the heterogeneity of the microscopic environments, uniquely from the diffusional patterns of immersed nanoparticles. More specifically, in terms of biomedical applications, intra-cellular microrheology [44] can allow to infer altered mechanical properties of the liquid environment which can define novel biomarkers to several pathologies [18]. On the other hand, in the context of material engineering, Brownian suspensions are generally characterized by a fluid's viscosity that drops, a response known as "shear thinning" [29]. That rheology is engineered into a range of consumer products, from shampoos and paints to liquid detergents, to make them to flow easily under a weak stress [43]. Opposite shear-thickening can also occur at very high shear-rates as effects of short-range lubrications [43] and contact forces [33].
When modelling Brownian particles immersed in a liquid medium, hydrodynamics, thermal fluctuations and mass diffusion needs to be consistently taken into account. There is a large number of numerical techniques that have allowed accurate simulations for these systems in the past, for which a complete review is out of the scope of this paper. Generally speaking, these techniques belong to two main classes: (a) implicit-solvent and (b) explicitsolvent methods [25]. In the first class, the solvent medium is not described, whereas the forces acting between suspended particles are obtained from Green's function theory in the over-damped Stokes limit, leading to far-field Oseen or Rotne-Prager-Yamakawa hydrodynamics approximations. Brownian Dynamics [11] or Stokesian Dynamics techniques [4] are two popular examples of implicit solvers. Although being highly efficient, the application of these methods to more complex systems, e.g. arbitrary-shape particles, non-Newtonian solvent media or interaction with complex boundaries, is difficult.
In the second class of methods, as the word suggests, the liquid behavior is explicitly described in terms of extra computational degrees of freedoms. Forces and torques are then simply computed on-the-fly from the fluid phase and exerted on the dispersed solid phase. In this class we finds standard direct numerical simulations approaches for the discretization of the liquid phase such as Finite Elements methods [5], Finite Volumes methods [34], Lattice Boltzmann method [16] or Lagrangian particle-based techniques such as the Smoothed Particle Hydrodynamics [22] and Dissipative Particle Dynamics [15]. In the latter class of explicit schemes, a crucial modelling feature is related to the algorithm's ability to accurately simulate fluid dynamics and thermal fluctuations at the same time, i.e the so-called Landau and Lifshitz Fluctuating Hydrodynamic (FH) regime [17]. Among several Brownian explicit solvers available on the markets, there is only a limited number of tools that have paid special attention to the exact (i.e. discrete 1 ) enforcement of thermodynamic properties. In fact, the direct incorporation of the FH equations in a standard discretization scheme, such as the finite element/volumes method [5] is not straightforward. In order to have the correct level of thermal fluctuations at all resolved scales, the FDT must be fulfilled at the discrete level (not only in the continuum limit) which is not generally the case [34]. In the Finite Volume-based scheme proposed by Donev et al. [7], for example, the Fluctuation-Dissipation Theorem was proved to be satisfied at the discrete level in virtue of the specific construction of the skew-adjoint relation between the discretized gradient and divergence operators.
Here we consider another Lagrangian stochastic particle solver, named Smoothed Dissipative Particle Dynamics (SDPD) [10,12]. In SDPD, Lagrangian elements (i.e. fluid particles) are considered to computationally map the fluid domain, which interact via pair-wise repulsive, dissipative and random forces. The SDPD method is built within the so-called GENERIC framework (General Equation for Non-Equilibrium Reversible-Irreversible Coupling) [14] which allows to easily construct discrete models with a minimal set of geometrical constraints that enforce the Laws of Thermodynamics. In [39] it was shown that a stochastic SDPD formulation can be developed which describes the Fluctuating Hydrodynamics and, at the same time enjoys a number of physical features: (a) it conserves exactly (at the discrete level) the total energy of the system; (b) it has an entropy which, by construction, is a monotonically increasing function of time; (c) it satisfies discretely a Fluctuation-Dissipation Theorem; and (d) it posses fluctuations associated to the hydrodynamics variables which are distributed according to the equilibrium Boltzmann distribution. Moreover, in [39] it was also shown that in the limit of vanishing fluctuations the SDPD method can be interpreted as a specific Lagrangian meshless discretization of the deterministic Navier-Stokes equations, i.e. linking it to the Smoothed Particle Hydrodynamics (SPH) method. The comparison of the GENERIC-derived particle model with SPH is precisely needed to determine the explicit expressions of the amplitudes of the thermal noise which, in turn, produce a consistent SDPD discretization of the FH with prescribed input transport coefficients. All these features make it a suitable numerical tool to study flow problems under Brownian conditions.
In relation to solid particles immersed in liquid, in [3] it was shown that hydrodynamic interactions between suspended solid Brownian particles are correctly captured by the SDPD technique as long as the flow is fully resolved. In more exact terms, by considering two moving solid particles of radii a, placed at a center-to-center distance d, forces and correlations mediated by the interstitial fluid are correctly captured as long as the the numerical SDPD resolution is sufficiently large (i.e. the SDPD fluid particle size dx s = d − 2a, being s the gap length between their surfaces). On the other hand, if the gaps size s → dx the problem is numerically under-resolved and hydrodynamic interactions are incorrect. In order to bypass this problem, in [2,26,41] we have incorporated a pairwise short-range lubrication contributions which "corrects" for the lack of numerical resolution when two solid particles becomes too close (i.e. for s ≈ dx). The method has been shown to accurately reproduce the dynamics of two-mutually interacting non-Brownian particles as well as the dissipative response of a many-particle system [37].
In the present paper, we generalize the above-mentioned methodology to Brownian regimes, by extending the lubrication corrections to the stochastic case. Since lubrication leads to irreversible dissipative forces, they need also to satisfy FDT at the discrete level. We show here how a general stochastic lubrication model can be formulated in a way that complies precisely with GENERIC. The coupled lubrication dynamics is integrated in time using a generalized semi-implicit splitting schemes for stochastic differential equations, in the spirit of [19]. The resulting stochastic SDPD model allows to correctly reproduce fluctuating hydrodynamic interactions in a dense Brownian many-particle systems, from far-field to short-range conditions, in a way that is fully GENERIC compliant. The model is finally validated in a number of Brownian flow conditions including particle diffusion and Brownian suspension rheology, and the results discussed in relation with established data form literature.

The GENERIC framework
In this section we present a brief review of the GENERIC formalism [14]. Let's consider a system with a vector of relevant variables x, a total energy E and a total entropy S. According to the GENERIC formalism, the deterministic part of the dynamics reads Both entropy and energy of the system are functions of the relevant variables. The first term of the right hand side is corresponding to the reversible part of the dynamics of the system, while the second term is describing the irreversible part. L and M matrices must satisfy some properties: L is antisymmetric and M is symmetric and positive-definite. The next degeneration conditions must be held in order to fulfill the first and second Laws of Thermodynamics (Ė = 0 andṠ ≥ 0), respectively: If there were other dynamic invariants I (x) of the system, two more conditions should be considered Furthermore, if thermal fluctuations are relevant, instead of Eq. (1), the system dynamics follow the next stochastic differential equation The stochastic term dx is a linear combination of independent increments of the Wiener process, and satisfies the Itô rule where k B is the Boltzmann's constant and Eq. (5) represents the Fluctuation-Dissipation Theorem. Note that this equation implies that matrix M is symmetric and positive-definite. Another way to write the stochastic term is where superindices refer to tensorial components, and where the noises d W hold the mnemothecnical Itô rule With this definition, matrix M can be written as When thermal fluctuations are considered, there is a stronger requirement than Eq. (2) to ensure that the total energy and the invariants do not change during the evolution of the system: These equations represent a minimal guide to embed 1st and 2nd Laws of Thermodynamics and incorporate thermal noise in new models, and they can be explained in geometrical language. For example, considering the first equation, E(x) = E T can be seen as the hypersurface of points of space defined by coordinates x, where the system has a total energy equal to E T . Derivative ∂ E ∂ x of the total energy E respect to the relevant variables x, is a normal vector to such a hyper-surface, and Eq. (9) requires that infinitesimal movements of the system in this space (e.g. due to thermal fluctuations) must be tangent to the hyper-surface in order to conserve total energy. The same holds for further invariants.

SDPD model for the liquid medium
In this section we summarize the main results of the SDPD method for the thermodynamicallyconsistent discretization of the Fluctuating Hydrodynamics. For a recent overview of the technique the reader is referred to [10]. A set of computational fluid particles labelled by Latin indices i, j = 1, . . . , N , will be considered here which is homogeneously distributed in the liquid domain. Positions are denoted by r i , momenta by p i and entropies s i , masses m i and temperatures T i . In [12] it was shown that the geometrical set of equations for the irreversible dynamics of a discrete particle system obtained via GENERIC can be expressed as: being dt a discrete time step, The term T i j is the difference of temperatures T i − T j and a i j , b i j and c i j are suitable coefficients that will be linked to fluid properties. The stochastic terms have been postulated to be of the form where dp i j = −dp ji , dq i j = −dq ji and d˜ i j = −d˜ ji represent spontaneous momentum and heat exchanges between particles i and j caused by a thermal fluctuation and are defined as being D the number of dimensions and 1 the identity matrix. Term dV i is an independent increment of the Wiener process and d W i j is a matrix of independent increments of the Wiener processes which satisfy the following Itô mnemotechnical rules where δ is the Kronecker delta and Greek superindices here refer to tensorial components. Finally, the matrix d W i j in Eq. (12) is a traceless symmetric matrix given by (14) and the terms A i j , B i j and C i j appearing in Eq. (12) are the thermal noise amplitudes. It is simple to show that, by defining the irreversible matrix according to the dyadic product Eq. (5) with the general noise terms in Eq. (12), the discrete irreversible dynamics obtained via GENERIC is precisely that reported in Eq. (10). However, this is is true only if the following matching of coefficients is considered which correspond to the discrete SDPD Fluctuation Dissipation Theorem. In virtue of the compliance with GENERIC, this geometrical structure of the particle equations enforces automatically: (a) the conservation of the total and local linear momentum; (b) the conservation of the total energy; (c) the monotonic time increase of the total entropy; (d) the Boltzmann distribution of the thermal fluctuations; and (e) the validity of the Fluctuation Dissipation Theorem at the discrete level. For further details the reader is referred to [10,12]. Note that Eq. (10) describes a general stochastic dissipative model for "particles" where no link has been established yet with fluid mechanics. Such a general structure can be compared to different particle-based discretizations of the non-isothermal Navier-Stokes equations, such as the Smoothed Particle Hydrodynamics [10,12] or the mesoscopic Voronoi method [32] in order to obtain specific expressions of the amplitudes of the thermal noise and, consequently, a related irreversible dynamics. If a SPH model is chosen, the governing equations are given by [12] where a and b are transport coefficients related to the shear and bulk viscosities η s and ζ s as The terms P i correspond to pressure contributions associated to the particle i and can be defined based on local fluid densities (see Eq. (21) below). Note that the pressure forces appearing in the first summation on the r.h.s. of Eq. (16) represent a reversible contribution and they can be also derived through the GENERIC route by considering the generalized derivative of the energy function in Eq. (1). For further details the reader is referred to [10]. For the moment, it is sufficient to notice that these pressure forces are anti-symmetric by swapping particle index (i ↔ j) and therefore they preserve the linear and angular momentum locally. The coefficient κ is the thermal conductivity coefficient of the fluid. The SPH kernel W i j = W (r i j ) is a normalized interpolation bell-shaped function with compact support of size r cut , and W i j = ∂ W (r , r cut )/∂r | r =r i j is its derivative. The quantity d i is the particle number density, related to the particle volume V i as The comparison with the GENERIC model allows to find the noise amplitudes which in this case read specifically Summarizing: Eq. (16) with the stochastic terms provided in Eqs. (11), (12) and the noise amplitudes given in Eq. (19), represent a GENERIC-compliant SDPD discretization of the the non-isothermal fluctuating Navier-Stokes equations [17]. It should be highlighted that Eq. (16) do not generally conserve the total angular momentum, due to the presence of the term proportional to v i j in the viscous inter-particle forces. Although it is possible to consider additional evolution equations for the particle spins [24], another simple way to restore angular momentum conservation is by an appropriate selection of the bulk viscosity to constrain the inter-particle forces to be central. This can be done by taking ζ s = 5 3 η s . With this specific choice, the final governing equation for the momentum yields In the isothermal conditions considered in this work, Eq. (20) will be considered, however more general non-equilibrium conditions can be modelled. Finally, the isotropic particle pressure P i is computed by using an equation of state [22] where ρ i = m i d i is the local mass density, and p 0 and ρ 0 are reference pressure and reference density respectively. Speed of sound c s = √ γ p 0 /ρ 0 is chosen sufficiently large compared to other characteristic velocities, therefore enforcing a weakly-compressible regime of operation.
Finally, temporal integration of the SPH equations for the matrix fluid is performed with a modified velocity-Verlet scheme [15]. For the weighting function W , the present work adopts a quintic spline kernel [23] with cutoff radius r cut = 4 dx (dx being the mean SDPD particle spacing) [9].

Fluid-particle interaction: long-range fluctuating hydrodynamics
Fluid-structure interaction with suspended inclusions of arbitrary shapes can be modelled using boundary particles located inside a solid region [3,41] (see Fig. 1). The no-slip boundary condition at the liquid-solid interface is enforced during each interaction between fluid particle and boundary particle by assigning an artificial velocity to the boundary particle, which satisfies zero interpolation at the interface [23]. Additional artificial penetration of SDPD fluid elements into the solid structure is avoided by means of specular reflection boundary conditions discussed in Sect. 2.6. Finally, once all fluid-boundary forces are defined, a total force F sph α and torque T sph α exerted by the surrounding fluid on a given solid sphere labelled α = 1, . . . , N c can be calculated and the corresponding coordinates updated as a rigid-body translation/rotation [37,38] asṘ where M α is the mass of the solid particle α, R α its center of mass position, V α its linear velocity, ω α the angular velocity and I α its moment of inertia. In this way, the interaction of a suspended solid particle with its surrounding fluid medium, as well as with other solid particles located at "sufficiently large" distances, is accurately reproduced by the explicit SDPD stochastic model. The terminology "sufficiently large" here is considered in the sense that the local flow in the gaps between solid particles needs to be fully-resolved numerically, that is, the surface-to-surface gap size s must be sufficiently larger than the resolution length dx in the SDPD model. Although this computational requirement does not represent an issue for dilute systems (where solid inter-particle distances are typically large), it might pose notable numerical difficulties in dense many-particle systems where particles are frequently in near-contact.

Solid particle-particle interaction: short-range stochastic lubrication
As mentioned in the previous section, long-range stochastic hydrodynamic interactions between suspended solid particles are mediated by the SDPD fluid and are accurately described by Eq. (22). As discussed in [2,37], in order to reproduce accurately the short-range hydrodynamic behavior however, viscous lubrication as well as short-range inter-particle repulsion models need to be considered. Short-range inter-particle repulsion forces prevent spurious inter-particle overlap and can be physically associated to molecular Van-der-Walls or steric surface interactions [21].
More tricky is the way in which short-range lubrication forces between suspended solid particles are formally introduced, since they are irreversible and therefore they need to follow the global GENERIC structure. As discussed above, when a fixed numerical resolution dx is considered , there will always exist a limiting small distance between close solid particles for which the shearing/squeezing flow in that liquid gap will be under-resolved. This problem can be solved in two ways: either (i) by dynamically and adaptively refining the local SDPD resolution in the gaps, as done in [36], with associated computational costs; or (ii) by introducing a sub-particle model handling the unresolved short-range hydrodynamic interactions implicitly. In the Stokes limit usually considered in Brownian systems, it is possible to work out analytical expressions for the lubrication forces acting between pairs of spheres in arbitrary relative motion. In particular, the normal and tangential lubrication forces acting between close spheres (i.e. at distances s a; a particle radius) read where e αβ = R αβ /R αβ is the vector joining the centers of mass of solid particles α and β, V αβ is their relative velocity and s = R αβ − (a α + a β ) is the distance in the gap between sphere-sphere surfaces and a α and a β are the sphere's radii. s n c and s t c are the normal and tangential cutoff distances below which the lubrication corrections are applied and they depend on the SDPD resolution length. Above those distances, lubrication corrections vanish since the hydrodynamics is correctly captured by the explicit SDPD fluid model. Analytical expressions for the scalar functions f αβ (s) and g αβ (s) are given by [37] Note that lubrication forces are irreversible and as such, in order to fulfill the Fluctuation-Dissipation Theorem, an associated fluctuating term needs to be associated in a GENERICcompliant way. In order to propose an expression for it, we rely on the general structure of the GENERIC particle model. Its irreversible part should be given by Eq. (10) which in the case of a solid particle considered here, can be rewritten as where a αβ and b αβ are suitable transport coefficients and, analogously, the stochastic term dP αβ is given by where the noise amplitudes A αβ , B αβ are related to a αβ , b αβ via In the specific case of the lubrication correction between solid particles, the change of momentum due to the lubrication interaction is given by where the last summation represents the overall contribution to the stochastic term on the dynamics of the solid particle α coming from short-range lubrication interactions with (nearly-touching) neighbors. 2 Note that recent experiments [27] showed that the lubrication interactions between two Brownian beads in near-contact follow strictly the expression given in Eqs. (23,24), which were previously validated in our deterministic model in [37]. For a cleaner comparison with Eq. (62), Eq. (28) can be re-written as where Eq. (23) have been used. As a result, the general terms a αβ and b αβ in Eq. (62) will correspond to in the case of solid particles interacting through lubrication. The corresponding thermal noise term will therefore read where a matrix W αβ of independent increments of the Wiener process has been introduced for every pair of solid interacting particles. The traceless symmetric part d W αβ , again, reads with the new lubrication-noise amplitudes A αβ and B αβ defined as The previous stochastic differential equations need to be integrated in time and the resulting force contributions updated (in addition to the forces coming from the SDPD solvent) in the rigid-body translational/rotational solid particle dynamics. Details on the derivation of the irreversible lubrication matrix and related GENERIC building block are provided in Appendix A.

Splitting scheme for the stochastic lubrication
Since the lubrication forces are singular (they diverge as 1/s or ln(1/s) for vanishing gaps, see Eq. 24), the integration problem becomes stiff and extremely small time steps are needed for tracking the dynamics of nearly-touching spheres. A very efficient semi-implicit splitting scheme similar to the one developed in [19,26,37] can be adapted here for taking into account the stochastic lubrication contribution in the many-particle system. The main idea is that for each solid particle α (instead of estimating the total stochastic lubrication force coming from all the neighboring solid particles and therefore integrating it explicitly on time), it is possible to split the problem in many simple two-body interaction problems and then solve them iteratively over all the pairs of solid particles. The efficiency of this strategy comes from two main ingredients: (i) the control on the overall threshold error and (ii) the possibility to solve the two-body problem analytically, in a way that no numerical matrix inversion is needed. The latter aspect makes the iteration procedure very fast [26].
More specifically, for every pair of solid particles (α, β) the velocity at the step n + 1 is given by where F det αβ = F lub,n α,β + F lub,t α,β is the deterministic part of the lubrication correction force. Given that we are splitting the time step in N sweep parts (so dt sweep = dt/N sweep ), the velocity V α , after a time dt sweep The two-body system to solve therefore reads (37) where V and V represent the old and new velocities respectively. The solution of this system is similar to the one discussed in our previous work [37], by replacing V αβ by where Given that linear momentum conservation must be held, individual velocities of the particles are finally calculated as The whole solution is obtained by iteratively solving the above equations for each interacting pair (α, β). To maintain the accuracy of the method and its robustness at all particle configurations, a pseudo time step dt sweep = dt/N sweep is used. Here, the number of subiteration N sweep is dynamically determined to ensure it is high enough to get convergence while simultaneously being sufficiently small to speed up the calculations. Beginning with a large default value for N sweep , at each time step say nth, two different sweeps are carried-out with N sweep = 2m and N sweep = 2m − 1. Then, the difference in the particles velocities are computed through an L 2 norm and compared against the value of a predefined tolerance . For further details, the reader is referred to [19,26].
Finally, the additional repulsive force acting between solid particles is introduced to prevent artificial particle overlap [2,37] F rep αβ = F rep τ e −τ s /(1 − e −τ s )e αβ where τ −1 determines the interaction range and F rep its magnitude. In order to model nearly hard-spheres, typically values of τ −1 = 0.001a and F rep = 21.15 have been used [37]. In this work τ −1 = 0.05 is chosen, corresponding to slightly longer repulsive forces. Since this represents a reversible contribution, it will not alter the GENERIC structure and therefore no particular treatment is needed here (beside enforcing its anti-symmetric form for particle swap which leads to local momentum conservation). All inter-particle interactions are implemented within the so-called Parallel Particle Mesh library (PPM) [31], a Fortran 90 software layer between the Message Passing Interface (MPI) and Client Applications for simulations of physical systems using Particle-Mesh methods with optimal scaling performance.

Specular reflection boundary conditions
Due to the presence of thermal noise, the SDPD particles can occasionally penetrate the solid regions (e.g. solid walls or regions occupied by solid suspended particles) within one time step. In order to prevent it, specular reflection boundary conditions are applied [28]: when a fluid particle i is penetrating into a solid boundary, its velocity projection normal to the solid boundary is inverted: where v i is the velocity before the application of the specular reflection, v i is the velocity after the reflection. In the case of rigid walls, we have being n wall the normal vector of the wall surface. Note that in the case of simple shear flow considered in Sect. 3, walls are moving tangentially to the fluid, so its velocity projected along its normal vector is identically zero.
In the case of a solid particle α, its normal velocity V α to the surface should be considered, in such a way that being n α the normal unit vector to the surface of the suspended solid particle α in the region where SDPD fluid particle i is penetrating. If the solid boundary is represented by a moving solid particle, both linear and angular momentum transferred to the solid sphere should be considered. The change of linear momentum produced by the specular reflection of the SDPD fluid particle i is given by while the change of angular momentum is being r i the position of the fluid particle i and R α the position of the center of the solid particle α. As a results of the overall possible reflections with SDPD fluid particles, the solid particle α will change its velocity V α according to where a summation over all penetrating fluid particles has been made. Similarly, the change on the angular velocity ω α of the solid particle should be corrected as where I α is its moment of inertia. In this way discrete conservation of linear and angular momentum is preserved.
In the case of penetration into the walls, there is no need to consider the transfer of linear moment to simulate the dynamics, since in the shear-flow simulations of Sect. 3 the wall velocity is fixed and prescribed. However, since the suspension viscosity will be inferred from tangential forces acting on the wall, an accurate force estimate is required. The extra force ΔF on the wall due to the specular reflection reads where the summation is over all penetrating fluid particles into the wall, dt is the time step, and Other possible choices for the velocity boundary conditions in SDPD include bounce-back and Maxwellian reflections. For a more detailed discussion on the effect of the several boundary conditions on the dynamics, the reader is referred, for example, to [28].

Numerical results
In this section we validate the single and multi-particle Brownian systems by checking their diffusional properties (e.g. velocity probability distribution function, particle mean square displacement and diffusion coefficient) as well as the rheological properties (i.e. the relative suspension viscosity) of the whole Brownian suspension. Results are compared with established data from the literature and discussed.

Isolated Brownian particle: diffusional properties
In this section, we validate the single solid particle diffusion. We consider a cubic simulation box of side L = 10 and a solid particle of radius a = 1 placed at the center. Physical parameters are: k B T = 10.0, fluid density ρ 0 = 1, fluid viscosity η s = 4. Solid particle mass is M = ρ 0 4 3 πa 3 = 4.188 and its thermal velocity is estimated as V th = √ k B T /M = 1.545. Resolution is taken as 5 SDPD particles per radius of the solid particle (see Fig. 1), so the mean fluid particle spacing is dx= 0.2a. SDPD particle mass is m = ρ 0 dx 3 = 0.008, leading to a SDPD fluid thermal velocity v th = √ k B T /m = 35.355. The speed of sound is chosen c s = 80 larger than V th to enforce a quasi-incompressible regime. With these parameters, the following time scales can be defined as follows: the sonic time scale t s = a/c s = 0.0125, the viscous time scale t ν = a 2 /ν = 0.25 and the diffusive time scale t d = a 2 /D 0 = 7.5, where D 0 = k B T /(6πη s a) is the solid particle diffusion coefficient. It is clear that the correct hierarchy of time scales is enforced, that is: t s t v t d , as physically required [40]. The total number of SDPD particles in the simulation box is N = 125, 000. Simulations are running until times t t d to ensure the correct achievement of the diffusive regime. The time step was dt = 5 × 10 −4 . Figure 2 shows the velocity statistics and mean square displacement of the solid particle defined as with the average being performed over multiple independent realizations. From Fig. 2 (left) it can be seen that the velocity probability distribution function match well the Maxwell-Boltzmann distribution in a given direction (x) at the prescribed temperature This is accurately reproduced both, for the SDPD solvent particles as well as for the suspended solid sphere. Figure 2 (right) shows the MSD with the two characteristic regimes: ballistic (∼ t 2 ) and diffusive (∼ t) at time scales, respectively, smaller and larger than the viscous time t ν . An estimate of the diffusion coefficient D for the particle can be inferred from the MSD, leading to a computed value D = 0.1 ± 0.01. Since we are using periodic boundary conditions, this value should be compared with the corrected diffusion coefficient D = D 0 λ obtained by taking into account the effect of the periodic images of the solid particles. Here D 0 is the limiting value obtained from Stokes-Einstein theory, whereas the factor λ reads [30]: where φ = 4πa 3 /(3L 3 ) is the effective spheres concentration in a periodic domain. This is also in good agreement with the correction proposed in [45] λ(a, L) = (1 − 2.84a/L). Considering that in the present simulation φ = 0.0041888, this gives an analytical estimate of the effective diffusion coefficient D = 0.09555 in good agreement with the numerical value.

Brownian multi-particle system: rheology
In this section we show results of the suspension viscosity for a mono-dispersed Brownian system and compare them with established data in a range of flow conditions and particle concentrations. To start, we define the concentration φ in the multi-particle system as where a is the solid particle radius and N c is the total number of suspended solid particles.
The system in this case is confined between two walls placed at a distance L z = 10 apart and moving in opposite directions ±V wall (see Fig. 3). In this way a uniform linear velocity profile is generated within the suspension, characterized by a constant shear rateγ = 2V wall /L z . The shear rate provides a typical time scale of the flow denoted as t f = 1/γ . In inertia-less Brownian systems, the other relevant physical time scale is the already defined viscous time t d = a 2 /D 0 associated to the diffusion of the solid particle. These two scales can be coupled together to define the so-called Péclet number At small temperatures or large particle sizes, the diffusion is an extremely slow process, i.e. t d t f → Pe 1 and the dynamics is governed by t f . On opposite, for Pe 1 the dynamics is governed by the particle diffusion. From the competition between these two time scales, a so-called "shear-thinning" rheological behavior for the suspension is typically observed where its viscosity generally shows a reduction trend for increasing Pe. At very large Pe, the viscosity shows a plateau corresponding to the deterministic limit, i.e. non-Brownian systems. However, depending on the specific short-range inter-particle forces, opposite shear-thickening can eventually sets in at large Pe. The effect is increased at large particle concentrations [13,21]. In the next section we report the results for the suspension viscosity as a function of Pe and for two typical values of φ = 0.1 and φ = 0.3, respectively, in the semi-dilute and concentrated regimes. They correspond to systems with N c = 24 and N c = 72 suspended solid particles. In this section identical parameters to Sect. 3.1 are used. The Péclet number is modified by varying V wall in the range ∼ 0.066 − 66.

Brownian suspension viscosity
In order to extract the viscosity from the simulations, we compute the wall shear stress from the total tangential force F wall exerted by the fluid and suspended solid particles on the walls. The component σ xz is related to the x component of F wall as σ xz = F wall x /A where A is the wall area. The total suspension viscosity for the given solid volume fraction φ finally yields The dimensionless relative suspension viscosity will be reported which is defined as η r (φ) = η(φ)/η s . Figure 4 shows the relative suspension viscosity as a function of Pe for φ = 0.1, 0.3. At low solid volume fraction (φ = 0.1, black points) no shear-thinning is observed with a constant value of η r ≈ 1.3 in excellent agreement with the theoretical result η r = 1 + 2.5φ + 6.2φ 2 calculated by Batchelor [1] in the semi-dilute regime and previous numerical results [13] (dotted black line).
At larger solid volume fraction (φ = 0.3, red points) a marked shear-thinning behavior is observed. By fitting the relative suspension viscosity data with the following function the estimated values at low and large shear rates are, respectively, η 0 r = 3.83 and η ∞ r = 2.51 which are in line with the values at φ = 0.3 obtained from the experimental correlations of de Kruif et al. [6] using a number of different Brownian samples, i.e. η 0 r = 3.24 ± 0.46 and η ∞ r = 2.77 ± 0.45. Note that a mild shear-thickening sets-in at Pe ≈ 100 in agreement with previous simulations results [13]. The present results show good agreement with experimental data and previous numerical calculations at the two given solid volume fractions and provide an indication on the accuracy of the present SDPD model for the simulation of Brownian multi-particle systems. More analysis is required to determine the effect of the rheometer gap-to-particle size ratio on the results as well as the repulsive parameter (F 0 ) and range of the repulsive forces (τ ), which is the matter of current numerical analysis. framework of GENERIC and therefore it complies with the First and Second Laws of Thermodynamics and Fluctuation-Dissipation Theorem at the discrete level. More specifically, the model splits the description of hydrodynamic interactions between suspended particles in two parts: a far-field stochastic hydrodynamic contribution is taken into account by an explicit fluid model, namely the Smoothed Dissipative Particle Dynamics method which is a thermodynamically-consistent discretization of the fluctuating Navier-Stokes equations. Short-range hydrodynamics interactions are modelled by a generalization of the Lubrication Dynamics method and solved by a semi-implicit splitting procedure. The lubrication model is generalized to the stochastic case by incorporating the effect of a lubricationinduced thermal noise into the dynamics of closely interacting solid particles in a way that is GENERIC-compliant. Preliminary validations include Brownian diffusion as well as rheological response of a multi-particle system in the semi-dilute and concentrated regimes, showing good results.
Diffusivity analysis as well as rheological response in more concentrated systems is currently under work and should allow to validate the model in the highly-concentrated limit. Unlike other techniques such as Brownian Dynamics or Stokesian Dynamics, one of the benefits of the present SDPD-based model for Brownian suspensions is that, for example, the non-Newtonian properties of the liquid phase can be straightforwardly included. Viscoelasticity has been recently included in a GENERIC-compliant multi-particle model in the deterministic limit [42], and can be extended in the future to the stochastic regime of fluctuating viscoelasticity in a similar way to that proposed here.
where the terms of spontaneous momentum transfer and heat fluxes between pairs of particles are postulated as follows being W αβ a matrix of independent Wiener processes associated to particles α and β. The traceless symmetric part d W αβ reads Note the similarity between these noise terms and the postulated ones for the model of liquid, Eqs. (11) and (12). The irreversible matrix M can be calculated via the dyadic product of the noise terms, Eq. (5), yielding The calculation of the irreversible matrix M blocks leads to next equations More details can be found in [12]. With this matrix M, the irreversible dynamics of the system can be finally calculated leading to the following equations of motion d P α | irr = − β a αβ V αβ + a αβ 1 − 2 D + b αβ e αβ · V αβ e αβ dt + β dP αβ .