Abstract
Despite extensive area of applications, simulation of complex wall bounded problems or any deformable boundary is still a challenge in a Dissipative Particle Dynamics simulation. This limitation is rooted in the soft force nature of DPD and the fact that we need to use an antipenetration model for escaped particles. In the present paper, we propose a new model of antipenetration which preserves the conservation of linear momentum on the boundaries and enables us to simulate complex and flexible boundaries. Finally by performing numerical simulations, we demonstrate the validity of our new model.
1. Introduction
Dissipative Particle Dynamics is a particle-based method for simulation of complex fluids at mesoscopic scales. DPD is a coarse-grained version of molecular dynamics (MD) in which DPD particles formed by a clusters of molecules. The coarse graining removes the hard core nature of interparticle forces and DPD particles interact with each other through a set of soft-core forces. DPD was introduced by Hoogerbrugge and Koelman, in 1992 [1] for predicting the behavior of isothermal complex fluids. The fluctuation-dissipation theorem by Español and Warren [2] provides the necessary condition for a system to reach thermal equilibrium state, and a revision of the method by Groot and Warren [3] has been developed for hydrodynamic problems. The method has proven to be capable of simulating a wide variety of systems, from simple fluids hydrodynamics [4–6] to colloidal and dilute polymer solutions [7, 8].
Because of the soft core nature of DPD forces, the length and time scales in DPD are much larger than MD, and this advantage makes DPD capable of simulation of larger length scales within larger time interval. However, because of these soft forces, in wall bounded domains, DPD forces cannot hold the fluid particles inside the domain. The absence of a general mechanism to prevent particle penetration into the solid boundary is one of the unresolved issues that extremely narrow down the range of problems that can be studied by the method.
In recent years, the major treatment for wall bounded problems has been developed by distribution of freezing particles in one or a few layers with combination of proper reflection procedure for escaped particles. To address the most important publications on this treatment of boundary conditions, we should mention the early works of Revenga et al. [9, 10], where they modeled the solid boundary with a single layer of stationary particles, and since the effective computed forces from the wall to fluid particles are not sufficient to prevent wall penetration, reflections were found necessary to reflect the escaped particles back into the domain. In [9], the effects of three different wall reflection methods, including specular, Maxwellian, and bounce-back reflections, were studied and their effectiveness in producing no-slip condition on walls was discussed. However, their approach has some essential problems. It cannot ensure the nonslip conditions in all cases, introducing conservative force leads to large density fluctuation and using the reflection mechanism narrows down the applications of method into simple cases.
The most representative works published so far following the Revenga et al. work [9] have been used different approaches to deal with the nonslip boundary conditions as well as density fluctuations near the wall [11–14]. However, since in all of the methods, a reflection mechanism is applied to reflect the escaped particles back into the domain, all of the methods have major difficulties for simulation of complex geometries as well as flexible boundary conditions.
To the best knowledge of authors, no systematic model has addressed the treatment of complex geometries or flexible boundary problems and, therefore, finding a general procedure for these problems still remains a challenge. In the present paper, we focus on this issue and propose a systematic method for imposing proper boundary conditions at the confined geometries with stationary or flexible boundaries. Here we start with a brief review on the method. Then a model is presented to impose the wall boundary conditions at the complex geometries. Finally, the proposed model is validated by presenting numerical results for a set of benchmark problems.
2. Dissipative Particle Dynamics
The motion of each particle in DPD is governed by Newton’s second law. The particle is characterized by its position , velocity , and mass , and the motion of this particle is governed by
In DPD, the inter-particle force is separated into three pair-wise contributions from conservative, dissipative, and random forces. All the forces between particles and vanish beyond some cut-off radius [2]:
Conservative force is defined as below: which is maximum repulsion force and is the unit displacement vector from particle to particle . Dissipative and random forces are, respectively,
In order for the DPD system to have a well-defined equilibrium state obeying Boltzmann statistics, the random and dissipative forces should follow the fluctuation-dissipation theorem [2]:
Substituting these three forces into the equation of motion (1), the particles’ time evolution can be obtained by means of a numerical integration.
3. Implementation of New Boundary Conditions
In this section, we propose a systematic method for imposing proper boundary conditions at the stationary or flexible boundaries of confined geometries. For this, we start our discussion for a general wall distribution, redefine the required properties of the wall, and develop a method to apply desirable boundary conditions. For any Lagrangian particle based method, one can classify the following properties for any correct wall boundary conditions:(1)the ability to impose accurate no-slip conditions, with the least fluctuations in flow variables, near the boundary,(2)having correct antipenetration mechanism, which calculates the exact force on the wall particles from the fluids.
The first issue has been discussed thoroughly in many literatures, such as [9–14]. For the second condition, we need to introduce a new antipenetration model. In this section, we focus on this problem and we try to develop a new method for the antipenetration procedure.
Due to the soft core nature of DPD forces, the penetration of fluid particles into the solid region seems unavoidable, and hence in the majority of DPD implementation, antipenetration mechanisms play an essential role in implementation of boundary conditions. As mentioned before, the conventional method is based on reflecting the escaped particles. However, for a general, complex, and flexible boundary, the entire proposed wall boundary conditions that are reviewed in the literature have two major problems: first, the conservation of linear momentum will be impaired by reflecting a penetrating particle back into the domain. While we correct the position and velocity of escaped particles, reactions of these changes in the momentum are not satisfied on the boundary particles. For the problems with fixed boundary, this has little effect son the simulation. However for the flexible boundaries, since the motion of wall is controlled by the fluid forces, we need a method to calculate these forces.
The second problem is that finding the penetrated particle by these methods is limited to simple geometries and for complex boundaries, it is difficult to recognize escaped particles and simply obtain their reflected positions.
To overcome these difficulties, we construct a surface gird on the boundary to mark the escaped particles, and then by properly adjusting the forces between an escaped particle and the wall particles, we try to force them back into the domain. The main steps of this procedure are as follows.(1)At the start of simulation, a triangular surface grid is generated on the solid wall using the DPD particles as its grid nodes. Each element of this grid consists of three DPD solid particles at positions of .(2)The normal vector of each surface element is calculated using its particles position by (3)For each of wall particles , the normal vector is defined as the average of normal vectors of all elements sharing that particle as a vertex.(4)For each fluid particles , that is, in the neighborhood of a wall particle , that is, , we check the possible penetration by calculating the internal cross of (see Figure 1). For a positive value of , the particle penetrated into the wall and a negative value indicates that the particle is within the domain: (5)After finding the escaped particle and its wall neighbor particle, instead of usual repulsive conservative force (3), an attractive antipenetration force is applied on this pair of particles. This force is acting along the centers of escaped particle and wall particle . The force should be defined in a way that escaped particle gets back into the domain at the end of the next time step.
The success or failure of this method is dependent on the correct definition of . In the present work, we propose the following form for the antipenetration force between the penetrated fluid particles and wall particle : where is the effective repulsive coefficient between wall and fluid particles and is a coefficient to adjust the magnitude of the force. The value of that gives good results in the simulations will be determined by numerical investigation in the next section.
4. Numerical Results
In this section, we explore the validity of the proposed model by presenting some numerical results. To find the proper value for antipenetration coefficient , we start our simulations by modeling a planar poiseuille flow. To examine the model for more complex geometries, we perform some simulations on poiseuille flow in a circular pipe. Finally, to ensure the validity of the model for flexible boundaries, we present the simulation of an elastic membrane under sudden internal pressure from a DPD fluid.
4.1. Plane Poiseuille Flow
All of DPD simulations for plane poiseuille flow were performed in a cubic cell with a size of . We use a conventional two-layer wall configuration that is widely used in the literature. In this configuration, a wall is constructed from two layers of particles distributed on a regular lattice with distance ; for example, see [13]. The following parameters were used in all of the cases: , , , and . The flow is generated by using a body force of on each particle.
For the first set of simulations, conservative force coefficient between the fluid-wall pair of particles has been set by the method introduced in [13]: . The simulation has been performed for different value of . The results show that for , there are some escaped fluid particles from the wall boundaries and the attractive wall particle forces are not sufficient to get them back into the domain. On the other hand, for very large value of antipenetration coefficient , the simulation broke down because of particle clustering near the boundaries.
In Figure 2, the non-dimensional density profile for different values of is plotted. As can be seen in this picture in all of the cases, fluid density outside the wall is not zero, but it does not mean that fluid particles are able to escape from the wall since they are pushed back to the domain in the next time step. Figure 2 shows that increasing leads to a reduction in the number of fluid particles penetrated into the wall.
Non-dimensional velocity profiles for different values of are displayed in Figure 3. The figure shows that the no-slip condition is not satisfied at the wall for some cases. Previous studies revealed that the original bounce-back reflection methods help to keep the no-slip condition [10] at the wall boundaries. By increasing the value of , slip velocity decreases, and as it can be seen in Figure 3, in this case, for a value of , a no-slip condition is set.
4.2. Poiseuille Flow in a Cylindrical Pipe
Next, we consider flow in a pipe with circular cross section. The pipe is modeled using two layers of wall particles, with the distance of from each other (see Figure 4(a)); simulation parameters are , , , , , and pipe radius is . The flow is generated by using a body force of on each particle.
(a)
(b)
Density profile for this case is plotted in Figure 4(b). Also in this figure, we compared the analytical velocity profile with the DPD results. A good agreement is obtained between the two solutions.
4.3. Spherical Membrane under Sudden Internal Pressure
Here, we try to simulate the deformation of an elastic membrane due to the internal pressure loading, where the internal pressure is applied from DPD fluid particles. The membrane is flexible and its position and shape can be changed in every step. To construct a membrane, we use uniformly distributed DPD particles on a sphere with the radius R0, where they linked to each other by triangular spring elements (Figure 5).
The complete dynamical description of such a system can be achieved by using the Helmholtz free energy of the system, which includes the linear link energy, energy of in-plane deformation of triangular elements, and bending energy between two adjacent elements (see Figure 6).
However, in the present problem, the uniform pressure on the membrane can only cause radial expansion or contraction. For this condition, the bending or planar deformation energies do not play any role in the final result. Here, only a linear spring element is used to form the links between the membrane particles as where is initial distance between particles and .
By this definition of membrane forces, we know that the change in the radius of membrane should be proportional to the internal pressure; that is, . In the DPD simulation, if the boundary forces are calculated correctly, we should reach a linear relation between radius of membrane and the pressure of its internal fluid.
We performed our simulations on a membrane constructed with 2966 particles on a sphere with initial radius. The membrane interior is filled with 13026 fluid particles, and the constant parameters in all of the simulations are set to , , , and .
To generate different internal pressure, different values are used for the conservative force coefficient . After sufficient time, the radius of sphere is measured, and by numerical investigation, the equilibrium values of temperature and density will be obtained. After that, we simulate a DPD system with the same density, temperature, and DPD parameters, and we calculate the pressure of the system from numerical investigation. The pressure calculation method is widely mentioned in the literature, such as [3]. Figure 7 shows the variation of sphere radius against the internal pressure. The linear relation between the two parameters and the fact that during the simulations all of the fluid particles are kept inside the membrane shows the validity of our new boundary models.
5. Conclusion
We propose a new model for identifying and reflecting escaped particles in DPD. The presented model automatically identifies the escaped particles, and unlike the traditional reflection mechanisms, it has the ability of preserving linear momentum. These advantages make our model suitable for simulation of moving boundary problems. The validity of our model has been verified by some numerical simulations.
Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.