- Email: [email protected]
Is the cutoff radius in DPD simulations with a fluid of constant density arbitrary?
Computer Physics Communications 183 (2012) 1714–1718
Contents lists available at SciVerse ScienceDirect
Computer Physics Communications www.elsevier.com/locate/cpc
Is the cutoff radius in DPD simulations with a fluid of constant density arbitrary? Isaac Gabriel Salib ∗ , Shimon Haber Technion – Israel Institute of Technology, Haifa 32000, Israel
a r t i c l e
i n f o
Article history: Received 1 September 2009 Received in revised form 22 February 2010 Accepted 15 March 2012 Available online 19 March 2012 Keywords: Coarse graining Dissipative Particle Dynamics Integration schemes Soft condensed matter
a b s t r a c t The DPD (Dissipative Particle Dynamics) method was used to simulate the flow of a fluid occupying the space inside a rotating infinite long cylinder. The fluid was initially at rest and rigid body rotation was expected to be reached. However, as the simulation proceeded, the DP’s (Dissipative Particles) began accumulating near the wall, leaving the cylinder’s core empty of DP’s, obviously, an aphysical solution. To solve this problem it was suggested to set the cutoff radius as a function of the DP’s number density. A conceptual explanation is suggested. © 2012 Elsevier B.V. All rights reserved.
1. Introduction There is a growing interest in simulating fluid flows in biological systems. Simulating blood flow in veins and arteries, or air laden particles flow in lungs is of great importance in medical and biotechnological researches. However, it is almost impossible to simulate these flows using classical simulation techniques: the fluids are complex, may include many species (such as blood flow consisting of red and white blood cells platelets, etc.) and bounded by intricate boundaries (such as soft tissues). Also nanoand micro-flows need a simulation technique which can capture their hydrodynamic behavior, taking in consideration their nanoand micro-scale complexity. A “coarse-grained” method for simulating fluid flows was suggested by Hoogerbrugge and Koelman [1,2]. This method, known as Dissipative Particle Dynamics (or DPD), is used to model a fluid by a finite number of particles, each of which represents a large number of molecules. The hydrodynamic behavior is predicted through tracing the particles, basically a Lagrangian approach. The force system exerted on the DP’s includes soft repulsive conservative, dissipative and random forces. The last two forces were termed by Hoogerbrugge and Koelman as the system’s “thermostat”. The dissipative term tends to “cool down” the system while the random term tends to “heat up” the system. The presence of both terms enables local momentum conservation and the prediction of hydrodynamic behavior in the macroscopic scale. The DPD model was modified by Español et al. [3,4]. They derived the Fokker–Planck
equation and showed how the DPD Langevin equations are related to statistical mechanics. Furthermore, in DPD, several no-slip boundary methods were employed [5]. Haber et al. [6] suggested a benchmark test of fluid flow between two rotating concentric cylinders. In their article, comparison was drawn between these methods and the bounceback reflection method (reflecting back the particles which pass the rigid boundary with an opposite velocity) was shown to yield the most accurate results. The two rotating concentric cylinders test was found to be a useful benchmark since this test case has an analytical solution and it represents a closed system where no artificial periodic boundary conditions are needed. Haber et al. [6] noted that the DP number density at the boundary was somewhat higher than between the cylinders. Surprisingly, longer runs revealed that after enough time steps, the DP’s tend to accumulate at the boundaries. The phenomenon was augmented at high cylinder rotating speeds for which fewer time steps were required. In order to investigate and resolve this issue, we used a simpler 2D model: a fluid, originally at rest, inside a rotating cylinder. The results are explained and a novel perspective of the DPD method is brought up, that is of particular interest for flows bounded by curved walls. 2. Method The basic and most common postulated DPD scheme introduced by Español et al. [3], takes the following form:
*
Corresponding author. E-mail addresses: [email protected] (I.G. Salib), [email protected] (S. Haber). 0010-4655/$ – see front matter © 2012 Elsevier B.V. All rights reserved. http://dx.doi.org/10.1016/j.cpc.2012.03.012
dri =
pi mi
dt
(1)
I.G. Salib, S. Haber / Computer Physics Communications 183 (2012) 1714–1718
1715
Fig. 1. DP’s accumulation near the boundary (each asterisk stands for a DP): Out of the 900 DP’s only a few remain inside the domain while the rest are at the boundary.
dpi =
fiCj (r i j ) +
j =i
+
The weight functions frequently used posses the following form:
−γ ω D (r i j )(vi j · rˆ i j )ˆri j dt
j =i
σ ω R (ri j ) dW i j rˆ i j
(2)
j =i
where ri j and vi j are the position and velocity of particle i relative to particle j, rˆ i j is the unit vector in the direction of ri j , pi and mi are the momentum and the mass of particle i respectively, and dt is the time increment. The conservative force is fiCj and the constants γ and σ are the coefficients of the dissipative and random forces respectively. The weight functions, ω D (r i j ) and ω R (ri j ), correspond to the soft matter interactions (all forces act up to a maximum distance range called cutoff radius rc ) and dW i j is an independent increment of the Weiner process for which Itô calculus rule implies:
dW ii dW j j = [δi j δi j + δi j δi j ] dt
(3)
From the corresponding Fokker–Planck equation, one gets Eq. (4):
ω2R (r ) = ω D (r ),
σ 2 = 2k B T γ
(4)
where k B is the Boltzmann’s constant and T is the equilibrium temperature. The conservative force commonly used is of the following form:
fiCj =
ai j (1 − r i j /rc )ˆri j
at r i j rc
0
at r i j > rc
(5)
where ai j is the maximum repulsion force constant. This latter parameter is taken to be uniform for all interparticle interactions and is simply replaced by a. It was established by Groot and Warren [7] through the fluid compressibility.
ω2R (r ) = ω D (r ) =
(1 − r i j /rc )2 at r i j < rc 0 at r i j > rc
(6)
The simulations were conducted using the bounce-back boundary condition. As the simulation proceeded, the DP’s molecules started to accumulate at the boundary (Fig. 1). This phenomenon occurred in case either the Euler or the Velocity Verlet [8] algorithms were used. It was also observed that the size of the time increment used had no effect on the DP accumulation at the boundary (except for the total number of time steps taken for it to occur). The DP’s which passed the boundary were brought back artificially by the bounce-back boundary conditions. These DP’s are usually reintroduced into the domain at positions close to the boundary (this may explain the fact that higher densities are obtained near the boundaries). But, as the simulation proceeded more DP’s were driven towards the boundaries. There are two mechanisms which, theoretically, hinder the DP’s boundary accumulation: 1. The fluctuating force which provides the DP with the random velocity and displacement which, by forcing the DP’s to intermingle randomly, could maintain a constant density over the whole domain. However, for simulations in which the rotational cylinder speed is much higher than the kinetic one, i.e. k B T /m, the random force is too small to achieve this balance. 2. The repulsive conservative force could also prevent the accumulation phenomenon. As noted earlier, since the DP number density at the boundary is higher than that in the rest of the domain, we would expect that the repulsive force exerted by
1716
I.G. Salib, S. Haber / Computer Physics Communications 183 (2012) 1714–1718
Fig. 2. Schematic sketch of two adjacent DP’s with interparticle distance that is (A) smaller than the sum of their radii and (B) larger than the sum of their radii.
the DP’s near the boundary would repel the particles back into the domain. This force is controlled by two parameters: cutoff radius and the force amplitude or a in Eq. (5). We tried first to adjust only the force amplitude. We found out that by increasing the conservative force amplitude, the accumulation phenomenon is just postponed, i.e. for the same magnitude of time increment, more time steps were needed to reach accumulation of DP’s near the wall. Apparently, only the adjustment of the cutoff radius and the magnitude of the conservative force could control the DP number density. 3. Theory The cutoff radius (coined rc ) designates the radius of influence, beyond which the DP cease to interact with other DP’s. This parameter is thought to be chosen arbitrarily. Different conclusions were drawn for the optimal cutoff radius, of which, the most acceptable, is that the cutoff radius has to define a spherical volume around the DP that contains at least three or more DP’s on average. It is worthy to note that setting the cutoff radius is of practical importance. For a given time increment, the smaller the cutoff radius, more time steps are needed to reach equilibrium. To our understanding, the cutoff radius also possesses a physical significance and not just a numerical one. The DP’s are, as noted earlier, lumps of molecules. Adjacent DP’s meddle. In the meddling zone, molecules bombard and momenta transfer. The weight function denotes that closer the DP’s, the larger is this zone and the greater is the resultant interactive forces (as illustrated in Fig. 2). Thus, the cutoff radius between two adjacent DP’s should equal the sum of their radii.
We shall limit our simulations here to constant density cases. The DP’s radii are uniform. Thus we expect that the cutoff radius is constant over the whole solution domain. 4. Results We simulated a fluid in a cylinder of radius R = 0.001 m that rotates with a speed of ω = 1 rad/s. Each DP is of m = 3.4907 × 10−12 kg mass (for water that means 1.16 × 1014 molecules in each DP). The simulation is conducted for T = 273 K temperature and 900 DP’s. Although, a wide range of numerical values for the dissipation coefficient were used, we choose here to present the results achieved for
γ m
=
1575 2π n2 rc5
μ
(7)
where μ is the shear viscosity and n is the initial uniform particle density. The conservative force amplitude used is taken as stated in Groot and Warren [7]. For the time increment, we followed the footsteps of Haber et al. [6]. In their choice of the time increment they took into account all the time scales in the test case rc / v T , m/γ and v T m/a, including, rc /ω R, that is linked to the boundary velocity, ω R, where R is the cylinder’s radius in terms of rc . The dimensionless time increment chosen is as follows:
dt = 10− f min 1, v T m/γ rc , v 2T m/arc , v T /ω R
where f is the numerical accuracy sought. For f we tried the following values: 2, 3, 4, and 5. Since the model used here is a 2D one, the area which each DP possesses is simply the cross section area of the cylinder divided
I.G. Salib, S. Haber / Computer Physics Communications 183 (2012) 1714–1718
1717
Fig. 3. DP’s inside the rotating cylinder after 100,000 time steps for rc = 0.0667 mm (each asterisk stands for a DP).
Fig. 4. DP’s inside the rotating cylinder after 100,000 time steps for rc = 0.4 mm (each asterisk stands for a DP).
by the number of the DP’s. Assuming that the DP is circular, the cutoff radius would be equal to twice this radius (the factor of 2 is the optimum for achieving the desired results):
R rc = 2 × √ = 6.67 × 10−5 m N
(8)
Simulation was conducted for 100,000 time steps. Fig. 3 illustrates the image of the DP’s at the end of the simulation. As we can see the DP’s are uniformly distributed and the hydrodynamic behavior is that of rigid body rotation as expected. In comparison, we have conducted a simulation with cutoff radius five times larger than that used in Eq. (7). The illustra-
1718
I.G. Salib, S. Haber / Computer Physics Communications 183 (2012) 1714–1718
tion of the DP’s locations after 100,000 time steps is presented in Fig. 4. As we can see the DP’s that were located near the boundary started accumulating at the boundary. As a result, unlike the uniform number density illustrated in Fig. 3, the DP number density is higher at the boundary. As the simulation proceeds, more and more DP’s accumulate at the wall and finally reach the distribution shown in Fig. 1. 5. Conclusions Most previous benchmark problems used to test the DPD method employed planar boundaries for Couette or Poiseuille flows where the flow does not experience centrifugal forces. Haber et al. [6] benchmark problem was the first to include a curved surface for which the continuous solution includes centrifugal effects. The DPD method is based on modeling lumps of fluid by DP particles. However, in the DPD model, the volume that each DP occupies is not set, only the DP mass is a priori assumed. Therefore, the density of the fluid may not be maintained. On using the rotating cylinder 2D simulation, this flaw of not keeping constant density is noted clearly; the DP number density could be much higher at the boundary if the cutoff radius is not properly picked. In order to avoid such a numerical artifact, we suggest that the cutoff radius is a measure of the DP’s extent. Adjusting the cutoff radius and the conservative force amplitude enable us to reach an almost constant density and prevent DP gathering near the boundary. Future research would attempt to relate the cutoff radius magnitude to the thermodynamical parameters in order to allow compressible flow simulations.
In this article, there is no account for using a different conservative force. Nevertheless, we also tried to seek a solution, by using a Lennard-Jones potential between the particles. However, when the interparticle distance tends to zero, this potential tends to infinity and so does the DP displacement, i.e. a very small time increment has to be used. In this case the DPD scheme is impractical. We also tried to replace the conservative force by a repulsive attractive spring force. In this case aggregates were formed throughout the domain, obviously an aphysical result. References [1] P.J. Hoogerbrugge, J.M.V.A. Koelman, Simulating microscopic hydrodynamic phenomena with dissipative particle dynamics, Europhys. Lett. 19 (1992) 155. [2] P.J. Hoogerbrugge, J.M.V.A. Koelman, Dynamic simulations of hard-sphere suspensions under steady shear, Europhys. Lett. 21 (1993) 363. [3] P. Español, P. Warren, Statistical mechanics of dissipative particle dynamics, Europhys. Lett. 30 (1995) 191. [4] P. Español, Hydrodynamics from dissipative particle dynamics, Phys. Rev. E 52 (1995) 1734. [5] D.C. Visser, H.C.J. Hoefsloot, P.D. Iedema, Comprehensive boundary method for solid walls in dissipative particle dynamics, J. Comput. Phys. 205 (2005) 626; I.V. Pivkin, G.E. Karniadakis, A new method to impose no-slip boundary conditions in dissipative particle dynamics, J. Comput. Phys. 207 (2005) 114; M. Revenga, I. Zúñiga, P. Español, Boundary model in DPD, Int. J. Mod. Phys. C 60 (1998) 1319. [6] S. Haber, N. Filipovic, M. Kojic, A. Tsuda, Dissipative particle dynamics simulation of flow generated by two rotating concentric cylinders. Part I: Boundary conditions, Phys. Rev. E 74 (2006) 046701. [7] R.D. Groot, P.B. Warren, Dissipative particle dynamics: Bridging the gap between atomistic and mesoscopic simulation, J. Chem. Phys. 107 (1997) 4423. [8] I. Pagonabarragga, M.H.J. Hagen, D. Frenkel, Self-consistent dissipative particle dynamics algorithm, Europhys. Lett. 42 (1998) 377.
Contents lists available at SciVerse ScienceDirect
Computer Physics Communications www.elsevier.com/locate/cpc
Is the cutoff radius in DPD simulations with a fluid of constant density arbitrary? Isaac Gabriel Salib ∗ , Shimon Haber Technion – Israel Institute of Technology, Haifa 32000, Israel
a r t i c l e
i n f o
Article history: Received 1 September 2009 Received in revised form 22 February 2010 Accepted 15 March 2012 Available online 19 March 2012 Keywords: Coarse graining Dissipative Particle Dynamics Integration schemes Soft condensed matter
a b s t r a c t The DPD (Dissipative Particle Dynamics) method was used to simulate the flow of a fluid occupying the space inside a rotating infinite long cylinder. The fluid was initially at rest and rigid body rotation was expected to be reached. However, as the simulation proceeded, the DP’s (Dissipative Particles) began accumulating near the wall, leaving the cylinder’s core empty of DP’s, obviously, an aphysical solution. To solve this problem it was suggested to set the cutoff radius as a function of the DP’s number density. A conceptual explanation is suggested. © 2012 Elsevier B.V. All rights reserved.
1. Introduction There is a growing interest in simulating fluid flows in biological systems. Simulating blood flow in veins and arteries, or air laden particles flow in lungs is of great importance in medical and biotechnological researches. However, it is almost impossible to simulate these flows using classical simulation techniques: the fluids are complex, may include many species (such as blood flow consisting of red and white blood cells platelets, etc.) and bounded by intricate boundaries (such as soft tissues). Also nanoand micro-flows need a simulation technique which can capture their hydrodynamic behavior, taking in consideration their nanoand micro-scale complexity. A “coarse-grained” method for simulating fluid flows was suggested by Hoogerbrugge and Koelman [1,2]. This method, known as Dissipative Particle Dynamics (or DPD), is used to model a fluid by a finite number of particles, each of which represents a large number of molecules. The hydrodynamic behavior is predicted through tracing the particles, basically a Lagrangian approach. The force system exerted on the DP’s includes soft repulsive conservative, dissipative and random forces. The last two forces were termed by Hoogerbrugge and Koelman as the system’s “thermostat”. The dissipative term tends to “cool down” the system while the random term tends to “heat up” the system. The presence of both terms enables local momentum conservation and the prediction of hydrodynamic behavior in the macroscopic scale. The DPD model was modified by Español et al. [3,4]. They derived the Fokker–Planck
equation and showed how the DPD Langevin equations are related to statistical mechanics. Furthermore, in DPD, several no-slip boundary methods were employed [5]. Haber et al. [6] suggested a benchmark test of fluid flow between two rotating concentric cylinders. In their article, comparison was drawn between these methods and the bounceback reflection method (reflecting back the particles which pass the rigid boundary with an opposite velocity) was shown to yield the most accurate results. The two rotating concentric cylinders test was found to be a useful benchmark since this test case has an analytical solution and it represents a closed system where no artificial periodic boundary conditions are needed. Haber et al. [6] noted that the DP number density at the boundary was somewhat higher than between the cylinders. Surprisingly, longer runs revealed that after enough time steps, the DP’s tend to accumulate at the boundaries. The phenomenon was augmented at high cylinder rotating speeds for which fewer time steps were required. In order to investigate and resolve this issue, we used a simpler 2D model: a fluid, originally at rest, inside a rotating cylinder. The results are explained and a novel perspective of the DPD method is brought up, that is of particular interest for flows bounded by curved walls. 2. Method The basic and most common postulated DPD scheme introduced by Español et al. [3], takes the following form:
*
Corresponding author. E-mail addresses: [email protected] (I.G. Salib), [email protected] (S. Haber). 0010-4655/$ – see front matter © 2012 Elsevier B.V. All rights reserved. http://dx.doi.org/10.1016/j.cpc.2012.03.012
dri =
pi mi
dt
(1)
I.G. Salib, S. Haber / Computer Physics Communications 183 (2012) 1714–1718
1715
Fig. 1. DP’s accumulation near the boundary (each asterisk stands for a DP): Out of the 900 DP’s only a few remain inside the domain while the rest are at the boundary.
dpi =
fiCj (r i j ) +
j =i
+
The weight functions frequently used posses the following form:
−γ ω D (r i j )(vi j · rˆ i j )ˆri j dt
j =i
σ ω R (ri j ) dW i j rˆ i j
(2)
j =i
where ri j and vi j are the position and velocity of particle i relative to particle j, rˆ i j is the unit vector in the direction of ri j , pi and mi are the momentum and the mass of particle i respectively, and dt is the time increment. The conservative force is fiCj and the constants γ and σ are the coefficients of the dissipative and random forces respectively. The weight functions, ω D (r i j ) and ω R (ri j ), correspond to the soft matter interactions (all forces act up to a maximum distance range called cutoff radius rc ) and dW i j is an independent increment of the Weiner process for which Itô calculus rule implies:
dW ii dW j j = [δi j δi j + δi j δi j ] dt
(3)
From the corresponding Fokker–Planck equation, one gets Eq. (4):
ω2R (r ) = ω D (r ),
σ 2 = 2k B T γ
(4)
where k B is the Boltzmann’s constant and T is the equilibrium temperature. The conservative force commonly used is of the following form:
fiCj =
ai j (1 − r i j /rc )ˆri j
at r i j rc
0
at r i j > rc
(5)
where ai j is the maximum repulsion force constant. This latter parameter is taken to be uniform for all interparticle interactions and is simply replaced by a. It was established by Groot and Warren [7] through the fluid compressibility.
ω2R (r ) = ω D (r ) =
(1 − r i j /rc )2 at r i j < rc 0 at r i j > rc
(6)
The simulations were conducted using the bounce-back boundary condition. As the simulation proceeded, the DP’s molecules started to accumulate at the boundary (Fig. 1). This phenomenon occurred in case either the Euler or the Velocity Verlet [8] algorithms were used. It was also observed that the size of the time increment used had no effect on the DP accumulation at the boundary (except for the total number of time steps taken for it to occur). The DP’s which passed the boundary were brought back artificially by the bounce-back boundary conditions. These DP’s are usually reintroduced into the domain at positions close to the boundary (this may explain the fact that higher densities are obtained near the boundaries). But, as the simulation proceeded more DP’s were driven towards the boundaries. There are two mechanisms which, theoretically, hinder the DP’s boundary accumulation: 1. The fluctuating force which provides the DP with the random velocity and displacement which, by forcing the DP’s to intermingle randomly, could maintain a constant density over the whole domain. However, for simulations in which the rotational cylinder speed is much higher than the kinetic one, i.e. k B T /m, the random force is too small to achieve this balance. 2. The repulsive conservative force could also prevent the accumulation phenomenon. As noted earlier, since the DP number density at the boundary is higher than that in the rest of the domain, we would expect that the repulsive force exerted by
1716
I.G. Salib, S. Haber / Computer Physics Communications 183 (2012) 1714–1718
Fig. 2. Schematic sketch of two adjacent DP’s with interparticle distance that is (A) smaller than the sum of their radii and (B) larger than the sum of their radii.
the DP’s near the boundary would repel the particles back into the domain. This force is controlled by two parameters: cutoff radius and the force amplitude or a in Eq. (5). We tried first to adjust only the force amplitude. We found out that by increasing the conservative force amplitude, the accumulation phenomenon is just postponed, i.e. for the same magnitude of time increment, more time steps were needed to reach accumulation of DP’s near the wall. Apparently, only the adjustment of the cutoff radius and the magnitude of the conservative force could control the DP number density. 3. Theory The cutoff radius (coined rc ) designates the radius of influence, beyond which the DP cease to interact with other DP’s. This parameter is thought to be chosen arbitrarily. Different conclusions were drawn for the optimal cutoff radius, of which, the most acceptable, is that the cutoff radius has to define a spherical volume around the DP that contains at least three or more DP’s on average. It is worthy to note that setting the cutoff radius is of practical importance. For a given time increment, the smaller the cutoff radius, more time steps are needed to reach equilibrium. To our understanding, the cutoff radius also possesses a physical significance and not just a numerical one. The DP’s are, as noted earlier, lumps of molecules. Adjacent DP’s meddle. In the meddling zone, molecules bombard and momenta transfer. The weight function denotes that closer the DP’s, the larger is this zone and the greater is the resultant interactive forces (as illustrated in Fig. 2). Thus, the cutoff radius between two adjacent DP’s should equal the sum of their radii.
We shall limit our simulations here to constant density cases. The DP’s radii are uniform. Thus we expect that the cutoff radius is constant over the whole solution domain. 4. Results We simulated a fluid in a cylinder of radius R = 0.001 m that rotates with a speed of ω = 1 rad/s. Each DP is of m = 3.4907 × 10−12 kg mass (for water that means 1.16 × 1014 molecules in each DP). The simulation is conducted for T = 273 K temperature and 900 DP’s. Although, a wide range of numerical values for the dissipation coefficient were used, we choose here to present the results achieved for
γ m
=
1575 2π n2 rc5
μ
(7)
where μ is the shear viscosity and n is the initial uniform particle density. The conservative force amplitude used is taken as stated in Groot and Warren [7]. For the time increment, we followed the footsteps of Haber et al. [6]. In their choice of the time increment they took into account all the time scales in the test case rc / v T , m/γ and v T m/a, including, rc /ω R, that is linked to the boundary velocity, ω R, where R is the cylinder’s radius in terms of rc . The dimensionless time increment chosen is as follows:
dt = 10− f min 1, v T m/γ rc , v 2T m/arc , v T /ω R
where f is the numerical accuracy sought. For f we tried the following values: 2, 3, 4, and 5. Since the model used here is a 2D one, the area which each DP possesses is simply the cross section area of the cylinder divided
I.G. Salib, S. Haber / Computer Physics Communications 183 (2012) 1714–1718
1717
Fig. 3. DP’s inside the rotating cylinder after 100,000 time steps for rc = 0.0667 mm (each asterisk stands for a DP).
Fig. 4. DP’s inside the rotating cylinder after 100,000 time steps for rc = 0.4 mm (each asterisk stands for a DP).
by the number of the DP’s. Assuming that the DP is circular, the cutoff radius would be equal to twice this radius (the factor of 2 is the optimum for achieving the desired results):
R rc = 2 × √ = 6.67 × 10−5 m N
(8)
Simulation was conducted for 100,000 time steps. Fig. 3 illustrates the image of the DP’s at the end of the simulation. As we can see the DP’s are uniformly distributed and the hydrodynamic behavior is that of rigid body rotation as expected. In comparison, we have conducted a simulation with cutoff radius five times larger than that used in Eq. (7). The illustra-
1718
I.G. Salib, S. Haber / Computer Physics Communications 183 (2012) 1714–1718
tion of the DP’s locations after 100,000 time steps is presented in Fig. 4. As we can see the DP’s that were located near the boundary started accumulating at the boundary. As a result, unlike the uniform number density illustrated in Fig. 3, the DP number density is higher at the boundary. As the simulation proceeds, more and more DP’s accumulate at the wall and finally reach the distribution shown in Fig. 1. 5. Conclusions Most previous benchmark problems used to test the DPD method employed planar boundaries for Couette or Poiseuille flows where the flow does not experience centrifugal forces. Haber et al. [6] benchmark problem was the first to include a curved surface for which the continuous solution includes centrifugal effects. The DPD method is based on modeling lumps of fluid by DP particles. However, in the DPD model, the volume that each DP occupies is not set, only the DP mass is a priori assumed. Therefore, the density of the fluid may not be maintained. On using the rotating cylinder 2D simulation, this flaw of not keeping constant density is noted clearly; the DP number density could be much higher at the boundary if the cutoff radius is not properly picked. In order to avoid such a numerical artifact, we suggest that the cutoff radius is a measure of the DP’s extent. Adjusting the cutoff radius and the conservative force amplitude enable us to reach an almost constant density and prevent DP gathering near the boundary. Future research would attempt to relate the cutoff radius magnitude to the thermodynamical parameters in order to allow compressible flow simulations.
In this article, there is no account for using a different conservative force. Nevertheless, we also tried to seek a solution, by using a Lennard-Jones potential between the particles. However, when the interparticle distance tends to zero, this potential tends to infinity and so does the DP displacement, i.e. a very small time increment has to be used. In this case the DPD scheme is impractical. We also tried to replace the conservative force by a repulsive attractive spring force. In this case aggregates were formed throughout the domain, obviously an aphysical result. References [1] P.J. Hoogerbrugge, J.M.V.A. Koelman, Simulating microscopic hydrodynamic phenomena with dissipative particle dynamics, Europhys. Lett. 19 (1992) 155. [2] P.J. Hoogerbrugge, J.M.V.A. Koelman, Dynamic simulations of hard-sphere suspensions under steady shear, Europhys. Lett. 21 (1993) 363. [3] P. Español, P. Warren, Statistical mechanics of dissipative particle dynamics, Europhys. Lett. 30 (1995) 191. [4] P. Español, Hydrodynamics from dissipative particle dynamics, Phys. Rev. E 52 (1995) 1734. [5] D.C. Visser, H.C.J. Hoefsloot, P.D. Iedema, Comprehensive boundary method for solid walls in dissipative particle dynamics, J. Comput. Phys. 205 (2005) 626; I.V. Pivkin, G.E. Karniadakis, A new method to impose no-slip boundary conditions in dissipative particle dynamics, J. Comput. Phys. 207 (2005) 114; M. Revenga, I. Zúñiga, P. Español, Boundary model in DPD, Int. J. Mod. Phys. C 60 (1998) 1319. [6] S. Haber, N. Filipovic, M. Kojic, A. Tsuda, Dissipative particle dynamics simulation of flow generated by two rotating concentric cylinders. Part I: Boundary conditions, Phys. Rev. E 74 (2006) 046701. [7] R.D. Groot, P.B. Warren, Dissipative particle dynamics: Bridging the gap between atomistic and mesoscopic simulation, J. Chem. Phys. 107 (1997) 4423. [8] I. Pagonabarragga, M.H.J. Hagen, D. Frenkel, Self-consistent dissipative particle dynamics algorithm, Europhys. Lett. 42 (1998) 377.