New modified weight function for the dissipative force in the DPD method to increase the Schmidt number

Dissipative Particle Dynamics (DPD) is a stochastic mesoscopic simulation approach [1,2], bridging the gap between atomistic and continuum fluid descriptions. The DPD particles represent coarse-grained molecules which move together in a Lagrangian fashion. Hoogerbrugge and Koelman [1] introduced this scheme for the first time and its basis in statistical mechanics becomes well known by means of the work by Español, Warren [3] and Marsh [4]. Due to the soft-particle interactions in DPD, time steps and length scales are much larger in comparison to the Molecular-Dynamics (MD) method. Besides that, because of the large number of particles in MD, this method is too inefficient for mesoscale simulations. Although there are simulation methods in this scale-like Lattice Gas Automata (LGA) [5] and Lattice Boltzmann Method (LBM) [6], it is difficult to use LGA or LBM to simulate complex fluids. Specifically, DPD appears as a successful approach in the simulation of complex fluids, such as suspensions of polymers, DNA, colloids and cells in blood [7–9].

Calculating the transport coefficients of a DPD fluid has become a challenge in this technique, which simulates fluids with Schmidt number about unity [2]. The Schmidt number, which is the ratio of the kinematic viscosity to the diffusivity of a solvent particle, is of order one for gases, and is not suitable for typical liquids such as water, with the Schmidt number on the order of one thousand. This represents the very small ratio of momentum transfer rate to that of mass transfer for liquids. Some attempts have been made to obtain the Schmidt number with the correct order of magnitude. For instance, Fan et al. [10] have modified the DPD technique to improve the Schmidt number and the dynamic characteristics of the DPD fluids. They consider the exponent of the weight function in the random force equation to be less than 1. It yields to a Schmidt number of about 10. Moreover, they increase the cut-off radius. So they have simulated polymers and DNA molecules with a higher value of the Schmidt number. Another model has been proposed by Lowe [11]. In this approach, Lowe uses a new thermostat and the same conservative forces, but with this method there is a limitation for wall-bounded systems as it is not clear for them how to impose the no-slip boundary condition accurately. This limits the use of the method for periodic domains.

In this paper, to increase the viscosity and the Schmidt number of DPD fluids, we have proposed a new formulation of the weight function for the dissipative and random forces. Then we compare the effect of the new modification on the viscosity and the Schmidt number with that of the standard one and the version proposed recently by Fan et al. [10]. By this comparison, we will show that our new modification increases the dissipative part of the Schmidt number to 200 times larger than the standard one. Also, we increase the viscosity and the Schmidt number by reducing the energy of the particles. To do this, we reduce the temperature of the system. To investigate our new formula, we imply the energy reduction to the modified weight function and show that by using the new function, the Schmidt number can be increased up to 200. The innovation of this work which motivated us to introduce it, is to promote the previous modification by Fan et al. [10] and increasing Sc which was an improvement for the method and the function is originally ours.

The rest of this paper is organized as follows: we first present an outline of the standard numerical scheme, Dissipative Particle Dynamics. Subsequently we propose a modified DPD formulation. Then we present the simulation results using our new modification with two levels of energy of the system. Finally the conclusion will be given.

Coarse-grained particles in the DPD method interact via pairwise-additive forces. These specified forces conserve momentum and are composed of three terms: repulsion conservative force $F_{ij}^{C}$, dissipation force $F_{ij}^{D}$ and random force $F_{ij}^{R}$, respectively. All of these forces vanish above a certain cut-off radius r_c, which defines the length scale. Since Newton's laws govern the motion of the DPD particles, for the i-th particle the following relation will be satisfied:

$$\begin{eqnarray} \begin{array}{rll} v_{i} &=\displaystyle \frac{\mathrm{d}r_{i}}{\mathrm{d}t},\\ f_{i} &=\sum_{j\ne i} \left( +F_{ij}^{D} +F_{ij}^R\right) =m_{i} \displaystyle \frac{\mathrm{d}v_{i}}{\mathrm{d}t}, \end{array} \end{eqnarray} \tag{ 1 }$$

where m_i, r_i and v_i are the mass, position and velocity of the i-th particle, respectively, and the sum runs over all other particles located in a sphere with a cut-off radius r_c. The soft conservation force traditionally is given by [2]

$$\begin{equation} F_{ij}^{C} =\begin{cases}a_{ij} \left(1-\displaystyle\frac{r_{ij} }{r_{c} }\right){\hat{r}}_{ij}& (r_{ij} <r_{c}),\\ 0& (r_{ij} \geq r_{c}),\end{cases} \end{equation} \tag{ 2 }$$

where $a_{ij} =\sqrt {a_{i} a_{j} }$ and a_i, a_j represent the maximum repulsion strength parameters for particles i and j, respectively. $\vec{r}_{ij} =\vec{r}_{i} -\vec{r}_{j}$ is the distance vector between particles i and j (pointed from particle j to i), $r_{ij} =\left| {\vec{r}_{ij}}\right|$ and $\hat{r}_{ij} =\vec{r}_{ij} /r_{ij}$. The dissipative and random forces have the following form [2]:

$$\begin{eqnarray} F_{ij}^{D} &=-\gamma \omega^{D}(r_{ij} )(\hat{r}_{ij}.v_{ij} )\hat{r}_{ij},\nonumber \\ F_{ij}^{R} &= \sigma \omega^{R}(r_{ij} )\xi_{ij} \Delta t^{-1/2}\hat{r}_{ij}, \end{eqnarray} \tag{ 3 }$$

where γ and σ are two coefficients, which characterize the amplitude of the dissipative and random forces, respectively. $\omega^{D}$ and $\omega^{R}$ are r-dependent weight functions which are zero for r > r_c. $v_{ij} =v_{i}-v_{j}$ is the relative velocity, $\Delta t$ is the interval of time and $\xi_{ij}$ is a symmetric Gaussian, normally distributed random variable, with zero mean and unit variance $(\xi_{ij} =\xi_{ji})$.

A fluctuation dissipation relation relates the random and dissipative forces, which implies,

$$\begin{equation} \omega^{D}(r)=[\omega^{R}(r)]^{2},\quad \sigma^{2}=2\gamma k_{B} T, \end{equation} \tag{ 4 }$$

where k_BT is the target Boltzmann temperature of the system. If the fluctuation-dissipation relation is satisfied, this method produces a correct $(N, V, T)$ ensemble [2,3]. The weight function is usually chosen as follows:

$$\begin{equation} [\omega^{R}(r)]= \begin{cases} \left(1-\displaystyle\frac{r}{r_{c}}\right) & (r<r_{c}),\\ 0& (r\ge r_{c}). \end{cases} \end{equation} \tag{ 5 }$$

In this paper, all particles have the same mass m. The particle mass, temperature and interaction range have been chosen as units of mass, energy, and length, respectively.

Here we use the velocity-Verlet algorithm of Allen and Tildesley [12] and also Groot and Warren [2] to integrate the equations of motion.

The Schmidt number is an important parameter to characterize the dynamic behavior of fluids especially complex fluids. Usually the standard version of DPD will be exerted, which simulates fluids with Schmidt number of order one [1]. The lack of non-central dissipative forces in addition to the soft conservative forces in DPD, yields to diminishing the viscosity and Schmidt number of the simulating fluid [10]. In reality, a single DPD particle in a flow is subjected to both translational and rotational movements while only central forces are considered in standard DPD simulations and non-central shear forces between dissipative particles are ignored. Español et al. [13–15] improved the method by developing a Fluid Particle Model (FPM) [14]. Compared to the standard DPD, additional non-central shear components of the dissipative and random forces are included in FPM. Therefore, we cope with both linear- and angular-momentum conservation in this method, which results in a higher computation cost, compared to DPD.

To improve the dynamic characteristics of a simulating system, recently Fan et al. [10] exerted a modification in the DPD method. They proposed another weight function for the dissipative and random forces which is given by

$$\begin{equation} [\omega^{R}(r)]= \begin{cases} \left(1-\displaystyle\frac{r}{r_{c}}\right)^{s}& (r < r_{c}),\\ 0& (r\ge r_{c}). \end{cases} \end{equation} \tag{ 6 }$$

In their paper they have shown that choosing $s \le 1$, yields a Schmidt number up to 10. Increasing the cut-off radius is one way to increase the Schmidt number. As a result of this new cut-off the computational expenses also increases, which is a considerable cost. Note that the computational cost increases proportionally to the cube of the larger cut-off radius.

In the present work, we propose the following weight function for dissipative and random forces in DPD to achieve a Schmidt number of larger values:

$$\begin{equation} [\omega^{R}(r)]= \begin{cases} 1-\left(\displaystyle\frac{r}{r_{c}}\right)^{k} & (r < r_{c}),\\ 0& (r\ge r_{c}). \end{cases} \end{equation} \tag{ 7 }$$

If k = 1, the conventional weight function in standard DPD systems, is achieved. Considering k > 1 yields to increasing the Schmidt number even more than the modification proposed by Fan et al. According to eq. (4) (fluctuation dissipation relation) the weight function of dissipative and random forces is related to producing an isothermal system. The idea behind this new function is to enhance the strength of the dissipative interactions of a particle with neighbors in a cut-off distance which yields to increasing the viscosity and hence the Schmidt number of the fluid. Since $r/r_{c}$ in the new proposed formulation is less than 1, increasing k leads to strengthening the dissipative forces. Figure 1 illustrates the effect of this new function vs. $r/r_{c}$ for various k.

Zoom In Zoom Out Reset image size

As shown in fig. 1, it is obvious that with larger values of k, the weighting function and as well as the dissipative interactions are stronger in a distance, then more viscosity is recovered. Figure 2 illustrates the comparison of this new modification with that of Fan et al. [10] and also the standard form which shows a higher level of the function in the range of $r/r_{c}$ for the new formulation relative to the standard one and the one proposed by Fan et al. [10].

Zoom In Zoom Out Reset image size

Using eq. (7) and similarly to Fan et al., we can produce an approximate deliberation for the dynamic properties of the new function. With a simplified assumption for the radial distribution function (a uniform density) [2], the viscosity of a simple DPD system has two parts: dissipative contribution due to dissipative forces between particles in different streamlines and kinetic contribution because of the diffusing particles across streamlines. The dissipative part of the viscosity is derived [2] as

$$\begin{equation} \eta^{D} = \frac{2\pi \gamma \rho^{2}}{15}\int\limits_0^\infty {r^{4}} \omega^{D}(r)\mathrm{d}r. \end{equation} \tag{ 8 }$$

With the new weighting function, we can represent $\eta^{D}$ by

$$\begin{equation} \eta^{D}=\frac{2\pi \gamma \rho^{2}r_{c}^{5} }{15}\left(\frac{1}{5}-\frac{2}{5+k}+\frac{1}{5+2k}\right). \end{equation} \tag{ 9 }$$

Increasing k will result in a stronger particle interaction and more dissipative viscosity than the conventional one in the standard DPD method. The self-diffusion coefficient of a DPD particle is obtained [2] as

$$\begin{equation} D=\tau k_{B} T, \end{equation} \tag{ 10 }$$

where

$$\begin{equation} \frac{1}{\tau }=\frac{4\pi \gamma \rho }{3}\int\limits_0^\infty r^{2}\,\omega^{D}(r)\,\mathrm{d}r \end{equation} \tag{ 11 }$$

and the kinetic contribution of the viscosity is proportional to D with the proportion coefficient $\rho/2$:

$$\begin{equation} \eta^{k}=\frac{\rho D}{2}=\frac{3k_{B} T}{8\pi \gamma \int\limits_0^\infty {r^{2}\,\omega^{D}(r)\,\mathrm{d}r}}. \end{equation} \tag{ 12 }$$

The kinematic viscosity of a DPD fluid is [2]

$$\begin{equation} \eta =\eta^{D}+\eta^{k}. \end{equation} \tag{ 13 }$$

Equations (8), (10) and (12), which specify the dynamic properties of a fluid, are changeable with different weight functions. The dynamic properties of the new modified system proposed in this paper with that of Fan et al. [10] and the standard DPD system, are presented in table 1.

Table 1:. Dynamic properties of DPD fluid with different weight function.

Properties	Conventional weight function (eq. (5))	Previous modified weight function, $s \le 1$ (eq. (6))	New modified weight function $k \ge 1$ (eq. (7))
Diffusivity (eq. (10))	$\displaystyle\frac{45k_{B} T}{2\pi \gamma \rho r_{c}^{3} }$	$\displaystyle\frac{3k_{B} T}{4\pi \gamma \rho r_{c}^{3} \left(\displaystyle\frac{1}{2s+1}-\frac{1}{s+1}+\frac{1}{2s+3}\right)}$	$\displaystyle\frac{3k_{B} T}{4\pi \gamma \rho r_{c}^{3} \left(\displaystyle\frac{1}{3}-\frac{2}{3+k}+\frac{1}{3+2k}\right)}$
Viscosity (eq. (13))	$\displaystyle\frac{\rho D}{2}+\frac{2\pi \gamma \rho^{2}r_{c}^{5} }{1575}$	$\begin{array}{l} \displaystyle\frac{\rho D}{2}+\frac{2\pi \gamma \rho^{2}r_{c}^{5} }{15}\biggl(\frac{1}{2s+1}-\frac{2}{s+1} +\,\displaystyle\frac{6}{2s+3}-\frac{2}{s+2}+\frac{1}{2s+5}\biggr) \\ \end{array}$	$\displaystyle\frac{\rho D}{2}\!+\!\frac{2\pi \gamma \rho^{2}r_{c}^{5} }{15}\left(\frac{1}{5}\!-\!\frac{2}{5+k}\!+\!\frac{1}{5+2k}\right)$
Schmidt number $(\mathrm{Sc}=\eta / \rho D)$	$\displaystyle\frac{1}{2}+\frac{(2\pi \gamma \rho r_{c}^{4} )^{2}}{70875k_{B} T}$	$\begin{array}{l} \displaystyle\frac{1}{2}+\frac{(4\pi \gamma \rho r_{c}^{4})^{2}}{90k_{B} T}\biggl(\frac{1}{2s+1}-\frac{2}{s+1}+\frac{6}{2s+3} -\,\displaystyle\frac{2}{s+2}+\frac{1}{2s+5}\biggr)\left(\frac{1}{2s+1}-\frac{1}{s+1}+\frac{1}{2s+3}\right) \\ \end{array}$	$\begin{array}{l} \displaystyle\frac{1}{2}+\frac{(4\pi \gamma \rho r_{c}^{4} )^{2}}{90k_{B} T}\biggl(\frac{1}{3}-\frac{2}{3+k}+ \displaystyle\frac{1}{3+2k}\biggr)\left(\displaystyle\frac{1}{5}-\frac{2}{5+k}+\frac{1}{5+2k}\right) \end{array}$

Figure 3 illustrates the ratio of the dissipative component of the Schmidt number for the new modified system denoted by $\mathrm{Sc}_{\mathrm{N}}^{\mathrm{D}}$ relative to that of the standard DPD system denoted by $\mathrm{Sc}_{\mathrm{S}}^{\mathrm{D}}$, vs. k given in eq. (7).

Zoom In Zoom Out Reset image size

Figure 3 shows the effect of the new weight function on the standard one on increasing the dissipative part of the Schmidt number in the DPD method. At $k =1$, the conventional weight function in a standard DPD system is obtained. It can be seen that the dissipative part of the Schmidt number for the new modification is approximately 200 times larger than that when using a conventional weight function.

Figure 4 illustrates the dissipative component of the Schmidt number for the new modified system relative to that of Fan et al. [10] with $s = 0.25$ (eq. (6)) which is denoted by $\mathrm{Sc}_{\mathrm{P}}^{\mathrm{D}}$.

It can be seen that k = 5 (eq. (7)), will result in the same value as using $s = 0.25$ (eq. (6)). For higher values of k, $Sc_{\mathrm{N}}^{\mathrm{D}}/Sc^{\mathrm{D}}_{\mathrm{P}}$ increases and becomes higher than one. It is 5.5 at $k= 120$.

Zoom In Zoom Out Reset image size

Several simulations have been done to evaluate the modified versions of the weight function in a Couette flow simulation. Using Lees-Edwards periodic boundary conditions [16], the viscosity of the fluid flow is measured. The Schmidt number, dynamic viscosity and self-diffusion as a function of the exponent for the weight functions are represented. According to table 1, Sc is proportional to $1/k_{B}T$. This means that lower energy level of the particles leads to an increase of the Schmidt number. In fact, reducing the energy of the particles leads to a decrease of the diffusion coefficient and also to an increase of the viscosity effects and thus increasing the Schmidt number. For different values of k, we calculate the Stokes-Einstein radius [17]. This parameter depends only on the repulsive coefficient in conservative force and determines the sphere, which is approximately impenetrable with other particles [18]. Then, since it is independent of the dissipative weight function, it must have an approximate constant value with different weight functions and exponents.

By means of the Stokes-Einstein equation and measurement of the transport coefficients, self-diffusion and dynamic viscosity, Stokes-Einstein radius R_SE can be calculated as [17]

$$\begin{equation} R_{SE} =\frac{k_{B} T}{CD\pi \mu}, \end{equation} \tag{ 14 }$$

where D is the self-diffusion coefficient of a particle which has Brownian motion in a domain, μ is the dynamic viscosity of the fluid. C is the resistance coefficient which depends on the hydrodynamic boundary conditions on the surface of a particle, which is assumed to be sphere. For a solid sphere, whose surfaces the solvent adheres to, Einstein obtained $C = 6$ and we use it in our calculations. To determine the self-diffusion coefficient, the Mean-Square Displacement (MSD) relation of Einstein in the long time is used which is given by [19]

$$\begin{equation} CD_{\infty } =\lim\limits_{t\to \infty} \frac{\left\langle {\left| {r(t)-r(0)} \right|^{2}} \right\rangle}{t}, \end{equation} \tag{ 15 }$$

where $\langle\cdot \rangle$ shows the ensemble average, and r(0) and r(t) are the particle position at origin time and time t, respectively.

The parameters used in our DPD simulations, are summarized in table 2. In this table, ρ is the number density of the particles, a is a repulsive parameter in the conservative force, γ and σ are amplitudes of the dissipative and random forces, respectively, T is the equilibrium temperature and r_C is the cut-off radius.

The properties μ, D, $\mathrm{Sc} =\mu/(\rho D)$ and Stokes-Einstein radius for different DPD fluids (various values of the exponent k) are obtained and plotted in figs. 5 and 6, respectively.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Figure 5 illustrates that increasing the exponent k will increase the Schmidt number. Similarly to fig. 3, the effect of increasing k is more pronounced at lower values. In fig. 6 the values obtained for the particle radius R_SE with parameters listed in table 1 remain considerably constant with different exponent k. Similarly to the results obtained by Pan [20], the repulsive parameters in the DPD conservative force are held constant. Based on this result, we can conclude that the Stokes-Einstein radius depends only on the DPD conservative parameter and hence the radial distribution function (RDF) is alike for fluids with different weight functions.

Figures 7 and 8 show the computed values of μ, D, Sc and Stokes-Einstein radius using the weight function of Fan et al. [10], as a function of s, respectively.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Lowering the energy level of the particles, which decreases the diffusion coefficient and increases the viscosity, will also increase the Schmidt number. Thus, new simulations have been done at a lower energy level, by using $k_{B}T= 0.2$ and also the new weight function in a Couette flow. The parameters μ, D, and Sc with lower energy and new weight function are plotted as a function of k in fig. 9. This figure illustrates that a Schmidt number of nearly 200 can be obtained.

Zoom In Zoom Out Reset image size

Figure 10 shows the Stokes-Einstein radius as a function of k with the new weight function and at lower energy level.

Zoom In Zoom Out Reset image size

In this figure, the values obtained for the particle radius R_SE with lower energy remain considerably constant as well with different exponent k.

Calculating the transport coefficients of a DPD fluid has become a challenge in this technique, which simulates fluid flows with the Schmidt number of order one. This value of the Schmidt number is not suitable for real liquids, which have the Schmidt number of order of one thousand. In order to simulate liquid flows with higher values of the Schmidt number, one needs to use an improved DPD method. In this paper, we proposed a new weighting function to achieve a solution for simulating fluid flows with higher values of Schmidt numbers.

To prove the revenue of our new function, we plotted all three versions of the weight function vs. the distance between particles. We have shown the power of the new weight function compared to other versions to achieve higher values of the Schmidt number. Strengthening the weight function leads to the increase of the dissipative force and hence the increase of the viscosity.

The dissipative component of the Schmidt number which is up to 200 times larger than the standard one has been obtained, using the new modification. Results are compared with that in the improved DPD method of Fan et al. [10] which is 35.5 times larger than the standard one. Reducing the energy of the particles is another approach that we exert to increase the Schmidt number of the simulating fluid. Eventually the new modified weight function is exerted in the Couette flow with two levels of energy. The diffusion and viscosity coefficients and the Schmidt number are calculated and compared with the standard and the improved DPD method of Fan et al. [10]. We will see, by using the new function and $k_{B}T=0.2$, that the Schmidt number values can reach up to 200.