A discretized integral hydrodynamics

Using an interpolant form for the gradient of a function of position, we write an integral version of the conservation equations for a fluid. In the appropriate limit, these become the usual conservation laws of mass, momentum and energy. We also discuss the special cases of the Navier-Stokes equations for viscous flow and the Fourier law for thermal conduction in the presence of hydrodynamic fluctuations. By means of a discretization procedure, we show how these equations can give rise to the so-called"particle dynamics"of Smoothed Particle Hydrodynamics and Dissipative Particle Dynamics.

Using an interpolant form for the gradient of a function of position, we write an integral version of the conservation equations for a fluid. In the appropriate limit, these become the usual conservation laws of mass, momentum and energy. We also discuss the special cases of the Navier-Stokes equations for viscous flow and the Fourier law for thermal conduction in the presence of hydrodynamic fluctuations. By means of a discretization procedure, we show how these equations can give rise to the so-called "particle dynamics" The inherent difficulties of the equations for fluid hydrodynamics has given rise to a variety of schemes that numerically simulate the behavior of a fluid at hydrodynamic scales [1]. Of renewed interest, there have appeared the simulation schemes based on the idea of substituting the fluid by "particles" that move under the influence of forces that, in a coarse-grain limit, reduce to some approximate form of the hydrodynamic Navier-Stokes equations; these schemes are called Smoothed Particle Hydrodynamics (SPH) [2,3] and Dissipative Particle Dynamics (DPD) [4]. Although their numerical implementation is somewhat different, the latter including a random force, we shall argue that their origin is essentially the same. The use of these particle-like simulations has been reported of being successful in different applications of fluid dynamics [5,6,7,8,9,10]. An attractive feature of these simulations is that one can use the enormous experience gained from the techniques of simulations of standars molecular dynamics; in particular, in dealing with rheological fluids.
In this article, we present an integral representation of the hydrodynamic conservation laws based on the concept of the interpolant of a function [2,3]; the interpolant is an integral representation of a function in terms of a weighting function. This function, in the appropriate limit, becomes a delta function and the interpolant yields an identity. For calculational purposes one does not take such a limit and, therefore, one ends with an approximating form of the corresponding function. The SPH simulation is based on a certain form of interpolant of the hydrodynamic functions, such as density, velocity, etc. Here, we base our scheme not on the interpolants of the functions, but rather on the interpolants of the gradient of the functions. This minor change proves to be very useful in writing integral equations that reduce, exactly, to the conservation laws in the limit. Moreover, we shall see that the conservation of mass, momentum and energy is exactly preserved at the level of the integral forms. At an approximation level, we shall show that the integral forms of the conserving currents are correct up to second order terms in the gradients of the corresponding fields. Since the integral expressions can be written down following the forms of the true laws, the phenomenological variables, such as viscosities and thermal conductivity, can be naturally included. In the same fashion, the extension to fluctuating hydrodynamics can be readily performed.
The resulting integral conservation laws may then be used as an alternative to the exact differential laws. Further, they can be utilized as the starting point for approximate solutions. In this context, we show that by an appropriate discretization of the integrals one can render the equations to look as the equations, not for the hydrodynamic fields evaluated at space-time points, but rather, for quantities pertaining to particles. For instance the field velocity v(r, t) becomes the velocity v i (t) of a particle at position r i (t); a law of motion for the latter must then be supplied. Within the present theory one readily finds the equations of motion for the particles that are fully consistent with the hydrodynamic equations. We shall show how the SPH and DPD equations may be found. As we shall see, in general, one can have additional terms arising from the convective non-linear terms of the hydrodynamic time derivative.
We organize the article as follows. In the next section we introduce the interpolant of the gradient of a function. With such an object, we write down the conservation laws, both in general and for the cases of Navier-Stokes and Fourier law. The extension to fluctuating hydrodynamics is also shown. In section III we present a discretization procedure that yields a particle-like simulation algorithm. We discuss the relationship with the SPH and DPD theories. We conclude and make additional remarks in section IV.

II. HYDRODYNAMICS IN INTEGRAL FORM
A. An interpolant of the gradient The idea of formulating an integral version of the laws of hydrodynamics is based on a limiting representation of the gradient of a function of position r. This representation we call it an interpolant, following Lucy [2] and Monaghan [3] SPH formulation. As we mentioned in the Introduction, we use the interpolant of the gradient of a function rather than that of the function itself.
First, we show that the following identity is correct, where the tensorial character of A is left unspecified, and W (|r − r ′ |; r 0 ) is a distribution or weighting function sharply peaked with width r 0 . We demand that all its moments exist although the function itself may not be integrable. We assume the following: with n ≥ 1. For n = 1 we require M 1 = 3 but for n > 1 we leave M n unspecified. These requirements are easily satisfied by noticing that the integrand of (1) can be related to the gradient of a distribution that tends to a delta function, namely where f (r; r 0 ) is such, that in the limit r 0 → 0 becomes lim r0→0 f (r; r 0 ) = 1 4πr 2 δ(r).
Clearly, any distribution that tends to δ(x) can be used. It is important to stress that W is a function of the magnitude |r| and not of the corresponding vector. This dependence is important for the use of approximations and, as we shall see later on, for the setting up of the conservation laws. For the validity of (1) we further require that A(r) is analytic everywhere inside the domain of the integral. This is not a stringent requirement since we are interested in hydrodynamic fields. Calling I(A(r)) the integral in equation (1), we first make a change of variable, r ′ → r + r ′ , and then, we perform a Taylor expansion of A(r + r ′ ) around r. We obtain where ∂ α = ∂/∂r α with r α the cartesian components of the vector r. Due to the spherical symmetry of W all the odd powers of r ′ (even powers in the derivatives) vanish identically. With the use of Eq. (2) the first integral in the rhs of (5) gives the Kronecker delta δ αβ , independently of r 0 , while the second integral yields r 2 0 (δ αβ δ γη +δ αγ δ βη +δ αη δ βγ )W 3 /15. Clearly the higher order terms are proportional to r 2n 0 times odd (2n + 1)-derivatives of A. We write, generically, In the limit r 0 → 0, the interpolant is the gradient of A. For approximation purposes, we note from (6) that the correction to the gradient is third order in the derivatives.

B. Conservation laws
The study of the hydrodynamics of a fluid is based on the conservation laws for mass, momentum an energy. [11] In the following we assume that all the fields are evaluated at a spatial point r and that all are time dependent as well. The conservation laws are where ρ is the mass density; v is the velocity of the fluid; j is the momentum density, j = ρv; P is the hydrostatic pressure;Π is the viscous stress tensor; e = ρv 2 /2 + u is the total energy density with u the internal energy density of the fluid; and J is the heat current. Integration over the whole volume of the fluid shows that the total mass, momentum and energy of the fluid are conserved.
With the use of the interpolant given by Eq.
(1), we can write down analogous conservation expressions (with implicit time-dependence): The above expressions, by the rule of the interpolant, are equal (in the limit r 0 → 0) to (minus) the divergence of the terms in square brackets evaluated at r ′ = r; in particular, the "kernels" P (r; r ′ ),Π(r; r ′ ) and J(r; r ′ ) must be chosen such that the interpolant equal the divergence of the actual pressure, viscous stress tensor and heat current when r ′ = r. But before we discuss how to choose these kernels we point out that all the terms inside the square brackets are symmetric with respect to the interchange of the variables r and r ′ . Therefore, integration with respect to r makes the rhs of all the equations vanish, thus yielding the conservation of the extensive variables, independently of whether the limit r 0 → 0 is taken or not.

C. Constitutive relations
The conservation laws must be provided with constitutive relations in order to have a closed set of equations. For discussion purposes we shall choose the mass density, velocity and temperature T (r, t) as the independent fields. Therefore, for the pressure and the internal energy we need to know the equations of state of the fluid in terms of ρ and T ; in particular, we assume we know the functional dependence P (ρ, T ). Thus, the kernel for the pressure may be chosen as Clearly, the kernel is not symmetric in its variables and hence we need the symmetrization in the conservation equations (10) -(12) above. Substitution of this form into, say, Eq. (11), gives to lowest order, the term where the rhs is evaluated at r ′ = r. The rhs is −∇P (r). Regarding the viscous stress tensorΠ, we use the usual one that gives rise to the Navier-Stokes equations, linear in the velocity gradients. The corresponding kernel may be written as where the viscosities η(r) and ζ(r) are, either, given functions of r, or they depend on r through a further dependence on density and temperature. Substitution of this equation into Eq.(11), yields to lowest order the familiar viscosity terms of the Navier-Stokes equations: For the form of heat current J we consider Fourier law in which the current is linear in the temperature gradient.
The kernel can be written as, Again, substitution into Eq. (12) yields, to lowest order, where the r-dependence of κ may be given through its dependence on temperature and density. Thus, we have shown a consistent way of presenting an integral form of the equations of hydrodynamics, such that, in the appropriate limit yield the true equations; we note that all the phenomenological coefficients are readily and unambiguously identified. This is an important point since, as we shall show in the next section, the choice of the functional forms of the stress tensor and the heat current are, by no means, unique; that is, we can prescribe different functional forms that in the limit also give similar expressions to the usual hydrodynamic laws. We shall defer further discussion of this point to section III.

D. Hydrodynamic fluctuations
We now turn our attention to the formulation of the study of hydrodynamic fluctuations. Following Landau and Lifshitz [12,13] we limit ourselves to small fluctuations around a given flow, solution to Eqs. (7) -(9), so that a linearization in the fluctuations is possible. In keeping with our assumption that the indepedent variables are the mass density, velocity and temperature of the fluid, we define the fluctuations as linear deviations from the flow, that is and analogous expressions for v(r, t) and T (r, t). The functions ρ 0 (r, t), v 0 (r, t) and T 0 (r, t) constitute a given flow, solution to the full non-linear integral equations (10) - (12), with the expressions (15) and similar linear equations for the (partial) time derivatives of momentum density and energy, δj and δe. Next, one identifies the source of the fluctuations as arising from spontaneous fluctuations of the stress tensor and the heat current. This is implemented in the usual way, [12] by adding to the linearized equations for the fluctuations of the momentum and energy densities, terms proportional to the divergence of a random stress tensor and to the divergence of a random heat current, respectively. To be precise, we add to the equation for the fluctuation of the momentum density, where the tensorΠ R (r, r ′ , t) is a gaussian random stochastic function, symmetric under interchange of r and r ′ , with zero mean and with its second moment obeying the usual fluctuation-dissipation relations, In the same fashion, we add to the equation for the fluctuation of the energy density, where J R (r, r ′ , t) is a gaussian random stochastic function, symmetric under interchange of r and r ′ , with zero mean and with second moment We recall that the quantities δρ(r, t), δj(r, t), etc. depend from those of the underlying flow, ρ 0 (r, t), j 0 (r, t), etc. but not the other way around. That is, one first solve for the latter and then one finds the fluctuations. It is understood that the flow is stable; that is, one should always have δρ/ρ 0 ≪ 1, etc. In the next section we shall also comment how one can include the fluctuations within a particle-like simulation.

III. A DISCRETIZED INTEGRAL HYDRODYNAMICS
The integral formulation presented in the previous section is simply an approximate representation of the usual hydrodynamic laws. Their usefulness resides on whether their solution may be easier to find than those of the actual equations or on their amenability for approximations. In this regard, we recall the approximate SPH and DPD schemes where a particle-like simulation, similar to a molecular dynamics simulation, represents the flow of a continuous fluid. In this section we present a particular discretization of the integral equations of section II that may be used as the basis for a simulation in terms of "fluid particles."

A. A particle-like scheme
The basic idea is first to divide space in cells of finite size ∆V and, then, to define the field variables for each cell. We call r i the position vector of the i-th cell, and the following list summarizes the variables for such a cell: Now, the kernels of pressure and viscous stress become "potentials" of force between the i-th and j-th cells while the heat density current is now a "current" of energy between such cells: P (r; r ′ ; t)∆V → P ij (t) pressure potential Π αβ (r; r ′ ; t)∆V → π αβ ij (t) stress potential J α q (r; r ′ ; t)∆V → J α ij (t) heat current (27) The integrals are then discretized by summing over cells directly and not over labels that localize the cell in a Cartesian grid, that is, where we have assumed there are N cells in the total volume. This discretization implies a careful choice of the discretized version of the weighting function W (|r − r ′ |). That is, we cannot simply change r by r i and r ′ by r j in the functional form of W since the equations for the moments, Eq. (2), would not be correct. This is due to the fact that those results make use of the spherical symmetry of W . Instead, we propose the following discretization that give rise to the correct moments: where we have defined r ij = |r i − r j |. This form also takes into account that W is always part of an integrand. As a particular example, using as a representation of a delta function, yields for W (r), and, correspondingly for W(r ij ), which shows that, in discretized form, all the moments but the zeroth are well defined. With the above reformulation, and definingê the conservation equations look as follows. Conservation of mass: Conservation of momentum: Conservation of energy: By construction, the total mass, total momentum and total energy are conserved. This can be seen by summing the above expressions over i = 1, 2, . . . , N .
So far, the equations are quite general and one needs constitutive relations for the kernels of pressure, viscous stress and heat current. For instance, the viscous stress tensor linear in the velocity gradients may be found by discretizing Eq. (16), while the heat current may be found from Eq.(18) One should keep in mind that in order to close the equations one still needs the equations of state for the pressure P ij = P(m i , T j ) and the internal energy per particle u ij = u(m i , T j ). Up to here it is simply a discretization of the equations. The interesting addition now [2,3,4] is to assume that the positions r i of the cells become the positions of particles that are allowed to move. This is certainly a bold assumption, since making the fluid particles move should be done in a Lagrangian formulation of the fluid dynamics rather than in an Eulerian one. The present are a different type of particles, however, since as we can see from the above equations, they not only change their momenta but they also have variable mass and carry with them their internal energy in addition to their kinetic one. Nevertheless, one may justify it by arguing that one is actually looking at every instant of time to a state of the fluid not on a grid but rather on a "fluidized grid"; the motion law of r i being used as an updating of the state of the fluid on such a grid. In any case one can asses the validity of such an assumption a posteriori; as we mentioned in the Introduction, successful simulations of actual flows with schemes like the present one have been reported, see Refs. [5,6,7,8,9,10]. Although the choice of the law of motion is arbitrary, it seems "natural" to consider the rate of change of the position of the cell as the velocity of the fluid at that point: We point out that this choice is not unique; see Ref. [3] for other forms used in SPH simulations. In principle the above scheme is complete and closed. However, for an actual implementation of a simulation based on it there are further questions to resolve, such as the boundary conditions in terms of the particles and the discretization of time. Since there are already in the literature a host of procedures [7,16] both for dealing with boundary conditions between particles and solid frontiers and for the time discretization, we shall only discuss the latter because of its relevance in the inclusion of hydrodynamic fluctuations.
A simple algorithm to simulate the dynamics consists of a two-step propagation in time [4]. First, there is a "collision" step in which one finds the values of m i , p i and e i at time t + ∆t from from knowledge of all the variables at time t, using Eqs. (34) -(36), with It is then followed by a "propagation" step in which the positions r i (t + ∆t) are computed using This combination is more accurate than if both steps were done with Eq. (40). This algorithm, however, is also useful to include the fluctuations as part of the evolution and not as a posteriori calculation, being thus helpful in determining the stability of the flow. This is an important point since the purpose of the simulations is to solve the equations by an actual propagation in time of the flow. (An analytical solution need not be done in this way; for instance, if the equations are linearized one may solve them using an integral transform technique.) Thus, one can include in the equations for the momentum and the energy, Eqs. (35) and (36), discretized versions of the random viscous tensorπ R ij , and random heat current J R ij , both symmetric in ij. Since their second moments must obey discretized versions of the fluctuation-dissipation expressions (23) and (25) these tensors may be added as π Rαβ ij (t) = (2kT i (t)η j (t) + 2kT j (t)η i (t)) 1/2 ∆X αβ ij (t) where ∆X αβ ij (t), ∆Y αβ ij (t) and ∆Z α ij (t), symmetric in ij, represent independent random increments (Wienner processes [14]) with zero mean and correlations From a practical point of view, any good commercial pseudo-random number generator suffices for these increments.
The important aspect we want to stress is that, in contrast to the continuous version (23) and (25), here the temperature, viscosities and thermal conductivity that appear in Eqs. (42) and (43) are evaluated at the current values of the full fluctuating quantities and not at the values of the variables at the underlying flow. Since it is assumed that the fluctuation are always small and do not make the flow unstable, this is a minor approximation. Moreover, as mentioned above, if the flow do become unstable by adding the random viscous tensor and heat current that may imply, barring numerical inaccuracies, that either the flow is indeed unstable or that the method itself does not faithfully describe the flow. It should be clearly understood that the particle-like representation of a continuous fluid flow depends on two different approximation; first, one approximates the true differential laws by integral expressions with a finite width r 0 of the weighting function; and second, the integrals are discretized. These approximations pose constraints on the length scales of the fluid. On the one hand, the density of point-particles must be such that the mean particle separation (∆V ) 1/3 is smaller than r 0 in order to have a good approximation of the integrals; and on the other hand, a typical hydrodynamical length, call it λ, must be larger than r 0 itself in order to have a good representation of the gradients in terms of the integrals (i.e. an independence of the parameter r 0 , see Eq.(6).) That is, one should always have, This way, the limit r 0 → 0 implies not only the equality of the integral and differential forms of the conservation laws, but also the continuum limit itself.

B. SPH and DPD as special cases
The purpose of this section is not to make a revision nor a comparison of SPH and DPD schemes with the present one, but rather, to show that they may be viewed as special cases of a more general scheme that reduces to the macroscopic conservation laws of fluids.
As we have seen the discretized conservation equations (34)-(36) are very general. One still needs to provide constitutive relations for the pressure, viscous tensor and heat current and, as long as ij-symmetrized forms are provided, the conservation laws are guaranteed. The particular expressions given in the previous section, such as (37) and (38), are just examples. But before we present other forms used, such as those of DPD and SPH, we want to mention some aspects of the time derivatives used.
In the schemes used in SPH and DPD, the time derivatives of the properties of the particles have been interpreted as already including the convective contribution. In the discretized version this is equivalent to identify This is a subtle point: One the one hand, one could argue that the derivative is following the motion of the fluid particle, like in a Lagrangian formulation; however, the right-hand-side of the corresponding conservation laws (34)-(36) should be accordingly transformed. Since the latter is not done in SPH and DPD, one may still say that those formulations correspond to not too large Reynolds numbers where one can neglect the convective contributions. It may be interesting to include those terms explicitely in a simulation.
With the above identification of the time derivatives, we can now see a closer resemblance to the equations of SPH and DPD. For purpose of exemplifying the relationship, we shall only discuss the equation for the momentum. Using the mass conservation equation and the fact that p α i = m i v α i , the equation for the momentum can be written as, In Refs. [2,3,5], SPH is formulated with forms for the pressure such as with a and b constants and with a given equation of state for P i in terms of m i . [15] The viscous stress tensor of SPH and DPD may be generally written with appropriate choices of A and B [4,7,8,9,10]. It is interesting to note that this form can give rise, in the continuum limit r 0 → 0, to terms proportional to ∇ 2 v and ∇(∇·v); however, one cannot independently identify the corresponding viscosity coefficients. This is to be contrasted with the expression of the tensor given by equation (37) where there is an independent identification of the viscosities. In the DPD simulations there is an additional ingredient. Namely, that the pressure term is taken to be stochastic. Within the present scheme this may be interpreted as including hydrodynamic fluctuations with white noise and with a particular temperature and viscosity as given by equation (42). In this regard we differ from the interpretation of DPD equations given by Español et al. [9] and Marsh et al. [10]. In that interpretation, the equations of DPD are taken as the "microscopic" dynamics of the particles of a fluid, from which the macroscopic laws are to be extracted, much in the spirit of Langevin and Boltzmann equations. Within that interpretation they have argued that the random part should be modified in order to account for the correct fluctuation-dissipation relation of Langevin-like equations. According to the present theory, since the equations of motion of the particles are only an approximation to the macroscopic equations of the fluid flow, there is no need to modify the DPD equations. Therefore, the random contributions of DPD may be seen to already refer to hydrodynamic fluctuations. Moreover, if one wishes to find the corresponding Fokker-Planck equations to the discretized hydrodynamic equations (34) -(36) one can follow the theory of van Saarloos et al. [13] of non-linear hydrodynamic fluctuations.

IV. FINAL REMARKS
In this article we have presented an integral form of the conservation laws of a macroscopic classical fluid in terms of an interpolant for the gradient of a given function of space. This form is amenable for a discretization of space and may be interpreted in terms of the dynamics of "fluid particles". To complete this discretized dynamics one must provide the law of motion of the position of the particles; one may prescribe that the field velocity equals the rate of change of the position of the particle. Within this scheme one can easily find out the corresponding Navier-Stokes equations of viscous flow and the Fourier law of heat conduction. Moreover, hydrodynamic fluctuations can also be readily taken into account.
We have argued that numerical simulations currently used, known by Smoothed Particle Hydrodynamics (SPH) [2,3] and Dissipative Particle Dynamics (DPD) [4], and which are based on a particle-like simulation of a continuum fluid, may be seen as special cases of the general formalism here presented. In that way, one the one hand, one guarantees that the simulations have a correct continuum limit, and on other, there is a clear route of how to represent, in the particle dynamics, known effects of macroscopic fluids; for instance, with the present theory one can see how to include thermal effects, absent in DPD simulations.
Finally, we want to stress the potential uses of this type of schemes. It does seem that a complication in any simulation of fluids is the discretization of space with its concomitant difficulties of boundary conditions; this is the more important if one is interested in rheological fluids, such as suspensions. That is, in order to simulate a simple flow a discretized differential scheme (e.g. finite differences) may appear to be better than a discretized integral version; this is because the latter makes use of the weighting function which in turn must resemble a delta function, and therefore, it appears that one needs more "particles" than points in a grid [6], cf. Eq. (45). However, having paid this price, there is a host of "tricks" and techniques, borrowed from standard molecular dynamics, that can be used to simulate moving boundaries and solid objects, e.g. Lee and Edwards shear boundaries [16], or "freezing" certain number of particles to simulate a rigid body [4], etc. Moreover, the inclusion of hydrodynamic fluctuations also seem to be much easier within a particle-like scheme than within a field-like.