Entropy stability analysis of smoothed dissipative particle dynamics

This article presents an entropy stability analysis of smoothed dissipative particle dynamics (SDPD) to review the validity of particle discretization of entropy equations. First, we consider the simplest SDPD system: a simulation of incompressible flows using an explicit time integration scheme, assuming a quasi-static scenario with constant volume, constant number of particles, and infinitesimal time shift. Next, we derive a form of entropy from the discretized entropy equation of SDPD by integrating it with respect to time. We then examine the properties of a two-particle system for a constant temperature gradient. Interestingly, our theoretical analysis suggests that there exist eight different types of entropy stability conditions, which depend on the types of kernel functions. It is found that the Lucy kernel, poly6 kernel, and spiky kernel produce the same types of entropy stability conditions, whereas the spline kernel produces different types of entropy stability conditions. Our results contribute to a deeper understanding of particle discretization.


Introduction
The smoothed dissipative particle dynamics (SDPD), which was proposed by Español(2003) [1], has attracted the attention of many physicists and engineers for over a decade. The SDPD is appropriate for simulations of mesoscale flows in which thermal fluctuations are non-negligible, and it has been acknowledged in many complex flow problems (e.g. cellular blood flows [2,3] and suspension flows including polymer molecules [4]). By introducing conservation of angular momentum, the accuracy of SDPD simulations of incompressible flows has been greatly improved [5,6]. Thus, the SDPD exhibits great potential for solving many types of thermal flow problems that exist around us.
As described in the literature [1], the governing equations of hydrodynamics with thermodynamic consistency are Here, ρ,κ,η, andζ are density, thermal conductivity, sheer viscosity, and bulk viscosity, respectively. The set of equation (1) and equation (2) represents the Navier-Stokes equations, and equation (3) represents the relationship among the entropy s, viscous heating field f, and temperature field T. Hereinafter, we refer to equation (3) as the 'Batchelor's equation', as in [7]. Because the SDPD is a kind of Lagrangian particle method, each of equation (1) to equation (3) is discretized using particles in a similar manner to smoothed particle hydrodynamics (SPH) [8,9]. According to equation (32) in [1], the discretized expression for the rate of change of the total entropy S is represented as Here, T i and T j are the temperature of the ith and jth particles, respectively. The viscous heating field f i , parameter d i , relative temperature T ij , and relative function F ij are given as where the relative vector r ij , unit vector e ij , relative velocity v ij , kernel function W, and its gradient function F are given as Additionally, the parameter h(h>0) is the kernel radius that determines the interaction range between particles. The function W (Lucy kernel function [10]) in equation (12) and the function F in equation (13) are the reposts of equations (9) and (10) in [1], which presents the original SDPD. Equation (4) represents the introduction of a weighted calculation depending on the positions of particles using the function F. Although computational physicists have discussed the validity of applying a particle discretization scheme to compute physical quantities in the case of density calculation of SPH [8,9] or that of moving particle semi-implicit(MPS) [11,12], the validity of applying a particle discretization scheme to the entropy calculation in equation (4) has not been corroborated in terms of theoretical thermodynamics, even though the mathematical derivation of equation (4) and its expression in the GENERIC framework [13][14][15] were discussed in the original research [1].
The primary purpose of this article is to discuss the thermodynamic validity of equation (4) in SDPD. Our strategies are as follows. First, we consider a quasi-static scenario with constant volume, constant number of particles, and infinitesimal time shift, assuming the moment of an explicit SDPD simulation of incompressible flows. Next, we derive a form of entropy from the discretized entropy equation of SDPD by integrating it with respect to time. We then examine whether the SDPD system exhibits physically reasonable behaviours by performing thermodynamic entropy analysis.
The remainder of this article is structured as follows. Section 2 derives the form of entropy from equation (4). In section 3, we examine the characteristics of a two-particle system of SDPD by examining entropy stability conditions. Finally, section 4 summarizes our results and concludes the article.

Derivation of entropy
Let us consider the case that time t varies from t 0 to τ, and temperature T i of the ith particle(i=1, 2, K, N) varies from T i 0 to t T i in Kelvin. Let the time τ be between t 0 and t 0 +Δt. In this paper, we focus on the simplest SDPD system, which is an explicit simulation of incompressible flows using a forward-Euler time integration scheme [16]. Under this premise, the left-hand side of equation (4) is approximated by the ratio of infinitesimal entropy ΔS to infinitesimal time Δt (dS/dt≈ΔS/Δt) in the simulation space. The positions and velocities of all particles are fixed during the interval of Δt. Also, the values of d i , d j , and F ij become constant in this interval. Accordingly, the viscous heating field f i of equation (5) becomes a constant field during Δt. To summarize, the time-dependant variables on the right-hand side of equation (4) become T i , and T j ; hereinafter, we denote the temperature T i as T i (t) when we need to expressly show the time-dependency.

General case of N particles
We derive the form of entropy S by integrating both sides of equation (4) over time: By changing the order of integration and summation, equation (14) can be rewritten as The first term of equation (15) can be rewritten using the definition of definite integration by substitution [17] as To integrate the second term, we introduce the concept of general topology. Let us consider the following relationship between a differentiable function ( ) G r from some open subset U ( n ) to  and a differentiable function r from some closed interval to U. Then, by the multi-variable chain rules [18], we get Subsequently, we can integrate both sides of equation (17) to get Let us discuss the physical interpretations of the small element dr and the domain of integration when we set r to be (T i , T j ). Because we perform the integration over a small time interval, we can regard dT i /dt and dT j /dt as constant values in the first approximation level. Due to the symmetricity of i and j in time dependence, the relationship of dT i /dt=dT j /dt=K=const can be established. Hence, equation (18) can be expressed using the vector = ( ) u 1, 1 as follows: Meanwhile, we can choose ( ) G r as , we can exclude this point from the domain of integration because temperature in Kelvin is always positive. We then obtain the following relationship: Here, we use the relationship of In the case that we integrate the left-hand side of equation (21) with respect to time t in the range [t 0 , τ], the corresponding domains of integration on the right-hand side with respect to T i and T j are , respectively. Besides, projection of the function G onto the T i T j -plane forms a finite area. Therefore, we regard the right-hand side as a surface integration with surface element vector dr as Here, the notation of S means that we perform the integration on the surface S on the T i T j -plane in the ranges Recall that we assume forward-Euler time integration in explicit simulations. Because all the known information we can use is defined at time t 0 , we use the surface element vector dr 0 as an alternative to dr, as follows: We then convert the right-hand side of equation (23) into a double integration on the T i T j -plane according to the formula for a surface integral of a scalar function G over a surfaceS [19,20], as Each element of the right-hand side in equation (28) represents the gradient of T i or T j on the T i T j -plane at time t 0 . Similarly to the aforementioned case of their time differential, we approximate these gradients as constant. Thus, equation (23) can be rewrittened using equation (27), (28), and a constant value M as follows: By substituting equation (7) and equation (29) into equation (16), we obtain Here, k k = M K. For simplicity, in this paper, we confine k to the case of k > 0. By integrating the left-hand side of equation (31), we obtain Here, ( ) S t 0 indicates the initial entropy at time t 0 . Finally, by performing integration of the second term with respect to T i and double integration of the third term on the right-hand side of equation (32) with respect to T i and T j , we obtain the form of entropy of SDPD at time τ for the general case of N particles: å å T  T  T   T  T  T  T  T  T  T T  T   Because the symmetries of d i and F ij (i=1, 2) can be confirmed from equations (6) and (8) in the case of N=2, the relationships d 1 =d 2 and F 12 =F 21 are true. Hence, the constant parameters of d, F, and α can be introduced as Likewise, the viscous heating fields of f 1 and f 2 become equal. The parameter f can be introduced as  T  T  T  T  T  T  T  T   T  T T  T   ln  Here, C α >0 because k > 0. Figure 1 shows the function of Λ(r) in equation (40). Here, Λ(r) is confirmed to become an increasing function after r/h>1/9, and it is larger than zero within the range of   r h 0 1. Hence, the resulting α becomes positive within   r h 0 1. Although the function Λ(x) becomes discontinuous for  r h 1, the region of  r h 1 is not referred to in the SDPD, so this is not a problem. In the next section, we discuss the characteristics of a two-particle system of SDPD regarding the thermodynamic validity by examining entropy stability conditions.

Derivation of entropy stability conditions
Denote the number of particles as N, the volume of the system as V, and the internal energy as U. Under the conditions of a constant volume process and fixed number of particles (N=V=const.), an entropy stability condition is derived from a thought experiment in which two identical systems transfer an infinitesimal internal energy ΔU between each other. If the entropy is stable, it does not increase regardless of the transfer(namely, S (U−ΔU)+S(U+ΔU)<2S(U)). By taking the Taylor series expansion of both sides of this relationship, we obtain the following stability condition [21]: Let us impose the condition of equation (41) on the entropy S in equation (39). The left-hand side of equation (41), i.e. the second-order partial derivatives of the entropy S by the internal energy U, is calculated by the following processes. First, we derive the first-order partial derivatives of S from the temperatures t T 1 and . 48 The heat capacity C V under an isochoric process is given by Meanwhile, the left-hand side of equation (41) can be rewritten as follows [21]: Here, the fourth and fifth terms vanish because of equation (49). Using equation (42) to equation (50), we obtain the following entropy stability condition:  f a f a a -+ -+ -+ -

Theoretical analyses with case studies
To reproduce a system with a constant temperature gradient, we impose a constraint on the domain of integration as follows: 0 Equation (53) suggests that a linear temperature relationship is imposed before and after the time evolution. By substituting the right-hand side of the two equations in equation (53) into equation (52), we obtain   In the case that the parameter b is sufficiently close to 1 (i.e. the temperature gradient between two particles is moderate), equation (56) works as a stabilizer of the system; it gives the maximum temperature limit as t T 1 increases and reaches the high-temperature area where  ( ) f b C , 0, whereas it gives the minimum temperature limit as t T 1 decreases and reaches the low-temperature area where Meanwhile, when the parameter b is sufficiently close to 0 (i.e. the temperature gradient between two particles is steep), f (b, C) becomes positive everywhere, and the upper part of equation (56) is therefore required as the stability condition.
To deepen the discussion on this issue, we provide a visual representation of the case of By replacing b 2 with the parameter X, we obtain Given that the signs of α and the function F are the same, we can subdivide the entropy stability conditions into eight different types, as listed in table 1. In Type2 and Type8, the system becomes stable because the temperature in Kelvin must be positive, while the sign of a ( ) D f b C , is negative. However, in Type3 and Type5, the entropy becomes unstable everywhere, and no scenario satisfies the condition that the temperature becomes less than or equal to D/α f (b, C) because it is negative.
The main point to emphasize from table 1 is that the entropy stability condition of the system depends on the function F, which indicates that the types of kernel functions influence the entropy stability condition of the system. In the classical SDPD model, because it uses the Lucy kernel, the sign of the function F becomes positive. Hence, the possible types are Type1, Type2, Type5, and Type6.
Let us consider same cases using other kernel functions. In the case that we use the spiky kernel [22] or poly6 kernel [23], each F is positive everywhere in the range of 0<r/h<1. Therefore, the possible types of entropy conditions are the same as those for the Lucy kernel. On the contrary, when we use Mao and Yang's spline kernel [24][25][26], the possible types of entropy conditions are Type3, Type4, Type7, and Type8 because the function F becomes negative in the range of 0<r/h<1. For reference, figure 5 shows the function F in different types of kernel functions.   Table 1. Eight different types of entropy-stability conditions depending on function f (b, C), the parameter D, and the function F. Figure 5. Comparison of function F in different types of kernel functions.
The fact that the kernel function contributes to the entropy stability conditions of a two-particle system can be extended to many-particle systems based on the concept of pair-wise particle methods [27][28][29][30]. In these methods, the total force f acted on a particle is broken down as = å f f ij , where f ij is the force between the ith and jth particles [27]. Namely, the dynamics of a total system can be described as superpositions of two-particle modes.
Consider a case that two persons carry out simulations using the same computational conditions except for their kernel functions. If one uses the spiky kernel and the other uses the spline kernel, the entropy stability conditions imposed on an identical pair of particles differ. In this case, it could happen that only one of the systems is judged to be unstable despite them considering the same physical scenario, which could lead to the emergence of different physical phenomena, such as turbulence or heat transfer. Consequently, the dynamics of the entire multi-particle system could change. How the difference between the judges of entropy stability conditions affects the system must be investigated in future studies.

Conclusion
In this article, we performed an entropy stability analysis of SDPD to evaluate the particle discretization of entropy equations exhibited in SDPD. First, we focused on the simplest SDPD system: a simulation of incompressible flows using an explicit time integration scheme, assuming a quasi-static scenario with constant volume, constant number of particles, and infinitesimal time shift. Next, we derived a form of entropy from the discretized entropy equation of SDPD by integrating with respect to time. Then, we examined the properties of a two-particle system for a constant temperature gradient.
Our theoretical analysis suggests that there exist eight different types of entropy stability conditions, which depend on the types of kernel functions. It was found that the Lucy kernel, poly6 kernel, and spiky kernel produce the same types of entropy stability conditions, whereas the spline kernel produces different types of entropy stability conditions. To summarize, our results suggest that computational parameters of kernel functions contribute to the physical conditions of particle discretization systems.
It is meaningful that we theoretically analyse the two-particle system because such a system with a small number of particles cannot be simulated owing to the lack of numerical accuracy. In that sense, our analysis in this study can be regarded as a 'thought experiment'. Our results contribute to a deeper understanding of particle discretization.
• When f=0, the relationship of β=γ=0 is true. Hence, the value of 4αβγ on the right-hand side of equation (