Simulation of individual polymer chains and polymer solutions with smoothed dissipative particle dynamics

: In an earlier work (Litvinov et al. , Phys.Rev.E 77, 066703 (2008)), a model for a polymer molecule in solution based on the smoothed dissipative particle dynamics method (SDPD) has been presented. In the present paper, we show that the model can be extended to three-dimensional situations and simulate effectively diluted and concentrated polymer solutions. For an isolated suspended polymer, calculated static and dynamic properties agree well with previous numerical studies and theoretical predictions based on the Zimm model. This implies that hydrodynamic interactions are fully developed and correctly reproduced under the current simulated conditions. Simulations of polymer solutions and melts are also performed using a reverse Poiseuille ﬂow setup. The resulting steady rheological properties (viscosity, normal stress coefﬁcients) are extracted from the simulations and the results are compared with the previous numerical studies, showing good results.


Introduction
The simulation of polymer dynamics represents a long-standing problem in computational rheology.The main challenges are related to the choice of an appropriate fluid-polymer interaction model as well as of the simulation method itself.An important fact which is generally known [1] is that, for a diluted solution, the hydrodynamic interaction (HI) has strong effects on the dynamics of the polymer.Experiments have confirmed that, under diluted conditions, polymer dynamics is best described by the Zimm model [2], where the effects of HI are taken into account by using the Oseen tensor formulation, and the solvent is assumed to behave as an incompressible Stokes flow.
A straightforward approach to take HI into account is to simulate solvents and polymers using molecular dynamics (MD).Such studies were done and provided valuable information for the theory of polymer dynamics [3].An obvious disadvantage of this approach is the computational cost; in fact, due to limited numerical resources, it is not possible to achieve the relevant length and time scales of experiments.Some modifications, such as coarse-grained molecular dynamics (CMD) [4], have been proposed in the past to reduce the required computational time, but the problem still remains.
More effective schemes include a stochastic thermal force explicitly applied on the suspended bead while the effect of HI among them is implicitly taken into account via Oseen tensor-based models or higher-order expansions of the incompressible Stokes equations.Afterwards, a resulting set of algebraical equations needs to be inverted.These approaches include, for example, Brownian dynamics (BD) [5][6][7][8][9] and Stokesian dynamics (SD) [10,11].Their main advantage is related to the good scalability of the number of equations with the number of suspended particles, allowing for the simulation of large systems [12].On the other hand, modelling complex external geometries is not an easy task in implicit-solvent methods and requires essential modification of the standard framework [13][14][15].
To remedy these problems, various researchers have focused in the past on mesoscopic methods which employ numerically effective coarse-grained models retaining the relevant hydrodynamics modes.Examples are multiparticle collision dynamics [16] and lattice Boltzmann methods (LBM) [17,18].Due to the presence of an underlying grid for the computation of the hydrodynamic interactions, LBM is computationally very effective; however, it generally requires hybrid schemes to describe the dynamics of immersed structures.Using different numerical approaches for solvent and suspended particles in fact requires special care in fluid-structure interaction models in order to ensure exact transfer of momentum and more generally formal thermodynamic consistency of the coupled-scheme.
A mesoscopic particle method called dissipative particle dynamics (DPD) which uses pairwise soft conservative, dissipative and random forces, has gained increasing popularity for its ability of simulating much larger physical length and time scales [19].The DPD model shows good performance on investigating the dynamic properties of polymers and conformation of DNA in a microchannel [20,21].A remarkable phenomenon called cross-stream migration of polymers was found recently in a non-uniform shear flow using the DPD model, and the effects of Reynolds number, Schmidt number, polymer chain length and polymer concentration were investigated [22][23][24][25].In addition, the so-called "reverse Poiseuille flow" with a DPD polymer model is a convenient virtual rheometer for the acquisition of steady state rheological data [26].
In this paper, a new type of DPD method, called smoothed dissipative particles dynamics (SDPD) [27][28][29], is used to investigate the static and dynamic behavior of polymers.This model is a mesoscopic extension of a popular particle-based method applied to continuum flow, i.e., smoothed particle hydrodynamics (SPH) [30].Due to its connection to SPH, SDPD allows using an arbitrary equation of state for the pressure and to specify transport coefficients directly without the need to perform kinetic analysis or preliminary runs.Thermal fluctuations can be included in a physically consistent way, that is, in agreement with a fluctuation-dissipation theorem (FDT) and with basic properties of statistical mechanics where fluctuations of hydrodynamic variables increase naturally by reducing the size of the problem down to the mesoscopic scale [29].These properties of SDPD make it a suitable tool to simulate physically relevant length and time scales as seen in experiments.If thermal fluctuations are not considered, SDPD reduces to a version of SPH where pressure/velocity profiles for incompressible flows are correctly reproduced [31,32].
In this paper, the thermodynamically-consistent SDPD model proposed in Reference [33] for a polymer molecule in solution is extended to three-dimensional situations allowing for effective simulations of polymer solutions under dilute and concentrated conditions.The paper is organized as follows: In Section 2, the SDPD model for the solvent and flexible polymer chain is outlined.In Section 3, the results for the static and dynamic properties of the simulated single polymer are presented.Steady shear-rate rheological properties of concentrated polymer solutions and polymer melts are provided and discussed in Section 4. Finally, in Section 5, the results are summarized and conclusions are drawn.

Mesoscopic Modelling of the Solvent
For the solvent, we assume an isothermal Newtonian model.This is described by the following Navier-Stokes equations written in a Lagrangian framework where d/dt is the Lagrangian derivative, ρ is the local fluid density, v the velocity, p the pressure and η the dynamic viscosity.The SPH discretization of the Navier-Stokes equations is obtained by considering a discrete number of fluid particles distributed over the domain and reformulating Equation (1) in terms of pairwise interparticle forces.The set of ordinary differential equations for the particle positions and velocities reads where the mass density is calculated as Here, m i is the mass of a particle i, r ij are the normalized vector and distance from particle i to particle j.To close the system, an equation of state for p i must be specified.In this work the following form is used: where p 0 , ρ 0 , b and γ are parameters which can be chosen based on a scale analysis to minimize density variations, i.e., by choosing the corresponding fluid speed of sound c s = p 0 γ/ρ 0 to be sufficiently large.Since a stiff equation of state is used, usually γ = 7, penetration between particles is prevented.
In the SDPD formulation [27,29,34], Equation (2) represents the deterministic part of the particle dynamics.Using the GENERIC formalism [35,36], thermal fluctuations on the fluid variables can be taken into account by postulating mass and momentum fluctuations to be where dW ij is the traceless symmetric part of independent increment of a Wiener process and B ij is In this way, the FDT is exactly satisfied at the discrete particle level [27].Some important properties of this particular set of equations are: (i) the total mass and momentum are exactly conserved; (ii) the linear momentum is locally conserved due to the anti-symmetric property of the discretization; (iii) because the formalism has been built within the GENERIC framework [35,36], for non-isothermal situations, it is possible to define an evolution equation for a particle internal energy and/or entropy such that the system conserves total energy, and the total entropy is a monotonically increasing function of time [37].

Mechanical Modelling of the Polymer Chain
In this work, a coarse-grain spring-bead chain model for the polymer molecule is used.The polymer beads have the same mass and interact with each other in the same way as the fluid particles.This assumption is justified by the fact that a polymer bead is a coarse-grained system obtained from a volume occupied by a large number of solvent molecules and a much smaller number of polymer monomers [33].
On the other hand, the effect of the chemical bonds between monomers is taken into account through an additional potential acting only between neighboring polymer beads.Many possible choices have been investigated [38] such as harmonic chain, Gaussian chain, Rouse chain, etc.A convenient model to simulate a realistic polymer molecule is the finite extensible non-linear elastic (FENE) potential [39] which is adopted in this work and reads where r is the distance between neighboring beads, R 0 is the maximum allowed spring length and k is the spring constant.Note that no other additional forces are introduced.Since the parameter R 0 is chosen as 1.4 for the average particle distance, the crossing of springs is prevented and polymer topology is preserved.The average SDPD particle distance is estimated based on the initial lattice positions and, being the solvent quasi-incompressible, remains approximately constant during flow.
In addition, the use of the stiff equation of state (Equation ( 3)) prevents beads from overlapping and therefore excluded volume effects are fully accounted for.Note that the pressure terms appearing in the discretized hydrodynamic equations produce already repulsive forces between approaching particles and therefore there is no need to introduce short-range Lennard-Jones forces between beads.
In the present model, since the solvent is also represented by particles, the hydrodynamic interactions are included explicitly by means of the interacting forces acting between the solvent particles and the polymer beads.These are the conservative forces caused by the pressure gradient, the viscous force caused by velocity gradients and the random force produced by thermal fluctuations.The final momentum equation for the polymer bead i reads, therefore, where m i is the mass of the bead i, is the total sum of the resulting SDPD forces acting on particle i given in Equations ( 2)-( 4), and i 1 , i 2 are the next two neighboring beads of i in the polymer chain which interact via the non-hydrodynamic forces derived by the potential Equation ( 6), i.e., Note that an important advantage of this particle-based approach is strict preservation of Galilean invariance, which can be violated in the presence of a computational lattice [40].

Simulation of a Single Polymer
In this section, the static and dynamic properties of a single polymer suspended in the Newtonian medium are investigated.A cubic periodic domain with total number of fluid particles 44 × 44 × 44 = 85, 184 is used for all the simulations in this section.The simulated single polymers have variable lengths with number of beads N = 5, 10, 20, 40, 60, 80 and 100.The shortest polymer N = 5 was not considered for some analysis since most of the theoretical models for the polymer are only valid in the asymptotic limit of a long chain.
The polymer and solvent particles are placed initially on a regular grid and a preliminary simulation is run until the system is fully thermalized.Afterward, the resulting configurations of the system are chosen as initial conditions for the subsequent simulations each running with a typical number of time steps equal to 10 5 .In order to obtain accurate statistics, the final results are averaged from 10 parallel runs.

Static Properties
The equilibrium static properties of the polymer are characterized by the mean square radius of gyration of the polymer chain.The radius of gyration is defined as where r ij = r i − r j is the distance between beads.The end-to-end distance is defined as where N denotes the number of beads of the polymer chain.The static exponent ν relates the characteristics size of polymer and the number of the beads by The calculated mean square radius of gyration is shown in Figure 1.The estimated static exponent is around 5% above the asymptotic value for a self-avoiding walk ν ≈ 0.588.
Another way to obtain the static exponent is to use the static structure factor This method has certain advantages since, even for a single polymer, it probes different length scales.In the scaling regime R −1 G k a −1 0 , where a 0 is the characteristic length of the neighboring beads distance, the relationship S(k) ∝ k −1/ν (13) holds.Figure 2 shows static structure factors for polymers with the number of beads N = 20, 30, 40, 60, 80.The results for the polymers with a small number of beads N = 5, 10 show deviations from the overall trend and are not presented.This can be explained by the fact that for short chains the effect of the polymer ends is non negligible.By fitting the power-law in logarithmic scale, one can obtain the values of ν listed in Table 1.As a comparison, with the renormalization group theory and Monte Carlo simulation, an accurate value of ν = 0.588 has been obtained [41,42].Our results are in good agreement with this prediction as well as with previous CMD, BD and DPD simulations [17,20,43,44].Since the obtained static exponent is typical for a self-avoiding chain, one may conclude that the excluded volume effects between SDPD particles are taken into account correctly.

Diffusion Coefficient
One of the most important dynamic properties of a polymer chain is the diffusion coefficient (D).This can be estimated from the mean-square displacement of the polymer center of mass: The relation between g/(6t) and 1/t was plotted and extrapolated to 1/t → 0 where g/(6t) → D. The resulting power-law dependence of D versus number of beads is shown in Figure 3.The fitted slope ν = 0.58 ± 0.05 is in good agreement with the prediction of the Zimm theory.Note that the same value of the exponent was obtained in the DPD simulations [44] by extrapolating the data from a much smaller size domain to the thermodynamic limit.

Longest Polymer Relaxation Time
The longest relaxation time of the polymer represents another important dynamical feature and is obtained by fitting the auto-correlation function of the radius of gyration where the brackets indicate averages over all parallel simulations at various time origins t 0 .For C(t) within the range [0.1 − 1.0] the curves were approximated by the expression C G (t) = exp(−t/τ Z ) and the fitted values of τ Z are shown in the Figure 4.One can verify whether the HI effects have been introduced in the current simulations by fitting the scaling of τ Z on polymer length with the static exponent.For the Zimm model, one expects τ Z ∝ N 3ν , while for the Rouse model τ Z ∝ N 1+2ν should hold.In Figure 4, the predicted static exponent resulting from the Zimm scaling is 0.57 ± 0.03 which is in good agreement with the independent value obtained above by the static analysis.On the contrary, application of the Rouse model delivers a much smaller value.Agreement with the Zimm model is consistent with the full description of hydrodynamic interactions between beads within the chain which is not taken into account in the Rouse model.

Rouse Modes
The Rouse coordinates for a discrete polymer are defined as where p = 0, 1, 2, . . .are the number of Rouse modes and r n is the position of bead n.R p are the eigenfunctions of the Gaussian model which provide important information on the static and dynamic behavior.Note that R 0 is simply the polymer center of mass whereas for p > 0, R p = 0, and R p (t) characterize the internal motion of the polymer on various characteristic lengths.
To test whether the HI is correctly reproduced by the SDPD model, we have analyzed relaxation times τ p of the Rouse modes.Figure 5 shows the scaling plot of τ p with N/p.The collapse of the curves for different polymers and Rouse modes is clear and the scaling exponent is found to be 1.78 ± 0.05, in good agreement with the prediction of the Zimm model (3ν ≈ 1.74) and substantially below the prediction of the Rouse model (1 + 2ν ≈ 2.16).When τ p is plotted against (p R p ) 1/2 , as shown in Figure 6, the fitted scaling exponent is in the range 2.79 ± 0.07, which is only slightly below the prediction of the Zimm model (3.0), whereas the Rouse model predicts a much larger value 2 + 1/ν ≈ 3.72.The last observation can be considered as an independent test for the presence and accurate description of HI since the prediction of the Zimm scaling (3.0) is not a consequence of the correct prediction of the static exponent [44].

Polymer Melt and Solution
The polymer melt is modelled as an ensemble of flexible polymer chains without solvent particles.The corresponding rheological properties are investigated by using the so-called "reverse Poiseuille flow" (RPF) which consists of two parallel Poiseuille flows driven by uniform body forces of equal magnitude but acting in opposite directions [26].DPD simulations using RPF have been successfully used in diluted polymer solutions [23], colloidal suspensions [45] and polymer fluids [26].
In the current simulation, the RPF is driven by a constant force F x acting on every SDPD particle in the x-direction defined as where L y is the size of the domain in the y−direction.A cubic periodic domain is used with a total number of fluid particles 29 × 73 × 15 = 31, 755.The domain size is L x = 7.21 × 10 −4 , L y = 1.80 × 10 −3 and L z = 3.71 × 10 −4 .The timestep is set to ∆t = 5.2 × 10 −6 and all simulations are run at least for 10 6 timesteps.For sake of completeness, the full set of parameters used in the simulations are listed in Table 2 based on SI units.For a pure solvent simulation under the considered driven force, the estimated Reynolds number based on maximal fluid velocity is Re = 5.8 and the Schmidt number is Sc = 770 which is in the same order as Schmidt numbers for real fluids [46].Figure 7 shows a typical snapshot of the polymer system.
For power-law models of non-Newtonian fluids, the RPF velocity profile reads [26] (19) where L = L y is the length of the simulation domain in the y−direction, p is the power-law index (e.g., p = 1 corresponds to a constant-viscosity Newtonian model), κ is the power-law shear stress coefficient and n is the number density. Figure 8 shows the velocity profile across the channel for polymer melts of 2, 5, 25-bead chains and, as a comparison, for a pure solvent.The solvent velocity profile is parabolic (p = 1), whereas the melt velocity profiles are increasingly non-Newtonian (for increasing number of bead chains) and are in good agreement with those reported for DPD simulations, with p being slightly higher for SDPD [26].
Table 2. Parameters used for the simulation of polymer melts.Here ∆x is the initial smoothed dissipative particle dynamics method (SDPD) particle spacing and c s is sound speed.Higher power indices indicate that the simulated melts are less viscous and exhibit a slightly weaker shear-thinning behavior compared to those reported in the DPD simulations.One possible reason for this discrepancy is that the value of parameters, spring constant k and maximum allowed spring length R 0 , employed in the FENE bond differs from those used in DPD.
Distribution of polymer beads across the domain is investigated in Figure 9 where a uniform bead density field is observed for the 25-bead melt.This is an essential phenomenon for the pure-chain melt as a result of incompressibility.As a further test, next we measure stress field present in the system [48].In the RPF, an off-diagonal component of the shear stress τ xy (y) = f n(L/4 − y) for y ∈ [0, L/2] and τ xy (y) = f n(y − 3L/4) for y ∈ [L/2, L] should be recovered.Figure 10 reports the imposed and the calculated normalized shear-stress distribution for the melt with N = 25 bead chains for comparison, showing good agreement.Figure 11 shows also the cross channel distribution of normalized normal-stress components whereas the normalized normal stress difference is shown in Figure 12.A coarser representation is used for normal stress and normal stress difference, because their averages appear to be more noisy than that for the velocity and the shear stress.The ability of the SDPD model for simulating a concentrated polymer solution is investigated by using a 50% concentration of 25-bead chains suspended in a solvent of SDPD fluid particles.The driving constant force is the same as in the case of a melt.The viscosity of the solution is smaller as suggested by the higher maximal velocity value (0.037) compared to the value (0.029) obtained for a melt composed by the same 25-bead chains, as shown in Figure 13.Moreover, the larger value of power index of the velocity profile (p = 0.9 with respect to p = 0.77 corresponding to the melt) suggests that the solution has an increasingly Newtonian-like behavior, with viscosity being only slightly shear-rate dependent.
Finally, Figure 14 shows the normalized bead density of the solution.The two large density peaks occur around at y = 1/4L and y = 3/4L positions where the velocity gradients are zero with maximal density around 1.3.This is a manifestation of the polymer cross-stream migration phenomenon observed in experiments [49] and simulations [24,26,50].

Conclusions
In this study, the static and dynamic properties of single polymer and polymer solutions were investigated using the smoothed dissipative particle dynamics (SDPD) method.By comparing our results with theoretical models and previous studies, the following conclusions can be drawn: (i) The obtained polymer static properties of radius of gyration and static structure factor for a single polymer molecule show the presence of excluded volume effects; (ii) The predicted polymer dynamics of diffusion coefficient, longest relaxation time and relaxation times of Rouse modes follow the Zimm model closely, which includes hydrodynamic interaction effects; (iii) The SDPD model captures the steady shear rheological properties of typical polymer melts; (iv) For the polymer solution with 50% polymer concentration, the method qualitatively reproduces the cross-stream migration phenomena reported in previous experimental and numerical works.Note that similar simulations with original DPD [44] have also recovered several static and dynamic properties of dilute polymer in solution without any special treatment of the possible non-physical penetration and crossing of polymer bonds.However, since the excluded volume effects are fully taken into account with the present SDPD method, we may expect more accurate predictions when these effects became important, as in the case of polymer-wall [21] interaction under the flow or transition of the polymer conformation [9], which are currently under study.

Figure 1 . 2 G
Figure 1.Scaling of the radius of gyration R 2 G for chain lengths corresponding to N = 5, 10, 20, 40, 60, 80 beads.The solid line represents the best fit of the static exponent.

Figure 2 .
Figure 2. Equilibrium static structure factor vs. wave vector magnitude k for polymers with length N = 20, 40, 60, 80.All the curves collapse on a master line within the region of R −1 G < k < a −1 0 .

Figure 3 .
Figure 3. Log-log plot of the diffusion coefficient D as a function of the chain length N. The solid line represents the best fit with the theoretical power law.

Figure 4 .
Figure 4. Log-log plot of the longest relaxation time, computed from the autocorrelation of the radius of gyration, as a function of the chain length N = 10, 20, 40, 60, 80.The solid line is the best fit to determine the static exponent according to the Zimm model.

Figure 5 .
Figure 5. Log-log plot of the relaxation time of Rouse modes τ p vs N/p for N = 5, 10, 20, 40, 60, 80, 100.The solid line is the best fit to determine the scaling exponent.

2 . 5 × 12 Figure 7 .Figure 8 .
Figure 7. Snapshot of the polymer melt simulation: x is a direction of the body force and y is the direction of the velocity gradient.The visualization was done with the Visual Molecular Dynamics (VMD) program [47].

Figure 9 .
Figure 9. Normalized bead density profile for the melt composed by 25-bead chains.

Figure 10 .
Figure 10.The calculated and imposed normalized shear-stress distribution for the melt composed by 25-bead chains.

Figure 11 .Figure 12 .
Figure 11.Normalized stress distribution for the melt of N = 25 bead chain.

Figure 13 .Figure 14 .
Figure 13.Velocity profile for the melt of N = 25 bead chain for bulk concentration 50%.Power-law index p is 0.90.

Table 1 .
Estimated static exponent from the fit of the linear part of S(k) within the region ofR −1 G < k < a −1 0 .