1 Introduction

The main purpose of this paper is to provide an overview of the smoothed dissipative particle dynamics (SDPD) method. The numerical technique was introduced in 2003 by Español and Revenga^[1] to remedy some deficiencies of the classical dissipative particle dynamics (DPD) ^[2-4]. Since then, it has been applied successfully to a wide range of problems in microfluidics^[5], nanofluidics^[6], colloidal suspensions ^[7-8], blood ^[9-10], tethered DNA^[11], and dilute polymeric solutions ^[12-14]. It has also been used for the simulation of fluid mixtures ^[15-18], mesoscopic viscoelastic flows^[19], and multiscale coupling strategies ^{[18, 20-21]}.

Rather than providing a complete overview on the state-of-the-art applications and challenges of SDPD (for a very complete recent review of DPD and related formulations, see Ref. [4]), we will focus here on a didactic overview on both fundamental and technical aspects of the methodology, which, in our opinion, have generated some misunderstanding in the past years or have not been fully recognized.

In order to make the presentation amenable to non-technical users, we propose a format for the paper in terms of simple questions/answers concerning SDPD. This should help to clarify some technical issues related to the original motivations behind its early development and to offer a simple conceptual framework for the computational scientists who would like to approach this technique for the first time. The answers represent the state-of-the-art in terms of the physical interpretation of the SDPD methodology as a "mesoscopic approach", and reflect the personal opinion of the authors. They are not meant to be definitive but, instead, the goal here is to stimulate further debate among the community as well as to explore, in more quantitative terms, the main SDPD advantages and limitations, also in relation to established methods.

2 What is SDPD?

SDPD is a Lagrangian particle method for the numerical solution of thermal flow problems able to address both, the macroscale and the mesoscale. Roughly, the mesoscale involves time scales from 1 ns to 10⁶ ns and space scales from 10 nm to 10⁴ nm. In SDPD, a set of fluid particles i=1, 2, ..., N are distributed homogeneously over a physical domain and move according to forces estimated from their local neighborhood. The following set of stochastic differential equations for the particle positions, velocities, and entropies define the model, and are solved numerically:

(1)

Instead of the entropy, we may use the internal energy as independent variable, evolving according to

(2)

Here, the mass density is calculated by , where is the volume associated with the particle i, and is the number density, where W_ij=W(r_ij, r_c) is a kernel function with compact support r_c, and

e_ij=r_ij/r_ij is the unit vector of the joining particles i and j. v_ij=v_i-v_j is the particle velocity difference, and T_ij=T_i-T_j is the particle temperature difference. and are random exchanges of momentum and energy per unit time between particles i and j. These random force and a random heat transfer, respectively, are specified below.

To close the system of equations (1) and (2), any arbitrary equation of state for P_i=P(ρ_i, T_i) can be specified. The above dynamic equations have the general equation for non-equilibrium reversible-irreversible coupling (GENERIC) structure^[22]. This implies that they are thermodynamically consistent, i.e., they respect the first and second laws of thermodynamics, and the thermal fluctuations lead to the correct thermal fluctuations governed by the Einstein formula for the equilibrium probability distribution. From the momentum equation, one can notice that the total force acting on a particle i is the sum of pairwise contributions of different nature: (ⅰ) a conservative force F_ij^C; (ⅱ) a dissipative force F_ij^D; and (ⅲ) a random force F_ij^R.

To interpret these equations, let us neglect the stochastic contributions in Eq. (2). The resulting equations represent a Lagrangian discretization of the non-isothermal Navier-Stokes equations (NSEs) using the so-called smoothed particle hydrodynamic (SPH) interpolation ^[23-28]. In a nutshell, the conservative force is a specific SPH discretization of the pressure gradient in the NSEs. The dissipative force is a specific SPH discretization of the viscous terms in the NSEs, where η is the fluid dynamic viscosity. The last equation represents the SPH discretization of the entropy equation, where κ is the thermal conductivity, and ϕ_i is the viscous dissipation function calculated on the particle i. In the energy equation (2), the first term is just the reversible work due to compression of the fluid, while the remaining terms have the same meaning as those in the entropy equation.

It is worth mentioning that the requirement that the equations should have the GENERIC structure imposes a very specific structure for F_i^C, F_i^D, and F_i^R. In this way, while the equations may be understood as a numerical discretization of the NSEs, it is a very special one that respects at the discrete level the first and second laws and the correct equilibrium distribution. In fact, even though many different ways of discretizing continuum equations exist, which converge to the thermodynamic compliant continuum equations, potential problems might arise in simulations, especially at the mesoscale scale, if the algorithm itself is not thermodynamic compliant.

The explicit form of the random terms is given in Ref. [1]. Note that an equivalent and simpler form (with less random number generations) has been given recently ^[29-30], and reads as follows:

(3)

where

V_ij(t)=-V_ji(t) are independent white noises with covariances given by

(4)

where the averages are over realizations of the Wiener process. Finally, to satisfy the fluctuation-dissipation theorem (FDT) that guarantees the correct equilibrium distribution, the amplitudes of the noise in Eq. (3) are given by

(5)

Note that the random force satisfies Newton's third law, i.e.,

and the heat transfer complies with energy conservation

The SDPD method represents a proper stochastic extension of the SPH method for the Lagrangian discretization of NSEs, with thermal fluctuations introduced in a thermodynamically consistent way. SDPD is therefore a stochastic Lagrangian discrete representation of fluctuating hydrodynamics^[31].

In the context of particle methods, another stochastic model aiming at the representation of fluctuating fluids is the DPD method ^[32-39]. The SDPD model developed in Ref. [1] was motivated by some technical drawbacks presented in the original DPD approach which are discussed in detail below.

3 Why do some scientists have trouble using classical DPD?

Classical DPD is a stochastic particle method, which has a structure very similar to SDPD with the force acting between particles having soft conservative, dissipative, and random contributions. The functional forms for these forces in classical DPD are, however, somewhat simpler than in SDPD, and they read

(6)

where a_ij is the strength of the repulsive potential, γ is a friction coefficient, and σ is the strength of the thermal noise.

is a weighting function for the conservative force, whereas ω^D(r_ij) and ω^R(r_ij) are two additional weighting functions. In order for the noise to satisfy exactly the FDT, the input parameters need to satisfy σ²=2γk_BT, and, at the same time, the weight functions need to be chosen such that

No energy or entropy equation is considered in the original DPD model.

This set of equations are manifestly simpler than those for SDPD in Eq. (1). The price for the simplicity of DPD is, however, that there is no direct connection between the model parameters and the physical parameters of the liquid system one tries to simulate. In particular, the following technical drawbacks are presented in the conventional DPD^{[2, 3, 40]}.

Isothermal model The DPD model can be understood as a thermostat, and can only deal with isothermal situations. Energy transport is precluded in the model.

Transport coefficients The expressions for the DPD transport coefficients, e.g., viscosity, can be obtained from kinetic theory ^[40-41]. A common expression considered by DPD users to "estimate" the fluid viscosity is

(7)

This expression is useful but of limited validity. In fact, it has been specifically derived from kinetic theory by assuming zero contribution from the conservative forces. It is therefore accurate in a so-called gas-like regime, where the particles are very weakly interacting with each other. When the particles interact significantly via the conservative forces, i.e., at high shear flow, low temperatures or in complex geometries, where compressible effects are more important, Eq. (7) gives a poor estimate of the viscous behavior of the model DPD fluid. In addition, recent results show that the measured transport coefficients under viscometric conditions are also highly sensitive to parametrization, integration method, and thermostating schemes^[42]. As a result, it is common among DPD users to map and calibrate the liquid properties under different choices of input parameters in a way not always systematic and/or unique ^[43-44].

Equation of state (EOS) for the pressure The EOS is fixed in classical DPD by the particular form of the conservative forces, and it reads

where α ≈ 0.101^[40]. This limits the range of the validity of DPD that cannot deal with the arbitrary thermodynamic behavior of real fluids.

Size of DPD particle It is difficult to specify in advance the scale at which a DPD simulation operates. In particular, there is no parameter in the model that sets the physical scale of the particle. This problem is crucial, for instance, in the cases of suspended colloidal particles or in microfluidics applications where the physical dimensions of the external objects determine whether and, more importantly, to what extent thermal fluctuations come into play.

In relation to this point, it should be remarked that no formal link between a DPD particle, understood as a coarse-grained entity, and the underlying atomistic dynamics has been established so far^[4]. The exact application of the Mori-Zwanzig projection formalism for the construction of coarse-graining (CG) dynamics has been limited so far to systems where the atoms are bonded together, e.g., membranes and molecules, but cannot be easily extended to unbonded atomistic systems, e.g., a liquid, where the very definition of CG variables is a challenging issue.

Resolution in DPD Because there is no specific parameter in DPD associated with the size of a particle, it is very difficult to check the effect of particle resolution on a given problem in DPD. To illustrate this point, consider the case of a flow around an obstacle, e.g., a solid sphere of physical size R. Ideally, for predictive calculations, one would like to have output simulation results, e.g., drag, which are independent of the number of particles N used in the simulation. However, if one changes the resolution in DPD, i.e., moving from r_c to r_c/2, i.e., N → 2^D N, where D is the physical dimension, the transport coefficients, e.g., viscosity, will change in a very complex way (see Eq. (7)), so that one does not have the same fluid anymore.

There have been attempts to address the issue of resolution, i.e., linking the values of the DPD model parameters with the number of used dissipative particles based on the notion that a dissipative particle represents N_m atoms packed into it^[45]. As a consequence of the increased effective repulsive forces, it was shown that there had already been disappointingly low N_m ( ≈ 10) atoms per DPD particle, which resulted in an artificial crystallization of the DPD fluid, making an atomistically-motivated choice of the interaction parameters in DPD useless for practical flow calculations.

The problem has been bypassed in Ref. [46], where it was shown that by assuming a specific scaling of DPD model parameters with the physical scale, the balance between conservative and random terms can be restored back to "safe" fluid levels. In that mapping, however, the link to the atomistic systems (and therefore to the real size of a DPD fluid element) was completely lost. In fact, it was shown in Ref. [46] that the conventional DPD dynamics could be made scale-invariant. In addition, as explicitly stated by the authors in the specific scaling proposed works, it is only for equilibrium phenomena, e.g., static properties and phase diagrams, and not for dynamic processes, e.g., transport coefficients.

In a nutshell, it is very difficult in classical DPD to assess the sensitivity of the simulation output to the specific choice of the adopted model parameters and particle resolution. This feature does not prevent DPD from nicely reproducing different flow behaviors, but it makes it almost incapable to predict the physical results starting from "a priori incomplete knowledge of the problem".

Since the original formulation of DPD, many different models have been devised in order to ameliorate the problems in the original model. In order to have arbitrary equations of state, the manybody DPD (MDPD) was introduced in Ref. [47], where the conservative forces were derived from a free energy function. In order to deal with non-isothermal situations, energy conserving DPD (EDPD) was introduced by including an internal energy for the dissipative particles ^[48-49]. In order to have better representations of the friction forces and angular momentum conservation, the fluid particle model (FPM) was introduced in Refs. [50] and [51]. SDPD can be regarded as the final outcome of these endeavors, grouping all the partial solutions to the above problems.

4 What is the size of an SDPD particle? Is the SDPD method a mesoscopic or a macroscopic method?

The physical size of an SDPD particle is determined by the user. The criteria are to use a sufficiently large number of particles to resolve the physical length scales of the system under study, in such a way that by increasing further the number of the particles, the results do not change. In fact, in SPH, as in every other computational fluid dynamics (CFD) numerical grid-based method solving continuum partial differential equations, e.g., finite difference method, finite volume method, and finite element method, the control fluid element represents a purely discrete volume, where mass/momentum conservation is locally enforced and it needs to be sufficiently small (compared with the geometrical length of the problem) to guarantee the numerical convergence.

In SDPD, the level of fluctuations associated with every particle is not-scale invariant, and this is consistent with the fact that the particle represents a physical thermodynamic system made by a finite (not infinite!) number N_c of constituents (molecules, atoms, etc.). According to statistical mechanics, the size of fluctuations should scale as , vanishing only in the continuum limit, i.e., when the size of the fluid element largely oversizes that of the contained constituents. Quite remarkably, the thermodynamic consistency of SDPD implies that the particles are subject to the fluctuations with the correct physical scaling. In fact, velocity fluctuations are distributed according to a Maxwell-Boltzmann function and their magnitude scales in such a way that

(8)

where D is the dimensionality of the problem. Now, because k_B, T, and ρ₀ are physical properties specified once for all at the beginning, the velocity fluctuations depend uniquely on the fluid particle volume . In other words, the size of the SDPD particle uniquely determines the level of fluctuations.

To illustrate this feature, let us consider the example of a solid sphere with radius R suspended in a liquid. In this case, to get accurate results, we would need to resolve the flow around the sphere with elements, say, 10 times smaller (for the sake of numerical convergence) than R. If we have a macroscopic sphere, i.e., R ≈ 1 m, then r_c^SDPD=0.1 m. This would give a thermodynamic volume ≈ 10^-3 m³. On this scale (for water with the density ρ₀=1 000 kg/m³ at room temperature T=298 K), the fluctuations on the SDPD particles will be uniquely determined to be 〈v²〉^1/2 ≈ 6 × 10^-11 m/s. Under common macroscopic flow problems, this fluctuating contribution (and corresponding random terms in Eq. (1)) can realistically be neglected, and therefore we end up with a deterministic SPH discretization of NSEs.

However, in the case of a flow around a nanosphere with the radius R=100 nm, r_c^SDPD=10 nm=10^-8 m, and ≈ 10^-24 m³. As a result, in this case (for the same fluid), 〈v²〉^1/2 ≈ 2 m/s, which is dominant under microfluidics conditions. Fluctuations lead to correct non-zero diffusional dynamics of the nanosphere.

In summary, whereas the level of thermal fluctuations in DPD is arbitrary (in fact, with certain prescribed scaling^[46], the method is scale-invariant), in SDPD, it is exactly prescribed by the scale of the problem. In the limit of large systems, thermal noise can be neglected, i.e.,

(in Eq. (1), we recover the deterministic SPH equations). On the contrary, under mesoscopic conditions, we have

the Brownian diffusion is critical, and we need the full stochastic SDPD model. This scale property is a built-in feature of SDPD, and no tuning is necessary.

Therefore, SDPD can be either a mesoscopic or a macroscopic flow solver, depending on the specific scale of the problem under study. For macroscopic flows, SDPD is just a particular, thermodynamically consistent, version of SPH.

5 Is there any link in the SDPD with the underlying atomistic system?

No link has been established so far. Equation (1) is derived irrespectively of the details of the underlying atomistic fluid, whose identity is kept only through the values of the viscosity and thermal conductivity and the equations of state. In this sense, the SDPD model is mesoscopic. However, it is understood as a top-down approach with a clear link to macroscopic continuum NSEs and correct incorporation of thermal noise. Moreover, there have been claims that DPD is obtained from the first principles in a bottom-up approach ^[52-53]. There is, however, a misleading concept on DPD here, at least with regard to its application to "fluids". Neither SDPD nor DPD has a direct formal link to the underlying microscopic, e.g., atomistic system. In fact, as discussed in Section 3, there have been no satisfactory derivations of DPD for fluid systems from the Mori-Zwanzig formalism. The problem is mainly related to the difficulty in the definition of a coarse-grained particle for a system made by unbonded diffusing atoms, i.e., a liquid. The best one can do, so far, is to aim at a correct description of the hydrodynamic regime with a properly introduced thermal noise.

It should be noted, however, that

(ⅰ) DPD remains a useful and accurate CG method for systems with bonded atoms, where the link to the atomistic level has been made ^[54-55].

(ⅱ) In the case of fluids, DPD remains a simple stochastic method which delivers "hydrodynamic behaviors" on larger length-time scales^[56] (albeit with the difficulties mentioned in Section 3).

(ⅲ) In both DPD and SDPD, one cannot attribute the physical meaning to the spatial processes such as correlations occurring below the SDPD or DPD particle scale r_c. However, it might still be possible to use/tune (unphysical) numerical statistics occurring on small space-time scales, i.e., smaller than r_c, to model different physical processes displaying similar statistical behaviors. This follows the spirit of the so-called implicit modeling^[57], where, for example, turbulent dissipation models are constructed, which start from the discrete properties of the chosen numerical scheme rather than a predefined continuum model to be discretized.

6 Does SDPD conserve total energy and linear/angular momentum?

Linear momentum in SDPD is conserved locally and globally. This comes from the fact that the pairwise forces acting between particles are anti-symmetric by the swapping particle index, i.e.,

(9)

This is consistent with Newton's third law that ensures that global momentum is also conserved, i.e.,

(10)

Total energy is conserved in SDPD at the discrete level. This follows directly from the GENERIC structure of the equations, but one can show directly from Eqs. (1) and (2) that

(11)

where we use the symmetric properties of the summation under the index swapping and the discrete momentum equation (1) to arrive at .

Angular momentum is not conserved by the original SDPD formulation given in Eq. (1). This is due to the term proportional to v_ij in the friction force which introduces tangential forces between particles. As the SDPD equations are a discretization of NSEs which conserve angular momentum, angular momentum is conserved in SDPD in the limit of high resolution. However, lack of conservation at the finite resolution discrete level is a severe problem in some applications including particle suspensions, and might lead to serious numerical artifact if it is not exactly enforced^[58].

Two routes have been considered to restore angular momentum conservation in SDPD: (Ⅰ) to change the viscous force; (Ⅱ) to introduce an additional spin variable. Let us discuss these two possibilities separately.

(Ⅰ) First of all, recall that SDPD is directly linked to the specific discretizations of the NSEs which are not unique. Therefore, any specific viscous force representing a possible discretization of NSEs and consistent with GENERIC is acceptable.

To illustrate this point, let us consider the general structure of the irreversible (viscous) part of the momentum equations as specified by GENERIC. This reads

(12)

where a_ij and b_ij are positive interparticle strengths. Now, the classical choice made in the original SDPD is

(13)

which leads to Eq. (1). However, other choices are possible. In fact, if one considers, for example,

(14)

where D is the problem dimensionality, this will lead to the angular momentum preserving SDPD formulation proposed by Bian et al.^[7] and Hu and Adams^[59], which is linked to a particular SPH discretization of the viscous dissipation in the NSEs ^[60-61]. In fact, with the choice made in Eq. (14), interparticle forces are uniquely directed along the center-to-center vector joining two particles, i.e., F_ij^D ∝ e_ij, which automatically ensures the conservation of angular momentum. The selection of Eq. (14) represents a common SPH representation of the viscosity term ^[59-60]. This expression has been mostly used in connection to the SPH/SDPD simulations of suspended particles ^{[7, 62-63]}.

It should be remarked that, for thermodynamic consistency, GENERIC also requires that the stochastic terms need to be modified in order to satisfy exactly the FDT with this new choice of the discrete irreversible particle dynamics. In other words, in SDPD, we have a variety of options that can be chosen for the interparticle forces, as long as they satisfy the GENERIC structure in Eq. (12).

(Ⅱ) A second possibility is represented by the introduction of an additional spin variable ω_i associated with each particle. In Ref. [64], it was shown that by considering an additional evolution equation for ω_i, which can be cast into the GENERIC framework, the overall angular momentum can be conserved by adopting the original SDPD formulation given in Eq. (1), i.e., preserving the term proportional to v_ij. The deterministic part in this case corresponds to a SPH discretization of a modified set of NSEs for fluids with spins^[65].

In conclusion, in SDPD, the total linear momentum, angular momentum, and energy are exactly conserved by the specific discrete set of equations. This is a remarkable property of the method because the specific discretization adopted satisfies the same momentum/energy conservation as the original continuum PDEs.

7 Does SDPD have an H-theorem?

The GENERIC structure of the SDPD equations automatically ensures that the second law is satisfied. At a deep level, entropy increase is structurally linked to the fluctuation-dissipation theorem. Indeed, the covariance of the noise (that gives the amplitude of the noise) is given by the dissipative matrix, which is, therefore, a positive definite matrix. As emphasized in the GENERIC framework, this latter property ensures the entropy increase in the evolution of the system. The total entropy is the sum of individual particle entropies If we neglect fluctuations and take the time derivative of S, we have

(15)

where the dissipation function is estimated on the particle i and the positivity of -W'_ij as follows:

(16)

This latter property is a constrain on the shape of the weight function that is satisfied by usual bell-shaped weight functions.

Therefore, SDPD has an H-theorem valid at the discrete level.

8 How do SDPD results depend on the particle number we choose?

As mentioned in Section 3, in SDPD simulations, the convergence criteria for the numerical results can be directly assessed. In other words, s in any CFD method, one can fix once for all the physical quantities characterizing a given problem, e.g., liquid speed of sound c_s, molecular viscosity η, thermal conductivity κ, temperature T, and physical geometry size R, and can change the number of SDPD particles N and their size r_c without having to worry about varying liquid properties or flow regimes.

It remains to understand what is the level of accuracy in SDPD, i.e., what is the typical resolution to be used for converged results. As SDPD is a stochastic version of SPH, it is helpful to resort to its classical error analysis. In SDPD (and SPH), any function f(r) defined over a domain Ω is represented through a two-step approximation process. In the first instance, as a result of the smoothing process with the interpolation kernel W, f can be approximated by

(17)

The so-called mollification error is the consequence of replacing the exact Dirac δ-function in an unapproximated convolution with a finite-size kernel W(r).

The second approximation comes from the quadrature of Eq. (17), i.e., by moving from a continuum space to a discrete space made by a random set of N points, i.e., particles, separated by a mean distance Δ, and it reads

(18)

where the usual replacement is done, i.e., , and , in which is the volume associated with the particle i introduced in Section 1.

In conclusion, the overall numerical error in SDPD (SPH) is the sum of two terms, i.e.,

(19)

which is a mollification error and a quadrature error^[66], where m and n represent the corresponding orders of convergence. In Eq. (19), r_c is made dimensionless with a specific problem length R. The first term, therefore, tells us how small the particle size must be with respect to the physical size of the geometry under consideration. The second term quantifies the effect of the average number of points used for the quadratures, i.e., the local number of neighbors N_neigh per particle. It can be shown that, if W(r) is a normalized and radially symmetric kernel, the exponent n=2. Therefore, in connection with the mollification error, one has the second-order spatial convergence, which is similar to standard grid-based techniques^[67]. It must be noticed, however, that the presence of the quadrature error introduces a residual which, if not properly taken into account, prevents the mathematical convergence. In fact, one would ideally require N_neigh to increase at a rate such that the second term is always smaller than or comparable to the first one. This might be very impractical to do. It should be noticed instead that the analytical estimates of the exponent m for irregular particle maps predict a value depending on the adopted level or regularity of the kernel W^[67]. For the classical quintic spline kernel W (4-times differentiable) normally adopted in SPH/SDPD, we have that m=3, which means that the truncation errors go quite rapidly to zero and, in most cases, the regime, where the residuals are dominant, might not be even accessible. In Ref. [68], for example, it was shown that for viscous-dominated flow problems, e.g., microflow conditions, by keeping fixed r_c/Δ, the resolution-independent results within few percents can be obtained, which is very satisfactory for practical engineering calculations.

Strictly speaking, the present discussion is valid only for SPH, i.e., SDPD without noise terms. The truncation errors in SPH are highly dependent on the particle map, being minimal for a perfect lattice (m=β+2, where β is a property of the kernel W, being the highest integer such that the βth derivative and all lower derivatives are zero at the edges of the compact support) and maximal for a completely random particle distribution (m=0.5). In SPH, the particles under flow are disordered but not randomly distributed. Therefore, its accuracy is expected to be between those two limits.

Moreover, it is remarkable that by relaxing the SPH particle map to that one commonly observed in SDPD, i.e., a liquid-like radial distribution function, a very high-order integration error exponent (m ≈ 8, when a C⁶ Wendland kernel is used) can be achieved in SPH, which is similar to the results obtained only with the particles located on a lattice^[69]. In this sense, thanks to the specific particle distribution enforced in SDPD by the stochastic forces, its accuracy is expected to be even larger than that in SPH where this homogenization mechanism is absent. This could allow to choose an even smaller number of neighbors N_neigh in SDPD compared to standard deterministic SPH formulations. More research in the direction of "a bit of noise increases numerical accuracy in SPH" would be welcome.

9 Is SDPD slower than the classical DPD?

There is a common belief that SDPD is a slower method when compared, e.g., with classical DPD. This judgment is motivated by two main arguments: (ⅰ) a more complex structure of the equations (i.e., interparticle forces); and (ⅱ) the fact that, in DPD, a small local number density is typically considered (n=N/V≤ 4), leading to an average number of neighbors N_neigh = 4πn/3 ≈ 17. In contrast, the modus operandi of SDPD requires a higher number of neighbors, typically on the order of 30 in 2D and 100 in 3D (for a typical choice r_c=3 Δ, where Δ is the average interparticle spacing). Let us discuss these two aspects separately.

(ⅰ) It is true that the algorithmic complexity of the SDPD equations is somewhat higher than in DPD. Note that for short-range interacting particle systems, the use of linked-list-cell routines reduces the total number of pairwise interactions to ~ (NN_neigh)^[70]. The algorithmic structure of the neighbors searching procedure is the same in the two methods, and thus the differences lie ultimately in the force calculation only. The complexity details of each force in the athermal DPD vs SDPD are listed in Table 1, where we count +, -, ×, and ÷ each as one floating-point operation and power as three floating-point operations^[70]. In relation to the random number generator, we assume that, in both methods, one single ξ_ij normally distributes with zero mean and unit variance counts as, say, k operations.

Note that in SDPD, the density can be absorbed in a previous pair loop, and therefore it does not add any extra complexity. Similarly, in SDPD, all other terms depending on single particle index i (d_i², P_i, ...) can be pre-estimated in a single particle loop before entering the pair-loop with negligible extra costs. The maximum load in SDPD comes therefore from the need to generate a tensorial Wiener process (6 components in 3D due to symmetry considerations, instead of 3 of DPD) in the random force and from an extra matrix-vector multiplication. Note that the new implementation ^[29-30] of the random force in Eq. (3) requires only four random numbers and the use of uniform random numbers instead of Gaussians further reduces the cost^[29]. In any case, it can be seen that, independently on the used algorithmic complexity k of the normally deviated random number generator, the overall ratio per time step between the two methods is less than 2.

(ⅱ) With regard to the average number of neighbors N_neigh, the algorithmic complexity factor between DPD and SDPD should be on the order of (N_neigh^SDPD/N_neigh^SDPD) which is typically on the order of 3 (2D) and 6 (3D).

Note that the request of a large N_neigh^SDPD in SDPD is dictated by convergence requirements. That is, for these large values of N_neigh^SDPD, measured local properties (velocity, pressure profiles, etc.) or integral quantities (drags, flow rates, etc.) do not depend on the adopted resolution. Unlike DPD, it is precisely that this convergence property of SDPD makes the model predictive, as the numerical output of the simulation is not particle-number dependent. There is no free lunch for this property, and the slowing factor of SDPD mentioned above is the price to pay for it.

On the other hand, we should stress here that this feature represents a "modus operandi" in SDPD rather than a method drawback. In fact, no one prohibits to use SDPD in a so-called "unresolved mode" with N_neigh^SDPD ≈ N_neigh^DPD. Similar to DPD, some kind of effective coarse-grained hydrodynamics will emerge but a calibration of the transport coefficients becomes necessary as in the original DPD model.

In a nutshell, SDPD might be slower than DPD when used in "resolved mode". That is the price to pay for having numerical convergence and accurate control on every transport coefficients. If one is interested in qualitative hydrodynamic aspects, one may relax this requirement, and similar CPU time is obtained in both methods.

10 How is the fluid complexity incorporated in SDPD?

Applications of mesoscopic methods in microfluidics often require the incorporation of new complex physics like, for example, non-Newtonian behavior. The standard SDPD equation (1) describes correctly a fluctuating Newtonian liquid, and thus the question is "how can we generalize the method to more complex fluids?". These systems are traditionally modeled in conventional DPD based on a mechanistic perspective. In fact, new forces (elastic, rigid, etc.) can be added between particles in order to describe long polymers ^[71-74], rigid or deformable particles in suspension^[75], porous media^[76], droplets^[77], cells ^[78-80], or even more complex thixotropic materials^[81] (see Fig. 1). This can be also viewed as a bottom-up approach where complexity is fed at the mesoscopic particle level and re-emerges as complex hydrodynamic behavior on the continuum scale. We denote this approach as mesoscopic viscoelasticity.

This way to simulate complex fluids has several benefits. First of all, it allows to incorporate physics based on a physically-grounded force-based approach. Secondly, but not less important, it allows to bypass the need of constitutive partial differential equations, which, for some complex fluids/materials, might not be theoretically derivable or of questionable accuracy. The drawback is that it might not be always easy to tune the mesoscopic physics (interparticle forces) to deliver a desired (e.g., experimental) target flow behavior on the continuum level.

Exactly, the same approach can be adopted to simulate suspended structures in SDPD, albeit with a better prescription of the thermal and physical properties of the solvent medium. As a matter of fact, SDPD has been successfully applied to model polymer molecules in suspensions ^{[11-12, 14]}, cells ^{[10, 82-85]}, mesoscopic multiphase flows ^[86-87], and colloidal or non-colloidal particles ^{[7, 62-63]}. We stress that, in SDPD simulations, the output flow behavior of the bulk systems is not tunable. To illustrate this point, let us consider the rheology of non-colloidal¹ rigid particle suspension studied in Ref. [63]. Output results depend uniquely on the physical model adopted for the solvent medium (a Newtonian inertia-less fluid in that case) and not on hidden details of the numerical implementation. As a matter of fact, converged rheological results for the suspension viscosity have been achieved in excellent agreement with previous results obtained by alternative techniques (i.e., Stokesian dynamics)^[89] and experiments^[88].

¹The term "non-colloidal" corresponds to a regime of flow (i.e., high Peclet number), where Brownian motion is irrelevant and a deterministic SPH description of the solvent can be used. This is the typical case in experiments with 40 μm beads suspended in glycerol mixtures under moderate shear flow^[88].

The fact that the obtained flow behavior is "numerically converged" and it uniquely relies on the physical model adopted with no tuning parameters (such as effective particle hydrodynamic radius, effective rigid-core radius, and effective bulk viscosity/solid volume fraction), precisely allows the technique to be predictive, as opposite to descriptive, i.e., being able to reproduce complex fluid rheology. This is one of major advantages in SDPD.

In the case of particle suspension with complex silicon oils (weakly shear-thinning fluids), it has recently been shown that, by exact parametrization of the rheological properties of the solvent, simulations can deliver flow behavior for the bulk suspension in quantitative agreement with experiments without any tuning^[90].

In the cases mentioned above, the accurate description of hydrodynamic interactions between suspended structures allowed for good results against experiments. It should be kept in mind that, in some cases, the exact nature of the hydrodynamic interaction might not be the essential aspect to capture. For example, in highly concentrated particle suspensions, short-range friction forces might be more important, or when modeling some melts, gels, or other thixotropic materials, the non-hydrodynamic nature of the interactions can be the relevant one. In all those cases, accuracy of SDPD (as well as of any other hydrodynamic solver) would no longer be required.

To conclude this section, we highlight the possibility in SDPD to incorporate complexity on a coarse-grained level, i.e., on scale larger than the suspended polymer molecule. We denote this alternative approach as coarse-grained viscoelasticity.

In Refs. [19] and [91], a coarse-grained fluid-particle model for a polymer solution was proposed. Instead of linking particles together as discussed above, here every SDPD particle is considered as a thermodynamic sub-system containing many polymer molecules (see Fig. 2). The state of the fluid particle is characterized by a configuration tensor c that describes their average molecular elongation and orientation. The specification of very simple physical mechanisms inspired by the dynamics of single polymer molecules allows one, with the help of the GENERIC formalism, to derive the equations of motion for this set of SDPD particles carrying polymer molecules in suspension which satisfy strictly thermodynamics consistency, i.e., the first and second laws of thermodynamics and FDT.

In a nutshell, if we consider a set of fluid particles labeled by Latin indices i, j=1, ..., N, in the most general case, the GENERIC-derived stochastic differential equations for the conformation tensor associated with each SDPD particle read

(20)

where we have restricted our attention to the isothermal case (for the general non-isothermal situation, the reader is referred to Ref. [19]). It can be shown that the resulting momentum equation is the same as Eq. (1), where the isotropic pressure tensor is replaced with a general anisotropic viscoelastic tensor taking into account the additional contribution due the microstructural state of the polymers, i.e., the conformation tensor c, in the form

(21)

where σ_i is a tensorial variable thermodynamically conjugated to c_i, i.e., , where T is a constant temperature and S_p(c) is the conformational-dependent entropy function for the sub-system depicted in Fig. 2. The first two terms in Eq. (20) take into account flow-induced stretching effects on the polymer molecules, the third term is an irreversible relaxation contribution, and the last term is a random contribution connected to the Brownian motion of the polymer molecules contained in a SDPD particle (see Fig. 2). Here, λ is the longest relaxation time of the polymer.

Expressions (20) and (21) are of general validity as no assumption has been made yet on the specific force law of the polymer. Compliance with the GENERIC structure ensures their thermodynamic consistency at the discrete level.

Polymer physics comes into play in this model with a proper definition of S_p(c). In the simplest case of a dilute suspension of Hookean dumbbells, this entropy reads^[92]

(22)

where k_B is the Boltzmann constant, and N_p is the total number of dumbbells per fluid particle. In Ref. [19], we have shown that this specific choice leads to σ_i=(c_i^-1-1) and the resulting expression for the polymeric stress is

(23)

where we have defined the polymeric viscosity as η_p=N_pd_ik_BTλ. Finally, the last term on the right-hand side of the evolution equation (20) for c_i reduces to

(24)

The deterministic part of the resulting equations in Eq. (20) can be therefore interpreted as a very specific SPH discretization of the classical Oldroyd-B constitutive model for a dilute suspension of Hookean dumbbells. This case is reached in the limit of large SDPD particle size by virtue of the correct scaling of thermal fluctuations^[19]. For small SDPD particles, stochastic terms become relevant, and Eq. (20) with this specific choice of S_p(c) delivers a stochastic generalization of the Oldroyd-B model.

It should be borne in mind that, in such a derivation, no reference to the target partial differential equations (PDEs) (i.e., Oldroyd-B) was considered. The fact that a SPH discretization of an Oldroyd-B equation is finally recovered represents a " posteriori" proof of the consistency of the approach, as it is the expected result for Hookean dumbbells in suspension. Generalization to more complex polymeric models, such as finitely extensible nonlinear elastic springs, with the proper introduction of thermal fluctuations is straightforward. In particular, coarse-grained thermodynamic consistent models can be constructed by physical specification of conformation-tensor-dependent entropy of the fluid particles appearing in Eq. (22) rather than by brute force discretization of existing continuum constitutive equations.

In summary, in SDPD, fluid complexity can be introduced both (ⅰ) at the mesoscopic particle-level (similar to conventional DPD) or (ⅱ) at a supra-particle coarse-grained level. In the latter case, supplemental coarse-grained variables associated with each particle (in addition to positions, momenta, energies, and entropies) need to be considered.

11 What is the relation between SDPD and classical CFD methods?

SDPD is a method for the description of fluctuating hydrodynamics characterized by finite inertia and compressibility. There is a large recent effort to extend classic CFD methods like finite differences, finite volumes, and finite elements in order to numerically solve the equations of fluctuating hydrodynamics ^[93-96]. In common with any other CFD method, SDPD faces the problem of disparate time scales. In fact, a number of time scales emerge associated with different processes in a fluid. For example, if we have a small particle of radius R in a fluid characterized by a speed of sound c_s and kinematic viscosity ν, i.e., moving with a typical velocity V, we may identify the following time scales: the sonic time it takes a sound wave to cross the particle's radius t_s=R/c_s, the vorticity time it takes for vorticity to diffuse across the particle's radius t_ν=R²/ν, the convective time t_V it takes for the flow velocity to travel the particle's radius t_V=R/V, and the diffusive time t_D that takes the particle to diffuse across its own radius t_D=R²/D. The ratio of time scales gives rise to the usual dimensionless numbers, e.g., Mach number Ma= t_s/t_V=V/c, Reynolds number Re=t_ν/t_V=R²V/ν, Schmidt number Sc=t_D/t_ν=ν/D, and Peclet number Pe=t_V/t_D. For a micron sized particle R=1μm in water, these numbers are either very small or very large, reflecting a huge separation of time scales, e.g., t_s ≈ 1 ns, t_V ≈ 1 μs, t_d=R²/D ≈ 1 s, typically leading to Ma ≪ 1, Re ≪ 1, and Sc ≫ 1. When we have such large/small dimensionless numbers, the governing equations become stiff. As the numerical time step is dictated by the fastest process, simulation of the longest time scales requires an impracticable number of time steps. Two different strategies have been considered in order to address the problem of separation of time scales in CFD. The first strategy is to consider limiting equations, the second one is to use telescoped equations.

In the first strategy, one takes rigorous mathematical limits and describes the fluid with a different set of equations, usually of elliptic character. For example, under microflow conditions, the Reynolds number could be typically on the order of (10^-5-10^-2) ≪ 1 in such a way that the inertia-less Stokes approximation (Re=0) is typically considered when solving the momentum equation. Moreover, liquids under normal flowing conditions can be considered as incompressible, i.e., the pressure field adapts instantaneously to a given flowing conditions, which is equivalent to taking the limit of the speed of sound c_s→∞ and Ma=0. Finally, in diffusion problems, the Schmidt number Sc in real fluids is typically on the order of (10³-10⁶)², and one can take an infinite Schmidt number limit and solve the corresponding equations^[97]. The limiting equations overcome the stiffness problem and the need to resolve very fast processes, allowing to use much larger time steps and hence much larger time scales. The limiting equations are, of course, approximations to reality that are computationally convenient.

²In polymer dynamics, for example, this condition is critical to guarantee that hydrodynamic interactions between monomers of a chain are fully developed on time scale of mass diffusion and to recover the correct Zimm scaling of polymer relaxation times^{[12, 14]}.

A different strategy to tackle the problem of separation of time scales is to use telescoped equations in which the values of the dimensionless numbers, while not complying with realistic values, still allow for a clear separation of time scales. Typically, ratios of (10-10²) are usually considered between different time scales. For example, one can choose a speed of sound c_s sufficiently large to reduce density fluctuations below few percents but still not realistically large to affect the efficiency of the time integrator. In the context of mesoscopic techniques, this strategy has been named telescoping time scales by Padding and Louis^[98].

Each strategy to address the problem of separation of time scales in hydrodynamics has its pros and cons. While the limiting equations admit time steps much larger than the telescoped equations, the solution of idealized elliptic equations (e.g., Poisson equation) needs sub-step iterations to get converged solutions within machine precision and therefore increase CPU time per step compared to explicit time integrators. In addition, elliptic problems are non-local and notoriously difficult to solve in combination with efficient parallel implementations. Scalability on massively parallel hardware is an important computational aspect to keep in mind when solving complex three-dimensional applications where efficient parallelization of the simulation approach is essential.

As remarked, SDPD is a CFD method to solve the stochastic equations of fluctuating hydrodynamics, and the above two strategies have been pursued in SDPD. For example, limiting equations have been considered by including a constrain on the density (thus changing the equations) that ensures incompressibility of the fluid^[99]. Although this approach was initially proposed for SPH, it can be straighforwardly applied to the stochastic SDPD. In fact, the Lagrange multipliers used to constrain the density enter uniquely the reversible part of the discrete dynamics without affecting irreversible/stochastic terms. Other possibilities to enforce incompressibility (zeroing the divergence of the velocity) in SPH include projection methods^{[60, 100]}. As the constraint here is on the velocity field itself, if one wants to use projection methods on SDPD, one needs to ensure that fluctuations obey the FDT, which is not obvious. Both approaches allow to bypass the Courant-Friedrichs-Levy stability constraint on the speed of sound (Δt≤ r_c/c_s). Also, novel splitting-implicit strategies have been introduced in SDPD for the description of the solvent which bypass the viscous stability time step condition (Δt≤ r_c²/ν), particularly strict under microflow conditions^[13]. In the context of SPH, highly viscous flows in the Stokes approximation have been recently solved by using conjugate gradient method^[101]. As mentioned, all these new approaches involve some kind of non-local effect which shows up numerically as a number of sub-time steps iterations required. This increases the CPU load per time step, and might introduce some technical difficulties in efficient parallelizations.

The strategy of telescoped equations has been also considered in SDPD. The weakly compressible SPH method to deal with the small Mach number problem is a popular method^[61] that translates itself into SDPD. The telescoped equation strategy in a particle method leads to local interactions that are highly suitable for high performance computing (HPC) calculations where optimal performance up to hundreds of thousands of cores can be easily obtained^[102].

We return now to the main question of the heading of this section that we can rephrase in the form: what are the benefits/drawbacks of particle-based Lagrangian as compared to grid-based Eulerian methods? In macroscopic flows, SPH has found a clear niche in the treatment of astrophysical flows where one solves only in the spatial domain where the fluid is or in violent flows where splashing, sloshing, and fragmentation are naturally described while Eulerian methods painfully suffer to capture these effects^[28]. In the mesoscopic realm, free-surface flow problems are less common, and these benefits do not seem critical. However, the Lagrangian character of a technique can still allow to track interfaces implicitly (i.e., without additional surface tracking algorithms) which is important for several problems in microfluidics such as dynamics of emulsions, deformable particles, and biological cells. In addition, when used in combination with viscoelastic models (the coarse-grained viscoelasticity approach), the microstructural state of the liquid does heavily depend on element flow history which is directly available in Lagrangian methods. Another appealing aspect of a particle method is the possibility to easily incorporate new physics through additional forces between the particles satisfying default standard conservation laws, general thermodynamics, and FDT at the discrete level. The use of the particle paradigm allows to benefit from the large know-how accumulated in molecular dynamics simulations. On the other hand, in Eulerian grid based methods, it is necessary to couple the hydrodynamic description of the solvent using a given discrete approach with a different treatment of the solid-liquid interface and fluid-structure-interaction (e.g., by coupling with molecular-dynamics-like approaches for suspended objects by using immersed boundary, smooth profile methods, etc.) for which the exact enforcement of thermodynamic consistency at the discrete level remains a challenge^[103].

In conclusion, the different simulation approaches have pros and cons in regard to efficiency and accuracy which are difficult to assess based on formal discussions. To our knowledge, there have been no direct comparisons between SDPD and Eulerian methods so far. The value of each technique should be judged a posteriori by the level of complexity reached by the target applications, the accuracy in describing them (by systematic comparison with experiments, whenever possible) and, more importantly, by the ability of each technique in predicting new physics, possibly quantitatively, in situations when no a priori knowledge of the problem is available. This is ultimately the "holy grail" of computational science.