Dissipative particle dynamics with energy conservation: isoenergetic integration and transport properties

TL;DR

This paper improves the dissipative particle dynamics with energy conservation (DPDE) method by introducing an energy-conserving integration algorithm and analyzing its transport properties, enhancing its accuracy for simulating fluids with temperature gradients.

Contribution

It proposes a local energy-conserving velocity-Verlet algorithm for DPDE and derives analytical expressions for thermal conductivity and fluid transport numbers.

Findings

Energy drift is eliminated with the new algorithm.

Analytical approximation for thermal conductivity is provided.

Transport properties can be tuned to match real fluids.

Abstract

Simulations of nano- to micro-meter scale fluidic systems under thermal gradients require consistent mesoscopic methods accounting for both hydrodynamic interactions and proper transport of energy. One such method is dissipative particle dynamics with energy conservation (DPDE), which has been used for various fluid systems with non-uniform temperature distributions. Despite the success of the method, existing integration algorithms have shown to result in an undesired energy drift, putting into question whether the DPDE method properly captures properties of real fluids. We propose a modification of the velocity-Verlet algorithm with local energy conservation for each DPDE particle, such that the total energy is conserved up to machine precision. Furthermore, transport properties of a DPDE fluid are analyzed in detail. In particular, an analytical approximation for the thermal conductivity coefficient is derived, which allows the selection of a specific value a priori. Finally, we provide approximate expressions for the dimensionless Prandtl and Schmidt numbers, which characterize fluid transport properties and can be adjusted independently by a proper selection of model parameters, and therefore, made comparable with those of real fluids. In conclusion, our results strengthen the DPDE method as a very robust approach for the investigation of mesoscopic systems with temperature inhomogeneities.

Equations

$$\frac{\text{d}{\text{r}}_i}{\text{d}t}={\text{v}}_i,\frac{\text{d}{\text{v}}_i}{\text{d}t}=\frac{1}{m_i}{\text{F}}_i,c_v\frac{\text{d}{T}_i}{\text{d}t}=q_i.$$

$${\text{F}}_i=\sum _{j\ne i}\left({\text{F}}_{ij}^C+{\text{F}}_{ij}^D+{\text{F}}_{ij}^R\right).$$

$${\text{F}}_{ij}^C=a_{ij}{\omega }^C(r_{ij}){\hat{\text{r}}}_{ij},$$

$${\text{F}}_{ij}^D=-{\gamma }_{ij}{\omega }^D(r_{ij})({\text{v}}_{ij}\cdot {\hat{\text{r}}}_{ij}){\hat{\text{r}}}_{ij},$$

$${\text{F}}_{ij}^R={\sigma }_{ij}{\omega }^R(r_{ij}){\xi }_{ij}\Delta t^{-1/2}{\hat{\text{r}}}_{ij}.$$

$${\gamma }_{ij}=\frac{{\sigma }_{ij}^2}{4k_B}\left(\frac{1}{T_i}+\frac{1}{T_j}\right)$$

$$\omega (r_{ij})=\{\begin{array}{ll}{\left(1-\frac{r_{ij}}{r_c}\right)}^{2s}& r_{ij}\le r_c,\\ 0& r_{ij}>r_c,\end{array}$$

$$q_i=\sum _{j\ne i}\left(q_{ij}^{HC}+q_{ij}^{VH}\right),$$

$$q_{ij}^A={\kappa }_{ij}^A{\omega }^H(r_{ij})\left(\frac{1}{T_i}-\frac{1}{T_j}\right)+{\alpha }_{ij}^A{\zeta }_{ij}{\left[{\omega }^H(r_{ij})/\Delta t\right]}^{1/2},$$

$$q_{ij}^B={\kappa }_{ij}^B{\omega }^H(r_{ij})(T_j-T_i)+{\alpha }_{ij}^B{\zeta }_{ij}{\left[{\omega }^H(r_{ij})/\Delta t\right]}^{1/2}.$$

$$\delta P_i$$

$$\delta K_i$$

$$J=\frac{1}{A}\frac{\mathrm{d}Q}{\mathrm{d}t}=\kappa \frac{\mathrm{d}T}{\mathrm{d}x},$$

$$q_{C,ij}^A=\frac{c_v^2}{k_B}{\kappa }_0^A(T_j-T_i){\omega }^H(r_{ij}).$$

$$q_{C,j}^A(h)=\rho {\int }_0^{2\pi }\text{d}\varphi {\int }_0^{{\theta }_c}\text{d}\theta {\int }_h^{r_c}\text{d}{r}_{ij}{r}_{ij}^{2}g(r_{ij})\sin \theta q_{ij}^A,$$

$$q_C^A=\rho A{\int }_0^{r_c}{q}_j^A(h)\text{d}h=\frac{\pi {\rho }^2{c}_v^2{\kappa }_0^A}{k_B}\frac{\mathrm{d}T}{\mathrm{d}x}AH(r_c),$$

$$H(r_c)=$$

$$=$$

$$\Delta ϵ=\overline{ϵ}(x+\lambda )-\overline{ϵ}(x-\lambda )\simeq 2\lambda \frac{\partial \overline{ϵ}}{\partial T}\frac{\mathrm{d}T}{\mathrm{d}x}=2\lambda c_v\frac{\mathrm{d}T}{\mathrm{d}x}.$$

$$⟨{\mathbf{v}}_i(t){\mathbf{v}}_i(0)⟩={\mathrm{e}}^{-t/{\tau }_c}{\mathbf{v}}_i^2(0).$$

$${\tau }_c^{-1}=$$

$$=$$

$${\tau }_c=\frac{3m{k}_{B}T}{2\pi {\sigma }^2\rho r_c^3{I}_2}.$$

$$N_J=\rho A{\int }_0^{\infty }f(v_x)v_x\mathrm{d}{v}_x=\rho A\sqrt{\frac{k_{B}T}{2\pi m}},$$

$${\bar{v}}_x=\frac{{\int }_0^{\infty }\mathrm{d}{v}_x{v}_{x}f(v_x)}{{\int }_0^{\infty }\mathrm{d}{v}_{x}f(v_x)}=\sqrt{\frac{\pi k_{B}T}{2m}}.$$

$$q_D^I=\frac{3(k_{B}T{)}^2{c}_v}{2\pi {\sigma }^2{r}_c^3{I}_2}A\frac{\mathrm{d}T}{\mathrm{d}x}.$$

$$\kappa \simeq \frac{\pi {\rho }^2{c}_v^2{\kappa }_0^A{r}_c^5}{k_B}{I}_1+\frac{3(k_{B}T{)}^2{c}_v}{2\pi {\sigma }^2{r}_c^3{I}_2}.$$

$$\text{Pr}=\frac{\mu C_P}{\kappa },$$

$$\text{Sc}=\frac{\mu }{m\rho D},$$

$$D=\frac{{\tau }_c{k}_{B}T}{m}=\frac{3(k_{B}T{)}^2}{2\pi {\sigma }^2\rho r_c^3{I}_2},$$

$$\mu =$$

$$=$$

$$\equiv$$

$${\text{v}}_i\left(t+\frac{\Delta t}{2}\right)$$

Full text

I Introduction

Simulations of mesoscopically structured systems, ranging from supra-molecular assemblies and artificial self-propelled microswimmers to the flow of biological cell suspensions in complex environments, have become important in studies of a broad variety of biophysical, biomedical, and engineering applications Pivkin et al. (2011); Winkler et al. (2014); Fedosov et al. (2014); Elgeti et al. (2015). This has driven a rapid development of mesoscopic simulation methods employed to advance our understanding of such systems Schiller et al. (2018); Gompper et al. (2009). Mesoscopic methods usually neglect molecular details, but retain features of suspended particles, such as deformability, inter-particle forces, or thermal fluctuations. Examples of mesoscopic methods are lattice Boltzmann method (LBM) McNamara and Zanetti (1988); He and Luo (1997); Succi (2001), multiparticle collision dynamics (MPC) Malevanets and Kapral (1999, 2000); Gompper et al. (2009); Kapral (2008), and dissipative particle dynamics (DPD) Hoogerbrugge and Koelman (1992); Español and Warren (1995, 2017).

Studies of many interesting phenomena require modeling of non-isothermal environments, where temperature gradients are significant, and transport of energy might be relevant. Examples include heat transfer Cleary (1998), thermodiffusion in binary mixtures Sarman and Evans (1992); Wold and Hafskjold (1999); Reith and Müller-Plathe (2000); Artola et al. (2008), and colloidal thermophoresis Galliero and Volz (2008); Volpe et al. (2011); Fedosov et al. (2015); Lüsebrink et al. (2012); Tan et al. (2017). Such problems generally rely on simulation methods which are able to represent the system in a micro-canonical ensemble where energy is exactly conserved. Various non-isothermal systems have been modeled using molecular dynamics Reith and Müller-Plathe (2000); Sellan et al. (2010), thermal LBM McNamara et al. (1997); Shan (1997), MPC Yang and Ripoll (2011); Yang et al. (2014); Yang and Ripoll (2016), and an energy-conserving version of the Monte Carlo method Langenberg and Müller (2016). The original DPD method Hoogerbrugge and Koelman (1992); Español and Warren (1995) is isothermal, and has been extended to account for energy conservation (DPDE) Español (1997); Bonet Avalos and Mackie (1997). DPDE has been successfully applied to simulate a number of thermal gradient problems, such as natural Abu-Nada (2010) and forced convection Yamada et al. (2011), temperature-dependent fluid properties Li et al. (2014), droplet flows Yamada et al. (2016), and thermophoretic Janus colloids Fedosov et al. (2015).

Despite the success of DPDE, current formulation of the method still requires essential improvements due to the difficulties in establishing an efficient and precise integration algorithm. Several algorithms have been employed in DPDE, including the Shardlow splitting scheme, stochastic velocity-Verlet, and velocity-Verlet Shardlow splitting algorithm Bonet Avalos and Mackie (1999); Mackie et al. (1999); Lísal et al. (2011); Larentzos et al. (2014); Li et al. (2014); Homman et al. (2016). Even though these algorithms have been developed as energy-conserving methods, a general problem of all the integration schemes is that they yield a significant net energy drift. This drift can be reduced to a large extent by decreasing the integration timestep. This makes the method less attractive, as it becomes more CPU-time consuming and some remaining drift does not completely disappear, rendering long-time integration difficult. Furthermore, it is not fully clear whether the DPDE method can faithfully represent important features (e.g. thermal diffusivity, momentum transport) of real fluids, as the transport properties of a DPDE fluid are not directly predetermined.

In this paper, we suggest a modification of the well-known velocity-Verlet algorithm, which results in energy conservation up to the order of machine precision, and completely eliminates energy drift. Within this algorithm, changes in kinetic and potential energies are counterbalanced by the internal energy at the level of each DPDE particle, such that the energy remains exactly conserved. We investigate the performance of the algorithm for several choices of conservative interactions between particles and two different formulations of the inter-particle heat conduction. The different heat-conduction terms exhibit nearly identical performance, while different choices for conservative interactions affect the fluid compressibility. Furthermore, an analytical approximation of the overall thermal conductivity $\kappa$ of a DPDE fluid is derived, which allows its direct estimation for selected simulation parameters. The thermal conductivity is governed by two contributions: (i) direct heat conduction and (ii) diffusive heat transport, whose relative ratio can be estimated from the analytical approximation of $\kappa$. Finally, we study the dimensionless Prandtl and Schmidt numbers characterizing fluid transport properties, and show that they can be independently adjusted to match those of real fluids.

The paper is organized as follows. Section II provides details about the DPDE method and the employed integration algorithm, with the details given in Appendix A. In section III, the DPDE method is validated and its performance is studied under thermal gradient conditions. In section III, we also derive the analytical approximation of the fluid thermal conductivity coefficient and discuss transport properties of simulated liquids in comparison to real fluids.

II Method

II.1 Mass and momentum evolution

In the same spirit as in the isothermal DPD method Hoogerbrugge and Koelman (1992); Español and Warren (1995); Marsh et al. (1997), DPDE particles are mesoscopic entities, which represent small fluid volumes containing numerous atoms or molecules. The $N$ constituent particles are characterized by their positions $\textbf{r}_{i}$ and velocities $\textbf{v}_{i}$. Furthermore, each of the DPDE particles is characterized by an additional variable accounting for the internal energy $\epsilon_{i}$ of the small fluid volume. This energy is connected to an internal temperature $T(\epsilon)$ via an entropy function $S(\epsilon)$ as $\partial S(\epsilon)/\partial\epsilon=1/T(\epsilon)$. A simple choice for the entropy function is that of an ideal solid, i.e. $S(\epsilon)=c_{v}\ln(\epsilon)+\mathrm{const}$, where $c_{v}$ is the heat capacity at constant volume. This choice results in a straightforward linear relation between internal energy and temperature as $\epsilon_{i}=c_{v}T_{i}$ Español (1997); Bonet Avalos and Mackie (1997); Ripoll et al. (1998); Qiao and He (2007).

Time evolution of particle characteristics (i.e. position, velocity, and internal energy) is governed by the Newton’s equation of motion and heat equation

[TABLE]

Here, $m_{i}$ is the particle mass, $\textbf{F}_{i}$ is the total force, and $q_{i}$ is the total heat rate. The force $\textbf{F}_{i}$ on particle $i$ is a sum of three pairwise interactions with neighboring particles $j$ given by

[TABLE]

The three contributions correspond to the conservative, dissipative, and random forces which take the form

[TABLE]

Here, $\textbf{r}_{ij}=\textbf{r}_{i}-\textbf{r}_{j}$, $\textbf{v}_{ij}=\textbf{v}_{i}-\textbf{v}_{j}$, $r_{ij}=|\textbf{r}_{ij}|$, $\hat{\textbf{r}}_{ij}=\textbf{r}_{ij}/r_{ij}$, and $\Delta t$ is the integration timestep. The conservative-force coefficient $a_{ij}$ controls fluid compressibility. The coefficients $\gamma_{ij}$ and $\sigma_{ij}$ represent friction and noise amplitudes that are connected to each other through the fluctuation-dissipation relation Español (1997); Bonet Avalos and Mackie (1997); Marsh and Coveney (1998) as

[TABLE]

with the Boltzmann constant $k_{B}$. The random force is determined by $\xi_{ij}$, a symmetric Gaussian random variable (i.e. $\xi_{ij}=\xi_{ji}$) with zero mean and unit variance. Generally, $\sigma_{ij}$ is selected as a constant and $\gamma_{ij}$ is calculated according to Eq. (4). The interaction strengths in Eq. (3) are further controlled by the weight functions which are most commonly chosen as

[TABLE]

where $r_{c}$ is the cutoff radius and $s>0$ is an exponent controlling the interaction strength. While the dissipative and the random functions are linked via the fluctuation dissipation theorem as \omega^{D}=\big{(}\omega^{R}\big{)}^{2}=\omega Español and Warren (1995), the choice of the conservative function $\omega^{C}$ is independent.

II.2 Energy evolution: isoenergetic integration

Time evolution of the internal temperature (or the corresponding internal energy) is governed by the third equality in Eq. (1). The amount of heat that each particle receives (or loses) from its close neighbors per unit time can be expressed as

[TABLE]

where $q_{ij}^{HC}$ corresponds to heat conduction, and $q_{ij}^{VH}$ to viscous heating.

DPDE particles are characterized by an intrinsic temperature, such that every pair of particles exchanges some heat by conduction, proportional to the interparticle temperature difference. Originally Español (1997); Bonet Avalos and Mackie (1997), this term has been proposed to be proportional the difference of the inverse temperatures,

[TABLE]

where $\kappa^{A}_{ij}$ and $\alpha^{A}_{ij}$ are heat-conduction coefficients of the deterministic and random terms, respectively. To guarantee energy conservation, each term in Eq. (7) needs to be antisymmetric under particle interchange. The deterministic term is antisymmetric by construction, and the fluctuating factor $\zeta_{ij}$ is defined as an antisymmetric Gaussian variable with zero mean and unit variance, such that $\zeta_{ij}=-\zeta_{ji}$. Note that $\zeta_{ij}$ is completely uncorrelated with $\xi_{ij}$. Detailed balance imposes the relation between the deterministic and fluctuating terms. Thus, the heat conduction coefficients are connected by the fluctuation-dissipation relation $\left(\alpha^{A}_{ij}\right)^{2}=2k_{B}\kappa^{A}_{ij}$, where $\kappa^{A}_{ij}=\kappa^{A}_{0}c_{v}^{2}(T_{i}+T_{j})^{2}/4k_{B}$ with the constant nominal strength $\kappa^{A}_{0}$. The weight function $\omega^{H}$ can generally be chosen independent of the other weight functions, although we select $\omega^{H}=\omega^{D}$ here for simplicity.

An alternative form for $q_{ij}^{HC}$ was proposed later Ripoll and Ernst (2005a, b), where the heat conduction is directly proportional to the temperature difference between particles as expected from the Fourier law,

[TABLE]

Here, the fluctuation-dissipation relation requires that $\left(\alpha^{B}_{ij}\right)^{2}=\kappa^{B}_{ij}T_{i}T_{j}$ with $\kappa^{B}_{ij}=c_{v}\kappa^{B}_{0}$. The correspondence between these two models for $q_{ij}^{HC}$ can be achieved for $\kappa^{B}_{0}=\kappa^{A}_{0}c_{v}/k_{B}$, if the local temperature differences are not large, i.e. $T_{i}\approx T_{j}$. Quantitative comparison between these two models has not been reported yet and will be presented later.

The second contribution to the heat transport in Eq. (6) is the viscous heating term $q_{i}^{VH}=\sum_{i\neq j}q_{ij}^{VH}$, which represents the work done by conservative and dissipative forces. Thus, it accounts for the variation in mechanical energy as $q_{i}^{VH}=\delta{E_{mec,_{i}}}=\delta P_{i}+\delta K_{i}$, where $P_{i}$ and $K_{i}$ are the potential and kinetic energies whose changes can be calculated as

[TABLE]

The second sum in the variation of kinetic energy is a consequence of the stochastic contribution. By considering the equation of motion (1) and the form of the DPD forces in Eq. (3), the change in mechanical energy can be expressed explicitly through the DPD force parameters Español (1997); Bonet Avalos and Mackie (1997). In fact, such non-trivial expressions for $q_{i}^{VH}$ are often integrated directly using various existing methods such as Shardlow or Trotter schemes Shardlow (2003); De Fabritiis et al. (2006); Serrano et al. (2006). However, we do not use such an expression, but calculate $q_{i}^{VH}$ by a direct tracking of changes in the kinetic and potential energies for every particle at each timestep (see Appendix A for details about the integration algorithm). In this way, viscous heating is accounted for and the energy is conserved locally as well as globally. In a study  Lísal et al. (2011) focused on the performance of different variants of the Shardlow splitting scheme, a related idea has already been outlined, but the energy conservation has been implemented at the level of particle pairs in contrast to the particle level proposed here, and the conservation of potential energy was not considered. A further important advantage of our integration approach is that it allows much easier parallelization than more sophisticated algorithms, such as Shardlow or Trotter schemes Larentzos et al. (2014); Homman et al. (2016), making it a powerful candidate for applications.

II.3 Simulation setup and parameters

Simulation units are selected to be the particle mass $m=1$, the length of the cutoff radius $r_{c}=1$, and the unit of energy $k_{B}T_{0}=1$ with the reference temperature $T_{0}$. In this way, $\tau=r_{c}\sqrt{m/(k_{B}T_{0})}=1$ corresponds to the unit of time. The default average temperature and number density of fluid particles are set to $\bar{T}=T_{0}$ and $\bar{\rho}=3/r_{c}^{3}$, respectively, unless specified otherwise. The DPDE forces employ $\omega^{C}=\omega$ with $s=0.5$ as used commonly in DPD, $\sigma=3k_{B}T_{0}\sqrt{\tau}/r_{c}$, and $\omega^{D}=\omega^{H}=\omega$ with $s=0.25$ whose value has been shown to increase the fluid viscosity in comparison to a frequently used value of $s=1$ Fan et al. (2006); Fedosov et al. (2008). The specific heat is set to $c_{v}=200k_{B}$ and the nominal strengths of interparticle heat conductivities are $\kappa^{A}_{0}=0.001/\tau$ and $\kappa^{B}_{0}=0.2/\tau$ for the two heat-rate models presented above. The dimensionless heat capacity $\beta=c_{v}/k_{B}$ characterizes the number of internal degrees of freedom of a DPDE particle, which in general should be large enough Ripoll and Ernst (2005a).

The conservative-force coefficient $a_{ij}$ is typically taken as a constant in DPDE, following a similar choice in the classical DPD method. However, a dependence of $a_{ij}$ on temperature has been suggested to represent better temperature-dependent fluid properties Li et al. (2014); Groot and Warren (1997). We consider three cases: (i) the ideal gas equation of state case with $a^{C}_{ij}=0$, (ii) the case with a constant coefficient $a^{C}_{ij}=a_{0}$, and (iii) the case with a temperature-dependent coefficient $a^{T}_{ij}=a_{0}(T_{i}+T_{j})/2\bar{T}$, for which the default constant parameter is taken to be $a_{0}=15k_{B}T_{0}/r_{c}$. Different combinations of the heat conduction term (i.e. $q^{A}_{ij}$ and $q^{B}_{ij}$) and the conservative-force coefficient (i.e. $a^{C}_{ij}$ and $a^{T}_{ij}$) employed in simulations will be denoted as e.g. $a^{C}q^{A}$. Note that the default combination is $a^{T}q^{A}$ if not stated otherwise. A default timestep corresponds to $\Delta t=0.005\tau$. Simulations are performed with periodic boundary conditions in a domain with dimensions $L_{x}=20r_{c}$, $L_{y}=L_{z}=10r_{c}$. All other parameters will be specified in the text whenever necessary.

III Results

III.1 Homogeneous temperature

Energy conservation.-

One of the main issues of the DPD method is the establishment of an efficient and consistent algorithm for the integration of the DPD equations, which has been, and still is treated in a large number of related studies Groot and Warren (1997); Nikunen et al. (2003); Thalmanna and Farago (2007); Litvinov et al. (2010); Leimkuhler and Shang (2015); Yamada et al. (2018). Although DPDE was proposed as an energy-conserving method, so far all investigations, in which the viscous heating term $q_{ij}^{VH}$ is calculated explicitly through the system parameters Shardlow (2003); De Fabritiis et al. (2006); Español et al. (2016), have reported different deviations from energy-conserving behavior for various timesteps and integration schemes. For example, Homman et al. Homman et al. (2016) report an energy drift within a range $1-2\times 10^{-5}k_{B}T_{0}$ for $\Delta t=0.006\tau$ and a total integration time of $10\tau$, while Lísal et al. Lísal et al. (2011); Larentzos et al. (2014) find an energy drift on the order of $10^{-4}$, even though the viscous heat rate was not explicitly calculated. Furthermore, Li et al. Li et al. (2014) demonstrate fluctuations in the total energy on the order of $10^{-3}k_{B}T_{0}$, even though no significant energy drift within the total integration time was detected.

As explained in Sec. II.2, the integration algorithm adopted here is fundamentally different from those used previously, since the variation of energy per particle is directly monitored and used for the calculation of the viscous-heating term, avoiding its approximation in terms of parameters such as the relative velocities or positions of all neighboring particles. We have verified the conservation of total energy for a range of timesteps $\Delta t/\tau\in[0.005;0.05]$ and for two system sizes, the one with default values and another one approximately $10$ times larger. Maximum error and fluctuations in the total energy are found to be on the order of $10^{-14}k_{B}T_{0}$, which demonstrates that our DPDE algorithm conserves the total energy by construction (see Appendix A), and therefore, the error in total energy is directly associated with the machine precision.

Note that in order to evaluate the total energy, the mechanical energy needs to be considered together with the internal energy $I_{E}$. The internal energy constitutes the dominant contribution to the total energy, since $I_{E}\approx Nc_{v}T_{0}$ and $c_{v}\gg 1$, and the kinetic energy is $K\simeq 3Nk_{B}T_{0}/2$ (i.e. $I_{E}/K\simeq 2c_{v}/3k_{B}$), while our simulations show that the potential energy is $P\simeq 2K$ for the default conservative interactions.

Temperature definition.-

In the DPDE method, three temperatures can be defined which, for consistency, should be equivalent. The internal temperature $T_{I}$ is computed by averaging over the internal particle temperatures $\langle T_{i}\rangle$; the kinetic temperature $T_{K}$ is calculated by averaging the kinetic energy of all particles as $\langle m_{i}v_{i}^{2}/3\rangle$, and both can also be compared to the reference temperature $T_{0}$. Comparison of the simulation results is shown in Fig. 1a,b for two different timestep values. $\Delta t=0.01\tau$ results in consistent internal and kinetic temperatures, while $\Delta t=0.05\tau$ leads to values of $T_{I}$ and $T_{K}$ which differ by nearly 50%, indicating that this timestep is too large.

Figure 1a also shows that the variance of kinetic temperature is markedly higher than that of the internal temperature, which can be rationalized by differences in temperature distributions. The variance of $T_{K}$ is determined by the Maxwell-Boltzmann distribution of particle velocities and by the number of particles used for $T_{K}$ averaging López-Ruiz and Calbet (2007). The variance of $T_{I}$ is given by $\sigma^{2}_{T_{I}}=\sigma^{2}_{\epsilon}/c_{v}^{2}$, where $\sigma^{2}_{\epsilon}$ is the variance of internal energy, which can be estimated from the distribution of $\epsilon$ in the micro-canonical ensemble Ripoll (2002).

Radial distribution function.-

The measurement of the radial distribution (RDF) function is another test to elucidate the appropriateness of the employed parameters, in particular the timestep. Figure 1c demonstrates that $\Delta t\lesssim 0.02\tau$ leads to timestep-independent RDF within the fluid. Even though Fig. 1 presents results only for the $a^{T}q^{A}$ model, other combinations of the conservative-force coefficient and the heat conduction term have also been tested, resulting in a similar conclusion. Note that the recommended values of the timestep are about one order of magnitude larger than those employed with previous integration methods Li et al. (2014).

III.2 Temperature gradients

Different strategies have been employed to simulate thermal gradients Sarman and Evans (1992); Reith and Müller-Plathe (2000); Lüsebrink and Ripoll (2012), which can be used in DPDE as well. One of these methods, the velocity exchange algorithm Reith and Müller-Plathe (2000), exactly conserves the overall energy by interchanging the velocity of the fastest particle within a pre-determined cold slab with the slowest particle of a hot slab. Alternatively, the average temperature of two pre-defined layers (or regions) can be set to two unequal values Pooley and Yeomans (2005); Yang and Ripoll (2013). Furthermore, simulations with periodic boundary conditions and temperature gradients often employ hot and cold slabs, in which heat is respectively injected and removed, mimicking the contact with reservoirs at different temperatures. In our simulations, heat is injected into the hot slab ($10<x/r_{c}<10.5$) at a constant heat rate $2J_{0}A$, where $J_{0}$ is the heat flux and $A=L_{y}L_{z}$ is the cross-sectional area of the domain, while heat is removed from the cold slab ($0<x/r_{c}<0.5$) at the same rate. Thus, a fixed amount of heat $Q=2J_{0}A\Delta t$ is uniformly added or removed every timestep from the internal energy of all particles in the volume of a given slab. The factor $2$ is due to periodic boundary conditions such that the injected heat $Q$ travels to the left and to the right from the hot slab. Hence, every particle in the hot (cold) slab receives (losses) a heat of $Q/N_{h}$ ($Q/N_{c}$) per timestep, where $N_{h}$ ($N_{c}$) is a time-varying number of particles in the hot (cold) slab. In this way, the total energy in the simulation domain is conserved and the temperature gradient can be indirectly regulated by changing the heat rate or flux.

Equation of state.-

First, we test the development of temperature gradient for the two models of heat conduction term, including $q^{A}$ in Eq. (7) and $q^{B}$ in Eq. (8), together with three different choices of conservative interactions, using $a^{T}$ and $a^{C}$ models with $a_{0}=15k_{B}T_{0}/r_{c}$, and a model with $a_{0}=0$. In the case $a_{0}=0$, the fluid has an ideal-gas equation of state with $p=\rho k_{B}T$, where $p$ is the pressure and $\rho$ is the density of the fluid. For comparison, isothermal DPD with the $a^{C}$ model yields an equation of state where pressure varies quadratically in fluid density Groot and Warren (1997). The six possible model combinations are tested for the same linear temperature profile, which is acquired by properly adjusting $J_{0}$. Thus, the profile $T(x)/T_{0}=0.018x/r_{c}+0.91$ with $T_{min}/T_{0}\approx 0.93$ and $T_{max}/T_{0}\approx 1.073$ is obtained for the models $a^{T}$, $a^{C}$, and $a_{0}=0$, with $J_{0}\tau r_{c}^{2}/(k_{B}T_{0})=2.7$, $2.75$, and $2.8$, respectively. Note that the temperature profile is independent of the choice of the heat conduction term. Figure 2 presents steady-state density profiles for the different models with the same temperature profile, which are nearly linear for all cases. The strongest gradient in fluid density is found for the $a_{0}=0$ model, while the $a^{C}$ model leads to the least variation in $\rho$. This means that the $a^{C}$ model yields the least compressible fluid, while the $a^{T}$ and $a_{0}=0$ models result in more compressible fluids. Note that there is no noticeable difference between the $q^{A}$ and $q^{B}$ models for the heat conduction term in DPDE.

Heat conductivity.-

The Fourier law of heat conduction in one dimension,

[TABLE]

allows the determination of the thermal conductivity $\kappa$ in simulations with fixed $J_{0}$ by analyzing the resulting temperature gradient. Figure 3 shows three temperature profiles for different $J_{0}$. All temperature curves exhibit a nearly linear dependence; however, for the highest $J_{0}$ value, the simulation data slightly deviate from the corresponding linear fit. The deviation of simulated temperature from a linear fit is quantified in the inset of Fig. 3. The temperature profile for $J_{0}\tau r_{c}^{2}/(k_{B}T_{0})=10$ can be fitted well by a quadratic function. Nevertheless, the deviation of temperature from the linear fit is small and remains within $\leavevmode\nobreak\ 2\%$. A large enough heat flux clearly leads to a strong temperature gradient with local changes in fluid density and structure characterized by RDF. This in turn may affect local heat conductivity, as it depends on inter-particle distances, resulting in a slightly non-linear temperature profile.

Figure 4a presents the dependence of heat conductivity $\kappa$ on the average fluid temperature and density for the reference parameters, using a small value of $J_{0}\tau r_{c}^{2}/(k_{B}T_{0})=0.25$. Thermal conductivity increases when average fluid density or temperature increase, which is qualitatively consistent with other experimental and theoretical investigations Chapman et al. (1990); Guildner (1975); Hanley et al. (1974); Ziebland and Burton (1955). We will discuss these dependencies in more detail in Sec. III.4, where an analytical expression for $\kappa$ is derived. Figure 4b shows the dependence of thermal conductivity $\kappa$ on heat flux $J_{0}$ for different models of the heat conduction term and conservative interactions. Since $\kappa$ depends on local temperature and density, it is obtained here by fitting temperature profiles within a small region where particle densities and temperatures are very close to their average values. Differences in $\kappa$ values for the $q^{A}$ and $q^{B}$ conductivity terms are not really significant, while a temperature-dependent conservative-force coefficient $a^{T}$ results in slightly smaller heat conductivities in comparison to the constant conservative force with $a^{C}$. Heat conductivity exhibits a slight decrease with increasing input heat flux $J_{0}$. However, this decrease is within $3-4\%$, indicating that $\kappa$ in the DPDE method is rather robust and nearly independent of the heat flux.

The fact that the DPDE model has a finite value of heat conductivity implies that only a limited amount of heat can be transferred, which determines a maximum heat flux $J_{max}$. This means that if lager heat fluxes are applied, some parts of the system may attain a negative internal temperature. This is clearly incorrect, as noted already in Ref. Ripoll et al. (1998), and leads to the simulation instability. The maximum heat flux $J_{max}$ depends very strongly on $\kappa_{0}$, especially in the cases where the energy transport by heat conduction dominates over the diffusive transport, and indirectly also on other system parameters, which are more important in cases when the diffusive transport is high (e.g. for low $\kappa_{0}$ values and moderate temperature gradients). The value of $J_{max}$ appears to be slightly larger for simulations with $q^{A}$.

III.3 Energy transfer

To investigate the mechanisms for heat transfer within a DPDE fluid in more detail, we define a reference plane to monitor the involved energy fluxes. This plane is perpendicular to the gradient direction, and can be placed at different positions between the hot and cold slabs. Given a heat flux $J_{0}$, energy conservation guarantees that the total energy flux through a plane at any position is constant, $J_{tot}=J_{0}$. Four different energy fluxes can be distinguished from two different types. One type is a heat-conduction flux, $J_{C}$, which corresponds to the exchange of internal energy between fluid particles located at different sides of the plane within distances smaller than the cutoff radius $r_{c}$, according to the $q^{HC}$ term in Eq. (7) or Eq. (8). Another contribution is the diffusive flux $J_{D}$ of energy, realized by fluid particles which actually cross the plane. Mass conservation enforces that the mass flux through the plane is on average zero, but due to the externally imposed temperature gradient, particles on the hot side are more energetic than on the cold side which results in a diffusive transfer of energy. In this way, diffusive energy transfer includes potential, kinetic, and internal contributions.

Figure 5 shows various contributions to the total energy flux through a plane in the hot area at $x=8r_{c}$. The major contributions correspond to the fluxes of particle internal energy, $J_{C}$ and $J_{D}^{I}$, while the diffusive fluxes from the potential and kinetic energies, $J^{P}$ and $J^{K}$, are very small. This can be intuitively understood, since the ratio between the particle internal energy and potential/kinetic energy is proportional to $c_{v}$, which is $c_{v}/k_{B}=200$ here. Note that in this example, the diffusive flux of internal energy $J_{D}^{I}$ is larger than the conductive flux $J_{C}$, and that $J_{tot}=J_{C}+J_{D}$, with $J_{D}=J_{D}^{I}+J^{P}+J^{K}$. Fluctuations of the total flux around the input value are due to the statistical error of the measurement procedure.

Figure 6a shows the conductive flux $J_{C}$ for various models through a plane in the hot area at $x/r_{c}=8$ as a function of the externally imposed $J_{0}$. The differences between the two implementations of the heat conduction term, $q^{A}$ and $q^{B}$, are systematic, but very small. However, the conductive flux decays faster for the $a^{T}$ model of the conservative-force coefficient in comparison to $a^{C}$ with increasing $J_{0}$. Note that for these two models, the temperature at the reference plane is the same, but the density is different, as can be seen in Fig. 2. This monotonic decay of the conductive flux with increasing $J_{0}$ does not occur though, when the reference plane is placed in the cold area at $x/r_{c}=2$, as can be seen by the increase of $J_{C}$ (open blue symbols) in Fig. 6b. Correspondingly, the diffusive flux $J_{D}$ increases with the injected heat through the plane at the hot area, while it decreases through the plane at the cold area. Therefore, local temperature strongly affects not only the overall heat conductivity as shown in Fig. 4, but also the ratio between conductive and diffusive fluxes of energy. In Fig. 6b, this means that on the hot side, the energy flux is mainly diffusive, while on the cold side, it is mainly conductive (for not too small applied heat fluxes).

III.4 Analytical calculation of energy transfer

In order to provide an analytical expression for the heat conduction coefficient $\kappa$, we analyze both the conductive and diffusive fluxes of the internal energy to total energy transfer, since kinetic and potential energies have a much smaller contribution (see Fig. 5). The conductive contribution to the energy transfer between two particles $i$ and $j$ at opposite sides of the reference plane is given by Eq. (7) or Eq. (8). These two expressions have been demonstrated to be quantitatively very similar, which means that the assumption $(T_{i}+T_{j})^{2}/(T_{i}T_{j})\simeq 4$ is valid for reasonably small temperature differences. The random part of heat conduction is zero on average, such that the average heat rate from Eq. (7) between two particles at the opposite sides of the reference plane becomes

[TABLE]

The linear temperature profile can be approximated as $T_{j}-T_{i}=\dfrac{\mathrm{d}T}{\mathrm{d}x}\Delta x=\dfrac{\mathrm{d}T}{\mathrm{d}x}r_{ij}\cos\theta$, where $\theta$ is the angle between $r_{ij}$ and the temperature gradient axis, see Fig. 7. Considering that particle $j$ is located at a distance $h$ from the plane, the total heat conduction rate between particle $j$ and its neighbors on the other side of the reference plane is given by

[TABLE]

with $\theta_{c}=\arccos(h/r_{c})$. For a non-ideal fluid, the RDF is $g(r_{ij})\neq 1$ and needs to be pre-computed for the numerical integration of Eq. (12). Total heat rate by conduction is then calculated by integrating $q^{A}_{C,j}(h)$ over a volume $Ar_{c}$ as

[TABLE]

where $A$ is the area of the reference plane and

[TABLE]

The second equality above makes use of the change of variables $s=h/r_{c}$ and $l=r/r_{c}$, while the last equality defines a numerical coefficient $I_{1}$, which in general depends on $g(r)$ and $\omega^{H}(r)$, and therefore on the system parameters. With the employed default parameters in Sec. II.3, $I_{1}=0.058$.

The second contribution to heat transfer by particle diffusion is calculated by considering particles crossing the reference plane at $x=$const. It is equal to $q^{I}_{D}=N_{J}\Delta\epsilon$, where $N_{J}$ is the number of particles crossing the reference plane from hot to cold per unit time. The internal energy transferred by each particle crossing from the hot to the cold side can be calculated as

[TABLE]

Here, the mean free path $\lambda$ can be approximated as $\lambda=\tau_{c}\bar{v}_{x}$, where $\tau_{c}$ is the collision time and $\bar{v}_{x}$ is the average velocity of particles crossing the plane in one direction. The collision time is related to the decay of the velocity autocorrelation function, which can be calculated in simulations by assuming an exponential decay of the velocity autocorrelation as

[TABLE]

In order to provide an analytical estimate for $\tau_{c}$, we follow the procedure in Ref. Groot and Warren (1997), where the friction from dissipative interactions defines $1/\tau_{c}$ and the sum over dissipative forces for neighboring particles is replaced by an integral. The random force vanishes on average and the conservative force determines the RDF, such that the expression in Ref. Groot and Warren (1997) can be generalized to

[TABLE]

The last equality defines the numerical coefficient $I_{2}$, which can be calculated for the default parameters in Sec. II.3 to be $I_{2}=0.131$. Substitution of $\gamma=\sigma^{2}/(2k_{B}T)$ (i.e. for $T_{i}\simeq T_{j}$) into Eq. (III.4) results in

[TABLE]

A good agreement between simulated collision times and the analytical approximation is shown in Fig. 8 for different values of the average temperature and density values.

Assuming that the mean free path is small, that the temperature gradients are not large, and that the velocity distribution function $f(v)$ can be approximated by the Maxwell-Boltzmann distribution, we can calculate both $N_{J}$ and $\bar{v}_{x}$. The number of particles, $N_{J}$, crossing the reference plane from hot to cold per unit time, accounts only for particles with positive velocities,

[TABLE]

and the average velocity of those particles is

[TABLE]

The results from Eqs. (18), (19), and (20) together yield the diffusive rate of internal energy

[TABLE]

The total heat rate is then the sum of the contributions from Eqs. (13) and (21). With Eq. (10), the heat conductivity can then be approximated as

[TABLE]

This analytical expression can be compared with the simulation results shown in Fig. 4, where a very satisfactory agreement is obtained without any adjustable parameters.

III.5 Prandtl and Schmidt numbers

A quantitative comparison of the properties of any simulated fluid with those of real fluids can be achieved by considering dimensionless numbers, such as the Prandtl number, $Pr$, or the Schmidt number, $Sc$. Prandtl number is the ratio of momentum and energy transport defined as

[TABLE]

where $\mu$ and $C_{P}$ are dynamic viscosity and specific heat capacity at constant pressure, respectively. For the large values of $c_{v}$ in DPDE, $C_{P}\simeq c_{v}$. The Schmidt number $Sc$ is the ratio of momentum and diffusive mass transport defined as

[TABLE]

where $D$ is the translational diffusion coefficient.

The Prandtl and Schmidt numbers in DPDE are functions of model parameters, and can be set to various values independently. The transport properties of a DPD fluid have been investigated in several studies Groot and Warren (1997); Fan et al. (2006); Noguchi and Gompper (2008). Expressions for $D$ and $\mu$ have been derived, and can be generalized in a similar way as previously shown for the collision time. This yields

[TABLE]

[TABLE]

Eqution (III.5) defines the numerical coefficient $I_{3}$, which can be calculated for the default parameters to be $I_{3}=0.073$, resulting in $Pr=2.0$ and $Sc=3.57$. Alternatively, we obtain slightly larger values, $Pr=2.1$ and $Sc=4.1$, by keeping all default parameters but considering the ideal case with $g(r)=1$. Furthermore, different expressions for the weight functions of the conservative, dissipative, and heat terms will lead to different $Pr$ and $Sc$ values. Liquids have typical $Pr$ and $Sc$ values within ranges $5-15$ and $100-1000$, respectively Lide (2004), while for gases, both $Pr$ and $Sc$ are typically close to unity Lide (2004); Kestin et al. (1984). Thus, our default system lies at the boundary between gaseous and fluid behavior. Using the expressions for transport coefficients above, model parameters (e.g. $\rho$, $r_{c}$, $\sigma$, $c_{v}$, $\kappa_{0}$) of a DPDE fluid can be selected such that the Prandtl and Schmidt numbers correspond more accurately to those of any given fluid. For example, an increase in the fluid density $\rho$ and/or the cutoff radius $r_{c}$ will lead to Prandtl and Schmidt numbers which are larger than those employed here, but at the same time will result in an increase of computational cost.

IV Summary and conclusions

In this paper, various aspects related to the performance of the DPDE method have been systematically investigated. A modification of the velocity-Verlet integration algorithm is suggested, which has a simple implementation and conserves local and total energies up to the order of machine precision. In this algorithm, the changes in local kinetic and potential energies are exactly counterbalanced by the internal energy every timestep at the level of single DPDE particles. The model is validated by verifying the equivalence of the two alternative definitions of the temperature (kinetic and internal), $T_{K}=T_{I}$, and by studying the behavior of particle RDF, which suggests that the timestep should be $\Delta t\lesssim 0.02\tau$ (here, $\tau$ is the unit of time) for a stable and consistent simulation.

Several choices for the heat-conduction term and conservative interactions in the DPDE method have been studied and compared. The interparticle heat conduction in DPDE has been postulated in two different ways. The expression $q^{B}$ is more intuitive since it is directly proportional to the particles’ temperature difference as expected from the Fourier law, while $q^{A}$ is proportional to the particles’ inverse temperature difference, which has been commonly used so far. Here, we show that when the prefactors in both approaches are properly related, the differences between the predictions of both models are minimal. Compressibility effects are especially relevant in the presence of a temperature gradient, due to the related differences in the density distribution, such that we investigate different choices for the conservative interactions. A DPDE fluid with temperature-dependent conservative interactions shows a compressibility lower than a DPDE fluid with an ideal-gas equation of state, but clearly larger than a DPDE fluid with constant conservative interactions. Nevertheless, the heat conduction coefficients are similar for both temperature-dependent and constant conservative interactions.

Thermal conductivity of the DPDE models is measured in simulations by fitting the Fourier law to an induced temperature gradient. Thermal conductivity describes the system in a very robust manner, showing minor changes for a wide range of temperature gradients. Energy in a DPDE fluid is transferred by the diffusive motion of DPDE particles, and by heat conduction due to the internal particle temperature, which depends on local average temperature. The inter-particle conductivity is directly controlled by $\kappa_{0}$, while the diffusive transport depends on the particle friction (or equivalently $\sigma^{2}$), such that depending on local conditions, conductive or diffusive transport of heat may dominate. Good agreement between the analytical expression and simulation measurements of the thermal conductivity is obtained for a large range of parameters. We have also presented analytical approaches to obtain the most important transport coefficients, and therefore the two most relevant non-dimensional fluid numbers, the Prandtl ($Pr$) and Schmidt ($Sc$) numbers. These expressions allow the selection of simulation parameters such that the corresponding $Pr$ and $Sc$ of a DPDE fluid approximate well those of liquids. In conclusion, our results provide a detailed guidance on how to properly employ the DPDE method in simulations of various mesoscopic systems with temperature inhomogeneities.

Appendix A Integration algorithm

For integration of Eq. (1), the velocity-Verlet algorithm Allen and Tildesley (1991) is adapted. In the first integration step, both particle velocities and internal energies are advanced half timestep, while particle positions are integrated full timestep as follows

[TABLE]

Note that the term $[\cdots]$ reflects a change in the particle’s kinetic energy $K_{i}(t)=\dfrac{m_{i}}{2}v_{i}^{2}(t)$, which is counterbalanced by a portion of internal energy in order to conserve the total energy locally. Before the second integration step, particle forces $\textbf{F}_{i}(t+\Delta t)$, potential energies $P_{i}(t+\Delta t)$, and energy rates $q_{i}(t+\Delta t)$ are computed based on $\textbf{r}_{i}(t+\Delta t)$ and $\textbf{v}_{i}\left(t+\Delta t/2\right)$. In the second step, particle velocities and internal energies are integrated half timestep as

[TABLE]

Here, changes in both kinetic and potential energies are offset against the internal energy to satisfy local energy conservation. Potential energy corresponds to the conservative interaction as $\mathbf{F}_{ij}^{C}=-\vec{\nabla}P_{ij}$, such that $P_{i}=\sum_{j}\dfrac{a_{ij}r_{c}}{2}\left(1-\dfrac{r_{ij}}{r_{c}}\right)^{2}$.

Acknowledgments

The authors gratefully acknowledge the computing time granted through JARA-HPC on the supercomputer JURECA Jülich Supercomputing Centre (2018) at Forschungszentrum Jülich.