Temperature-dependent properties of liquid-vapour coexistence system with many-body dissipative particle dynamics with energy conservation

Kaixuan Zhang School of Aerospace Engineering and Applied Mechanics, Tongji University, Shanghai 200092, China Jie Li School of Aerospace Engineering and Applied Mechanics, Tongji University, Shanghai 200092, China Shuo Chen schen_tju@mail.tongji.edu.cn School of Aerospace Engineering and Applied Mechanics, Tongji University, Shanghai 200092, China Yang Liu Department of Mechanical Engineering, The Hong Kong Polytechnic University, Hong Kong, China

Abstract

The dynamic properties of fluid, including density, surface tension, diffusivity and viscosity, are temperature-dependent and can significantly influence the flow dynamics of mesoscopic non-isothermal systems. To capture the correct temperature-dependence of a fluid, a many-body dissipative particle dynamics model with energy conservation (mDPDe) is developed by combining the temperature-dependent coefficient of the conservative force and weighting terms of the dissipative and random forces. The momentum and thermal diffusivity, viscosity, and surface tension of liquid water at various temperatures ranging from 273 K to 373 K are used as examples for verifying the proposed model. Simulations of a periodic Poiseuille flow driven by body forces and heat sources are carried out to validate the diffusivity of the present model. Also, a steady case of heat conduction for reproducing the Fourier law is used to validate the thermal boundary conditions. By using this mDPDe simulations, the thermocapillary motion of liquid water nanodroplets on hydrophobic substrates with thermal gradients is investigated. The migration of the droplet is observed on flat substrates with gradient temperature. The velocity of the migration becomes larger for higher temperature difference, which is in agreement with the present theoretical analysis and DVDWT simulations. The results illustrate that the modified model can be used to study Marangoni effect on a nanodroplet and other heat and mass transfer problems with free interface.

I introduction

Dissipative particle dynamics is a particle-based model which was proposed by Hoogerbrugge and Koelman 1992Hoogerbrugge and has been developed a suitable method for micro- and nano simulations. By combining the advantages of large timescale in lattice-gas automata (LGA), and the mesh-free property in molecular dynamics (MD), DPD is more efficient for simulating larger fluid system than MD. Moreover, DPD is defined as a coarse-grained model, which means every DPD particle represents a cluster of atoms/molecules and the computation is greatly reduced. Due to the larger spacial and time scale, DPD becomes a flexible method for investigate simple fluids 2006_Chen and complex fluid problems such as polymer 2009_Li and DNA suspensions 2006_Fan , red blood cells dynamics 2013_Peng ; Ye2014Dissipative , and biofluids Xu2011Dissipative .

The interaction between DPD particles determines the basic properties of fluid system, which involving conservative force, dissipative force and random force. The conservative force is a soft repulsive force, which makes contribution to the fluid compressibility. The dissipative force generates friction between DPD particles, which can describe the viscosity of the system. The random force makes up the defect due to coarse-graining treatment by introducing a stochastic force on each DPD particle and meanwhile, these two forces satisfy the fluctuation-dissipation theorem and act as a thermostat to control the system at a constant temperature 1997Groot . Since there is no attractive force between DPD particles, the equation of state does not match van der Waals curve and can not simulate vapour-liquid coexistence system. By introducing attractive force term and adding local density in the repulsive term, the modified conservative force produces a $many-body$ interaction with $long~{}range$ attraction and $short~{}range$ repulsion 2003_Warren . The pressure of the system becomes a cubic relationship with density, which can generate a free interface. This many-body dissipative particle dynamics method has been widely used to simulate droplet dynamics on substrates, multiphase flow in complex geometry and phase transition problems 2011_Arienti ; 2012_Li_Dissipative ; Mo2017Hypernetted ; Zhao2018Rheology ; 2018_Pan ; Yuxiang2019Droplet ; 2019_kai .

Both the classic DPD and its $many-body$ version, i.e. mDPD, can only describe isothermal system due to the thermostat in the models. To simulate the process of heat transfer, an extra temperature quantity is introduced on each DPD particle and generate energy-conserved dissipative particle dynamics (eDPD) Espa_ol_1997 , which has been adopted to investigate heat conduction 2007_Qiao , natural convection Abu2010Natural , and solidification 2016Modeling . Ripoll $et~{}al.$ 1998_Ripoll_Dissipative confirmed that the heat conduction of eDPD system obeys Fourier law by simulating the 1D heat conduction problem. He $et~{}al.$ He2008Self adopted eDPD to investigated the process of heat conduction between nano-fluid and nano-materials. Abu-Nada $et~{}al.$ Abu2010Natural considered Neumann and Dirichlet boundary conditions and simulate the 2D heat conduction problem by using eDPD method. They also modified the parameters in eDPD model to study natural convection and then validated their results by comparing with CFD simulations. Li $et~{}al.$ Li2014Energy proposed analytical formula for determining the mesoscopic heat friction and validated the prediction by reproducing the experimental data for Prandtl number of liquid water at various temperatures. Johansson $et~{}al.$ 2016Modeling adopted eDPD to simulate phase transfer problem, and they captured the process of solidification which is agreement with theoretical analysis. They also suggested that a modified model is needed to reproduce the release of the latent heat during the process of the transition from liquid to solid phase.

To extend the isothermal mDPD equations to modeling heat transport in non-isothermal multiphase fluid systems, many-body dissipative particle dynamics with energy conservation (mDPDe) was developed by combining mDPD with eDPD. Yamada $et~{}al.$ Yamada2016Dissipative simulated the process of heat conduction between droplets and solid substrates by using mDPDe. The droplet density at different temperature seems to be constant in their simulations, which does not match that of most liquids in the real world. Wang $et~{}al.$ 2020_Wang_zhang adopted mDPDe to investigate ice crystal nucleation leading to droplet freezing on flat substrates. They captured the nucleation and shape deformation process of a water droplet during freezing which are in agreement with the experimental data. Their work illustrated the modified mDPDe can be used to study freezing of water droplet. However, the corresponding temperature-dependent properties of the system at various temperature are still not clear in the numerical model, which are also important to illustrate the validity of the model.

The objective of the present work is to propose a model for capturing the correct temperature-dependent properties of fluid system with free interface, such as Marangoni effect on a droplet. Specifically, liquid water is used as an example for verifying the mDPDe model we propose. The density, surface tension, diffusivity and viscosity of liquid water as well as its Schmidt and Prandtl numbers in the range of 273K to 373K are reproduced with the present mDPDe model. The relationship between the surface tension and the expressions of the conservative force have been analyzed. The temperature-dependent density and surface tension are obtained, which can be important for investigating fluid problems with complex mass transfer. The relation of the weight function with temperature proposed by Li $et~{}al.$ is adopted Li2014Energy . Results for the diffusivity and viscosity as well as Schmidt number at various temperatures are presented and compared with the available experimental data of liquid water. Furthermore, a fitting formula for correcting the mesoscopic heat friction is obtained, which reproduces more accurate Prandtl number of liquid water at various temperatures. This proposed model is not limited to liquid water, and it can be readily extended to other fluids for modelling the correct dependence of thermal hydrodynamic properties in terms of temperature.

The paper is organized as follows. In Section 2 we describe details of mDPDe formulations and parameters, including analytical formulas for determining the surface tension and the mesoscopic heat friction. Section 3 presents our validation of the mDPDe model, and the performance of the present mDPDe model in reproducing correctly the temperature-dependent properties consistent with the experimental data. Furthermore, we simulate the thermalcapillary motion of a droplet on a hydrophobic substrate with temperature gradients and compare the velocity with theoretical analysis and DVDWT simulations. Finally, we conclude with a brief summary in Section 4.

II Method

Many-body dissipative particle dynamics model with energy conservation(mDPDe) is developed based on many-body dissipative particle dynamics (mDPD) and energy conservation dissipative particle dynamics(eDPD). Each mDPDe particle interacts with other particles through distance- and velocity-dependent forces and energies within a certain cutoff radii 1992Hoogerbrugge ; Espa_ol_1997 ; J1997Dissipative . The momentum and energy transfer between mDPDe particles obey the momentum and energy conservation, respectively. The evolution of the particle motion is controlled by Newton’s second law:

m_{i}\frac{d\vec{v}_{i}}{dt}=\sum_{j\neq i}\left(\vec{f}_{ij}^{C}+\vec{f}_{ij}^{D}+\vec{f}_{ij}^{R}\right)+\vec{f}_{ext}

(1)

where $\vec{f}_{ext}$ is the external force. $\vec{f}_{ij}^{C}$ , $\vec{f}_{ij}^{D}$ , and $\vec{f}_{ij}^{R}$ denote the conservative, dissipative and random forces among particles, respectively. Their specific expressions are as follows:

\vec{f}_{ij}^{C}=A_{ij}\omega^{C}\left(r_{ij}\right)\vec{e}_{ij}+B_{ij}\left(\tilde{\rho}_{i}+\tilde{\rho}_{j}\right)\omega^{d}\left(r_{ij}\right)\vec{e}_{ij}

(2)

\vec{f}_{ij}^{D}=-\gamma_{ij}\omega^{D}\left(r_{ij}\right)\left(\vec{e}_{ij}\cdot\vec{v}_{ij}\right)\vec{e}_{ij}

(3)

\vec{f}_{ij}^{R}=\sigma_{ij}\omega^{R}\left(r_{ij}\right)\zeta_{ij}\Delta t^{-1/2}\vec{e}_{ij}

(4)

where $A_{ij}$ and $B_{ij}$ are the amplitudes of attractive and repulsive forces, respectively. $\omega^{C}\left(r_{ij}\right)=\left(1-r_{ij}/r_{C}\right)$ and $\omega^{d}\left(r_{ij}\right)=\left(1-r_{ij}/r_{d}\right)$ represent the distance-dependent weight functions which vanish if the distance is larger than the corresponding cut-off radii, $r_{C}$ and $r_{d}$ , respectively. $\tilde{\rho_{i}}$ and $\tilde{\rho_{j}}$ are the weighted local density functions of particles $i$ and $j$ , respectively, which are described as $\tilde{\rho}_{i}=\sum_{j\neq i}\omega^{\rho}\left(r_{ij}\right)$ . And in the description, the weight function of the local density is expressed as $\omega^{\rho}\left(r_{ij}\right)={105}/{16\pi r_{d}^{3}}\cdot\left(1+{3r_{ij}}/{r_{d}}\right)\left(1-{r_{ij}}/{r_{d}}\right)^{3}$ , and $\int_{0}^{\infty}d^{3}\mathbf{r}\omega^{\rho}\left(r_{ij}\right)=1$ as the weight function is normalized by the factor ${105}/{16\pi r_{d}^{3}}$ . Another choice is $\omega^{\rho}\left(r_{ij}\right)={15}/{2\pi r_{d}^{3}}\cdot\left(1-{r_{ij}}/{r_{d}}\right)^{2}$ , which can provide smaller density for the system. In this work, we adopt the former one. $\zeta_{ij}=\zeta_{ji}$ is a random number with zero mean value and a variance of unity, which keeps the momentum of the interacting pair of particles to be conserved. The dissipative coefficient $\gamma_{ij}$ and random coefficient $\sigma_{ij}$ satisfy the fluctuation-dissipation theorem, which requires:

\sigma_{ij}^{2}=\frac{4\gamma_{ij}k_{B}T_{i}T_{j}}{\left(T_{i}+T_{j}\right)},\quad\omega^{R}\left(r_{ij}\right)=\left[\omega^{D}\left(r_{ij}\right)\right]^{1/2}

(5)

where $k_{B}$ and $T$ are Boltzmann constant and system equilibrium temperature, respectively. $\omega^{D}\left(r_{ij}\right)$ and $\omega^{R}\left(r_{ij}\right)$ are the weight functions for the dissipative and random forces, respectively, in which $\omega^{D}\left(r_{ij}\right)=\left(1-r_{ij}/r_{C}\right)^{s_{v}}$ for $r_{ij}<r_{C}$ . $s_{v}$ is related with the viscosity of fluid and here, we adopt the modified expression proposed by Li $et~{}al.$ Li2014Energy to obtain temperature-dependent viscosity.

As the method of MDPDe satisfies the energy conservation equation, the heat transfer between particles is realized by the exchange of internal energy, which is an additional property and expressed as follows Espa_ol_1997 :

C_{v}\frac{dT_{i}}{dt}=\sum_{j\neq i}\left(q_{ij}^{V_{h}}+q_{ij}^{C_{h}}+q_{ij}^{R_{h}}\right)

(6)

where the heat flux, $q_{i}$ , is the sum of collision-caused heat flux $q_{ij}^{C}$ , mechanical energy-caused viscous heating $q_{ij}^{V}$ and thermal fluctuation-caused heat flow $q_{ij}^{R}$ within the truncation radius. $C_{V}$ is the heat capacity at constant volume, which is normalized by $k_{B}$ . The specific expressions are given by

q_{ij}^{V_{h}}=\frac{1}{2C_{v}}\left[\omega^{D}\left(r_{ij}\right)\left\{\gamma_{ij}\left(\vec{e}_{ij}\cdot\vec{v}_{ij}\right)^{2}-\frac{\sigma_{ij}^{2}}{m_{i}}\right\}-\sigma_{ij}\omega^{R}\left(r_{ij}\right)\left(\vec{e}_{ij}\cdot\vec{v}_{ij}\right)\zeta_{ij}\right]

(7)

q_{ij}^{C_{h}}=\kappa_{ij}\omega^{CT}\left(r_{ij}\right)\left(\frac{1}{T_{i}}-\frac{1}{T_{j}}\right)

(8)

q_{ij}^{R_{h}}=\alpha_{ij}\omega^{RT}\left(r_{ij}\right)\zeta_{ij}^{e}\Delta t^{-1/2}

(9)

where $k_{ij}$ and $\alpha_{ij}$ represents the strengths of conductive and random heat fluxes respectively, which are given by

\kappa_{ij}=k_{o}k_{B}T_{eq}^{2}\left(\frac{\epsilon_{i}+\epsilon_{j}}{2k_{B}T_{eq}}\right)^{n_{k}}

(10)

\alpha_{ij}=\sqrt{2k_{B}\kappa_{ij}}

(11)

where $\epsilon_{i}$ is the interal energy. $k_{o}$ is a positive constant which determines the thermal conductivity of mDPDe particles. $T_{eq}$ is the equilibrium temperature of the system and $n_{\kappa}$ is the constant which was chosen as 2 in the previous work 2007_Qiao ; Abu2010Natural . $\kappa$ is a constant related to the thermal conductivity of mDPDe particals, and $\omega^{CT}\left({r_{ij}}\right)$ and $\omega^{RT}\left(r_{ij}\right)$ are the weight functions, respectively and their relation can be described as $\omega^{CT}\left({r_{ij}}\right)=\left[\omega^{RT}\left(r_{ij}\right)\right]^{2}$ and $\omega^{CT}\left({r_{ij}}\right)=\left(1-r_{ij}/r_{CT}\right)^{s_{T}}$ . $s_{T}$ is assigned as 2 in the previous paper Li2014Energy . $\zeta_{ij}^{e}=-\zeta_{ji}^{e}$ is the random number having similar properties as $\zeta_{ij}$ , which also keeps the energy of the interacting pair of particles to be conserved Espa_ol_1997 .

III Results and Discussion

III.1 Validation of Fourier law

To validate the hydrodynamic property of the present non-isothermal mDPDe model, we simulate periodic Poiseuille flow with a simple mDPDe fluid and compare the results with isothermal mDPD model. Two equal force in reverse direction $\left(F=0.02\right)$ are applied on each mDPDe particle to drive the flow. The mDPDe parameters for the periodic Poiseuille flow are listed in the caption of Fig. 1. The computational domain is divided into 50 bins along the $z$ -direction. The equilibrium velocity profiles is shown in Fig. 1, which are obtained by averaging enough sampled data. The results are in agreement with that of mDPD model, which illustrates that we can use the modified mDPDe model to obtain the same hydrodynamic property with mDPD.

Refer to caption — Figure 1: The velocity profile of the periodic Poiseuille flow along the $z$ -direction with the attractive parameter $A=-40$ , the repulsive parameter $B=25$ , the temperature $k_{B}T=1.0$ , the dissipative parameter $\gamma=8.0$ , and the external force $F=0.02$ .

The next test case for validation is heat conduction between a cold wall and a hot wall. The schematic of the geometry is shown in Fig. 2, in which the stationary fluid is confined between a hot wall of $T_{H}=1.2$ and a cold wall of $T_{C}=1.0$ . We test different parameters, including constant attractive parameter $A=-40$ and repulsive parameter $B=25$ ; constant $A=-40$ and temperature-dependent $B=25*k_{B}T$ ; and temperature-dependent $A=-40*k_{B}T$ and constant $B=25$ . Fig. 3 (a) shows the temperature profile of the fluid from all these tests has a linear spatial distribution for the steady-state heat conduction, which obeys the Fourier law. However, as shown in Fig 3 (b), the density in these systems change differently with the increasing temperature in the system. Referring to the real density profile of water, the density decreases with increasing temperature in most range from 273 K to 373 K, a temperature-dependent repulsive parameter $B=25*k_{B}T$ and a constant value of attractive parameter $A=-40$ is more proper for simulating the heat conduction between liquid water and substrates.

III.2 Temperature-dependent properties

The above two cases validate the modified mDPDe model. For an isothermal fluid system, the present model captures the correct behaviours of fluid as mDPD does. For a non-isothermal fluid system, the present model can reproduce Fourier law in heat conduction between a hot wall and a cold wall. Our next objective is to construct a model for capturing the correct dynamic properties of common fluids in non-isothermal systems.

The liquid water is employed as an example fluid in the present study. We test the self-diffusivity, kinematic viscosity, surface tension, density, and thermal diffusivity at various temperature, which are all output properties instead of input parameters. And their exact values can be only obtained by simulating mDPDe system. The compressibility of liquid system can be described as $\kappa_{c}^{-1}=(\partial p/\partial\rho)_{T}/k_{B}T$ . As the pressure of the mDPD fluid system can be approximated as $p=\rho k_{B}T+\alpha A\rho^{2}+2\alpha Br_{d}^{4}\left(\rho^{3}-c\rho^{2}+d\right)$ , in which $\alpha=0.101\pm 0.001$ , $c=4.16\pm 0.02$ , and $d=18\pm 1$ 2003_Warren . The compressibility of the fluid system can be derived as $\kappa_{c}^{-1}=1+1/k_{B}T\left[2\alpha A\rho+2\alpha Br_{d}^{4}\left(3\rho^{2}-2c\right)\right]$ . On the other hand, the dimensionless coefficient corresponding to the compressibility of liquid system is $\kappa_{c}^{-1}={[L]^{3}}/{\rho k_{B}T\beta_{T}}$ 1997Groot . For liquid water at 300 K, the thermal term is $k_{B}T=4.142\times 10^{-21}\mathrm{kgm}^{2}\mathrm{s}^{-2}$ and $\beta_{T}=4.503\times 10^{-10}\mathrm{ms}^{2}/\mathrm{kg}$ . For a mDPDe system with the density $\rho=6.8$ , $A=-40$ , and $B=25$ , the corresponding scaling length is $[L]\approx 1.18~{}nm$ . The dimensionless coefficient is approximated as 130, which is larger than that from the pressure (76.3). Thus, the compressibility of the fluid system could be underestimated, as we adopt $A=-40$ and $B=25*k_{B}T$ to test the basic properties including density $\rho$ , surface tension $\sigma$ , and kinematic viscosity $\nu$ and scaling the units by comparing these parameters with those of water at different temperatures. The underestimation can be accepted in some degree as the error from the density and pressure measurement could be magnified by the coefficients in the relationship between the compressibility and the parameters of the mDPDe system. As for the diffusivity of the system, We determine the self-diffusivity of the mDPDe system by the mean-square displacement

D=\lim_{t\rightarrow\infty}\frac{1}{6t}\left\langle|\mathbf{r}(t)-\mathbf{r}(0)|^{2}\right\rangle

(12)

where $|\mathbf{r}(t)-\mathbf{r}(0)|^{2}$ is the mean-square displacement (MSD). Fig. 4 shows the MSD of the system at $T=1.0-1.25$ . The MSD has been scaled by 6.0 thus the slope of the line is the self-diffusivity for each case. It can be found in this Figure that the diffusivity increases when the temperature increases.

$T_{i}$ (K)	kBT	$\sigma$	$\nu$	$\rho$
273	0.91	9.93	7.71	7.00
283	0.94335	9.63	5.49	6.92
293	0.97667	9.34	4.04	6.81
313	1.0433	8.76	2.91	6.62
333	1.11	8.09	1.89	6.43
353	1.17667	7.33	1.41	6.25
373	1.2433	6.93	1.01	6.08

Table 1: The basic properties of mDPDe model at various temperatures.

The temperature unit was characterized as $T_{R}=300K$ . The mass, time, and length units were characterized as 2011_Arienti ; 2019_kai $M_{DPD}=L_{DPD}^{3}{d^{*}}/{d}$ , $T_{DPD}=\left(M_{DPD}{\sigma}/{\sigma^{*}}\right)^{1/2}$ , and ${L_{DPD}^{2}}/{T_{DPD}}={v^{*}}/{v}$ . The basic properties of water are marked as density $\rho$ , surface tension $\sigma$ , and viscosity $\nu$ , while the corresponding parameters in mDPDe system are those with asterisk. For every group of properties including density, viscosity, and surface tension, we can obtain corresponding length [L] $\left(T_{i}\right)$ , mass [M] $\left(T_{i}\right)$ , and time [T] $\left(T_{i}\right)$ units according to the systems of equations mentioned above. Based on the test results which are shown in Table 1, a series of these basic units can be obtained. Due to all the length units [L] $\left(T_{i}\right)$ fluctuates within one order at around $10^{-9}~{}m$ , we adopt the average value of the series results derived from the table to describe the length unit of the system at the temperature of 273 K to 373 K, as well as the time and mass units. Thus, the length, time, and mass units of the system are respectively $[\bar{L}]=1.18\times 10^{-9}~{}m$ , $[\bar{T}]=5.58\times 10^{-12}~{}s$ , and $[\bar{M}]=2.73\times 10^{-25}~{}kg$ . The kinematic viscosity of the mDPDe system is computed from the periodic Poiseuille flow method. The velocity profile obtained for the periodic Poiseuille has been shown in Fig. 1. The kinematic viscosity can be determined by fitting the velocity profile with the analytical solution $u(z)=g_{x}z(d-|z|)/2\nu$ in which $\nu$ is the kinematic viscosity. $g_{x}=0.02$ is the body-force applied on mDPDe particles and $d=10.0$ is the half-length of the computational domain in $z$ -direction. Based on the previous scaling, we can obtain the suraface tension from mDPDe simulations, as shown in Fig. 7, which is in good agreement with that in experimental data. Moreover, From the simulations we find the corresponding Schmidt number are also in good agreement with that of the real water, as shown in Fig. 8 (a). Overall, the proposed mDPDe model reproduces the correct temperature-dependent dynamic properties including the diffusivity and the kinematic viscosity as well as Schmidt number comparable to the experimental data of liquid water for the temperatures ranging from 273 K to 373 K. The relative errors of Schmidt number are generally less than 10 $\%$ of experimental data.

For the thermal conductivity in the liquid system, we adopt the $heat~{}conduction~{}analog$ of periodic Poiseuille flow method which was proposed by Li $et~{}al.$ Li2014Energy to obtain the thermal diffusivity of mDPDe fluid ranging from 273 K to 373 K. The schametic of the geometry and the temperature profile are shown in Fig. 5 and 6. And the result is represented by Prandtl number $Pr=\nu/\lambda$ , which is the ration of momentum diffusivity $\nu$ to thermal diffusivity $\lambda$ . Fig. 8 (b) shown the comparison of the temperature-dependent Prandtl number between the results of mDPDe simulations and the experimental data of liquid water. It can be observed that for different the parameter $s_{T}$ , the error or difference between the numerical results and experimental data is different. In particular, for $s_{T}=2.0$ , which is adopted as the default value, the error at high temperature is larger than 15 $\%$ . While for $s_{T}=1.0$ , the values of Prandtl number at low temperature is lower than the experimental data over than 18 $~{}\%$ . Thus, we test a series of values of $s_{T}$ at various temperatures and then obtain an empirical formula to describe accurate Prandtl number of mDPDe model, i.e. $s_{T}=-0.75*T+1.78$ . It can be observed that the Prandtl number of mDPDe fluid are consistent with the experimental Prandtl numbers. The relative errors is less than 10 $\%$ of experimental $Pr$ for the entire temperature range of 273K to 373K.

III.3 Thermal-capillary motion of a droplet

Based on the modified mDPDe, we investigate thermalcapillary motion of a droplet on a hydrophobic substrate with a temperature gradient, as shown in Fig. 9. This phenomenon has been investigated by other theoretical, experimental and numerical methods in the previous work 1989_Motions ; 2012_Xu , which suggested that the moving velocity is dependent on the the temperature difference between the ends of the droplets along the temperature gradient. The interaction between liquid-liquid particles in the conservative force is $A=-40,B=25*k_{B}T$ and the interaction between solid-liquid particles are $A_{sl}=-15$ and $B_{sl}=12.5$ , for which the static contact angle is around 120 degree.

In the simulations, we firstly set a droplet on a hydrophobic substrate with constant temperature $k_{B}T=1.0$ . After the droplet get its equilibrium state on the substate, i.e. forming a spherical cap (around 10000 timesteps), the temperature gradient is applied on the substrate and then the droplet can spontaneous move from the hot to the cold area along the substrate due to the Marangoni effect (around 90000 timesteps). The timestep is set as 0.005 and then the total time of Marangoni effect is around 450 in mDPDe units. We test three different cases with various sizes of droplets on the same substrate. Based on the scaling as mentioned above, the sizes of the water droplets in the simulations are around $8.3~{}nm$ including 8061 mDPDe particles, $9.5~{}nm$ with 14000 mDPDe particles, and $11~{}nm$ with 22408 mDPDe particles, respectively. The length of the substrate is around $70~{}nm$ along the $x$ -direction and $45~{}nm$ along the $y$ -direction and $1.2~{}nm$ along the $z$ -direction which includes 14400 mDPDe particles. The temperature of the susbtrate is from $360~{}K$ to $300~{}K$ along the $x$ -direcion. And the gradient of the temperature is around $0.857~{}K/nm$ . Bounce-back boundary condition is used in this model which has been used in the previous works Yuxiang2019Droplet . We measure the distance of droplet motion during the process of Marangoni effect and then calculate the velocity of the droplets which shows the velocity of droplets increases with larger temperature difference along the temperature gradients. It is noting that the process of Marangoni effect can be affected by the random number selection in the mDPDe simulations, which was also illustrated in the previous work (2011_Arienti, ). By comparing with the theoretical analysis and numerical results with DVDWT method 2012_Xu , as shown in Fig. 10, this modified mDPDe model can effectively capture the thermalcapillary motion of nanodroplets on hydrophobic substrates with a temperature gradient. In the comparison, the velocity $V=400m/s$ and the temperature $T=600K$ for water, which are used to scale the variables 2012_Xu .

IV Conclusions

A modified many-body dissipative particle dynamics with energy conservation (mDPDe) model for reproducing correctly the temperature-dependent properties including density, surface tension, Schmidt and Prandtl number has been proposed, which can be adopted to investigate the droplet motion on a substrate induced by Marangoni effect. The relationships between liquid-vapour surface tension of mDPDe fluid and the parameters in the conservative force are analyzed. Combining with the weighting terms of the dissipative and random forces, the temperature-dependent self-diffusivity and thermal diffusivity are obtained. The liquid-vapour surface tension, viscosity, momentum and thermal diffusivity of liquid water at various temperature ranging from 273 K to 373 K were used as a benchmark for verifying the model. The results show that the present model reproduces the correct temperature-dependent properties including liquid-vapour surface tension, Schmidt number, and Prandtl number consistent with the available experimental data of liquid water. Moreover, an corrected formula for obtaining more accurate Prandtl number at high temperature has been obtained from mDPDe simulations.

Hence, this work proposes a non-isothermal mDPD model with energy conservation (mDPDe) which is able to reproduce the correct temperature-dependent liquid-vapour surface tension, Schmidt number and Prandtl number. Furthermore, this mDPDe model is used to simulate the thermalcapillary motion of a droplet on a hydrophobic substrate with a temperature gradient. The results show that the velocity of droplets increases with larger temperature difference between the advancing and receding triple-phase contact area, which is in good agreement with the theoretical analysis and DVDWT simulations 2012_Xu . Although we test our model with liquid water, the method proposed in the present work for modelling the correct temperature-dependent properties is not limited to water only and it can be readily extended to other fluids.

V Acknowledgement

This work was supported by the National Natural Science Foundations of China (Grant No. 11872283). K. Zhang would like to acknowledge Dr. Jiayi Zhao and Mr. Damin Cao for their helpful discussions and Dr. Chensen Lin for his technological help.

References

[1] P.J. Hoogerbrugge and J.M.V.A. Koelman. Simulating microscopic hydrodynamic phenomena with dissipative particle dynamics. Europhys. Lett., 19:155–160, 1992.
[2] Shuo Chen, Nhan Phan-Thien, Boo Cheong Khoo, and Xi Jun Fan. Flow around spheres by dissipative particle dynamics. Phys. Fluids, 18(10):103605, 2006.
[3] Xuejin Li, Jiayi Guo, Yuan Liu, and Haojun Liang. Microphase separation of diblock copolymer poly(styrene-b-isoprene): A dissipative particle dynamics simulation study. J. Chem. Phys., 130(7):074908, 2009.
[4] X. J. Fan, N. Phan-Thien, S. Chen, X. H. Wu, and T. Y. Ng. Simulating flow of dna suspension using dissipative particle dynamics. Phys. Fluids, 18:063102, 2006.
[5] Zhangli Peng, Xuejin Li, Igor V. Pivkin, Ming Dao, George E. Karniadakis, and Subra Suresh. Lipid bilayer and cytoskeletal interactions in a red blood cell. Proceedings of the National Academy of Sciences, 110(33):13356–13361, 2013.
[6] Ting Ye, Nhan Phan-Thien, Boo Cheong Khoo, and Chwee Teck Lim. Dissipative particle dynamics simulations of deformation and aggregation of healthy and diseased red blood cells in a tube flow. Phys. Fluids, 26(11):111902, 2014.
[7] Zhijie Xu, Paul Meakin, Alexandre Tartakovsky, and Timothy D. Scheibe. Dissipative-particle-dynamics model of biofilm growth. Phys. Rev. E, 83(6):066702, 2011.
[8] R.D. Groot and P.B. Warren. Dissipative particle dynamics: Bridging the gap between atomistic and mesoscopic simulation. J. Chem. Phys., 107(11):4423–4435, 1997.
[9] P. B. Warren. Vapor-liquid coexistence in many-body dissipative particle dynamics. Phys. Rev. E, 68:066702, 2003.
[10] M. Arienti, W. Pan, X. Li, and G. Karniadakis. Many-body dissipative particle dynamics simulation of liquid/vapor and liquid/solid interactions. J. Chem. Phys., 134(20):204114, 2011.
[11] Z. Li, Z. W. Zhou, and G. Hu. Dissipative particle dynamics simulation of droplet oscillations in ac electrowetting. J. Adhes. Sci. Technol., 26(12-17):1883–1895, 2012.
[12] Chao Jie Mo, Li Zi Qin, and Li Jun Yang. Hypernetted-chain-like closure of ornstein-zernike equation in multibody dissipative particle dynamics. Physical Review E, 96(4):043303, 2017.
[13] Jiayi Zhao, Shuo Chen, and Nhan Phan-Thien. Rheology of polymers in many-body dissipative particle dynamics simulations: Schmidt number effect. Molecular Simulation, 44(3):1–18, 2018.
[14] D. Pan, G. Zhao, Y. Lin, and X. Shao. Mesoscopic modelling of microbubble in liquid with finite density ratio of gas to liquid. EPL, 122(2):20003, 2018.
[15] Yuxiang Wang, Jiayi Zhao, Dingni Zhang, Meipeng Jian, and Xiwang Zhang. Droplet sliding: The numerical observation of multiple contact angle hysteresis. Langmuir, 35(30):9970–9978, 2019.
[16] K. Zhang, Z. Li, M. Maxey, S. Chen, and G. E. Karniadakis. Self-cleaning of hydrophobic rough surfaces by coalescence-induced wetting transition. Langmuir, 35(6):2431–2442, 2019.
[17] P Español. Dissipative particle dynamics with energy conservation. EPL, 40(6):631–636, 1997.
[18] R. Qiao and P. He. Simulation of heat conduction in nanocomposite using energy-conserving dissipative particle dynamics. Molecular Simulation, 33(8):677–683, 2007.
[19] Abu-Nada and Eiyad. Natural convection heat transfer simulation using energy conservative dissipative particle dynamics. Phys. Rev. E, 81(5):056704, 2010.
[20] E. O. Johansson, T. Yamada, B. Sunden, and J. Yuan. Modeling mesoscopic solidification using dissipative particle dynamics. International Journal of Thermal Sciences, 101:207–216, 2016.
[21] M. Ripoll, P. Espaol, and M. H. Ernst. Dissipative particle dynamics with energy conservation: Heat conduction. International Journal of Modern Physics C, 9:1329–1338, 1998.
[22] P. He and R. Qiao. Self-consistent fluctuating hydrodynamics simulations of thermal transport in nanoparticle suspensions. Journal of Applied Physics, 103(9):783–788, 2008.
[23] Zhen Li, Yu Hang Tang, Huan Lei, Bruce Caswell, and George Em Karniadakis. Energy-conserving dissipative particle dynamics with temperature-dependent properties. Journal of Computational Physics, 265:113–127, 2014.
[24] Yamada, Toru, Johansson, Erik, O., Yuan, Jinliang, Sunden, and Bengt. Dissipative particle dynamics simulations of water droplet flows in a submicron parallel-plate channel for different temperature and surface-wetting conditions. Numerical Heat Transfer Part A, 70:1–18, 2016.
[25] Liwei Wang, Jinzhao Dai, Pengfei Hao, Feng He, and Xiwen Zhang. Mesoscopic dynamical model of ice crystal nucleation leading to droplet freezing. ACS Omega, 5(7):3322–3332, 2020.
[26] J. B. Avalos and A. D. Mackie. Dissipative particle dynamics with energy conservation. EPL, 40:141–146, 1997.
[27] F. Brochard. Motions of droplets on solid surfaces induced by chemical or thermal gradients. Langmuir, 5(2):432–438, 1989.
[28] X. Xu and T. Qian. Thermal singularity and droplet motion in one-component fluids on solid substrates with thermal gradients. Phys. Rev. E, 85(6 Pt 1):061603, 2012.