# Polymer rheology simulations at the meso- and macroscopic scale

###### Abstract

We show that simulations of polymer rheology at a fluctuating mesoscopic scale and at the macroscopic scale where flow instabilities occur can be achieved at the same time with dissipative particle dynamics (DPD) technique. We model the visco-elasticity of polymer liquids by introducing a finite fraction of dumbbells in the standard DPD fluid. The stretching and tumbling statistics of these dumbbells is in agreement with what is known for isolated polymers in shear flows. At the same time, the model exhibits behaviour reminiscent of drag reduction in the turbulent state: as the polymer fraction increases, the onset of turbulence in plane Couette flow is pushed to higher Reynolds numbers. The method opens up the possibility to model nontrivial rheological conditions with ensuing coarse grained polymer statistics.

## Introduction

The plethora of intriguing phenomena that can be observed in flows of complex fluids is attracting increasing attention among physicists. The study of polymer fluids and melts has since long occupied a central position within this field [1, 2]. Typically, one is either interested in the response properties of polymeric fluids or the way they flow. The former set of problems concerns viscoelasticity of polymers and biomaterials, behaviour of single long flexible molecules in flows, etcetera [2, 3, 4]. The latter usually focuses on the differences between macroscopic flows of Newtonian and polymeric fluids in the same geometry. One such striking example is the recently discovered chaotic flows of dilute polymer solutions at very small Reynolds numbers – the so-called purely elastic turbulence [5, 6]. Another is the phenomenon of drag-reduction – the observation that even minute amounts of polymer can significantly suppress Newtonian turbulence and hence reduce turbulent drag [7]. In this Letter we develop a mesoscopic simulation method that is capable of addressing both classes of problems.

Simulation methods of polymers essentially fall into two classes. On one hand there are many mesoscopic coarse-grained approaches that have mainly been developed to study the thermo- and hydrodynamic properties of polymers in equilibrium and in weakly non-equilibrium situations such as an imposed small shear. Such results can, e.g., be compared with recent experimental results for the orientation statistics of single DNA molecules in solution [3]. However, these models typically cannot be scaled up to simulate flow at macroscopic rheological scales. On the other hand, numerical studies of polymer rheology at macroscopically relevant scales are based almost exclusively on numerical implementations of continuum constitutive equations like the Oldroyd-B or FENE-P models [1, 2]. By their very nature, these deterministic approaches only give the average behaviour, so they cannot give insight into the coupling between the macroscopic flow behaviour and the molecular properties. In this Letter, we for the first time bridge the gap between these two approaches and scales by introducing a coarse-grained Dissipative Particle Dynamics (DPD) model [8, 9, 10] for visco-elastic flows which is a solution of elastic springs (dumbbells). Unlike the previous DPD simulations of polymeric fluids [11], we do not attempt to resolve internal dynamics of long polymer molecules. Instead we view the dumbbells as collective elastic degrees of freedom (normal modes) that render the solution viscoelastic. We show that the model is powerful enough to exhibit both stretching and tumbling of dumbbells at mesoscopic scales, and proper flow behaviour at hydrodynamic scales like the dramatic polymer drag reduction of turbulence.

Our approach holds an additional promise. Computational rheology still turns out to be a major challenge [12]. The main difficulty lies in the fact that the convective nonlinear term in the constitutive equation can lead to a local exponential growth of the components of the polymer stress tensor [13]. As a result, it is difficult to come up with robust numerical schemes. Although progress can be made by going to special variables [13], there is still a great need for an easy-to-use and robust method for simulating polymer rheology in complex geometries. As is well known, mesoscale models like DPD and Lattice Boltzmann [14] are versatile methods to simulate laminar and turbulent Newtonian flows in rheometric flows [15] and in complicated geometries [16]; we expect our extension of the DPD model to have the same advantage.

## DPD model for a visco-elastic fluid

DPD [8, 9] is an off-lattice method in which one simulates the dynamics of particles which we can intuitively think of as representing mesoscopic blobs of fluid. The interparticle interactions consist of conserved, dissipative and random forces tuned to reproduce hydrodynamic behaviour on the scale of a few particles. Particles are assumed to have mass and velocity at positions so that their equations of motion are

(0) |

For the first three terms we follow [15] and take the standard choices [8, 9]: the force is a soft repulsion,

(0) |

Here is the distance between the particles, the interparticle unit vector, and the coefficient measures the strength of this interaction. The conserved force increases linearly when particles come within the range . From here on we take , so that distances are in units of the particle radius. The dissipative force acts to equalise velocities of neighbouring particles, while the random force represents a coupling to a heat bath. We use the standard form [8, 9]

(0) | ||||

(0) |

where is a friction coefficient, is the temperature, is the relative velocity vector, is the integration timestep, and are independent Gaussian random variables with zero mean and unit variance. The form of is chosen to ensure that the fluctuation-dissipation theorem holds in the absence of applied flow [9, 17]. Note that the best convergence of the DPD model to hydrodynamic behaviour is achieved when all terms in Eq.(DPD model for a visco-elastic fluid) are roughly of the same order.

The dumbbell force is our extension of the DPD model [18]. It is motivated by the observation that the so-called Oldroyd-B constitutive equation is exact for non-interacting dumbbells with linear elastic springs [1, 2]. Linear Hookean springs, however, can lead to diverging stresses since the two partners can separate infinitely far. We therefore build our DPD model for polymer flows by using nonlinear dumbbells: we assign a unique partner to a finite fraction of the particles and introduce, in the spirit of the Finitely Extensible Nonlinear Elastic (FENE) spring constitutive equations [1, 2], a FENE-force between each particle and its partner,

(0) |

The FENE spring behaves as a Hookean spring with stiffness at small extensions and as a stiff rod when the extension is close to its maximum value .

While this force works well in smooth flows, large local stresses in turbulent flows can, because of the finite timestep, sporadically lead to a blow-up due to configurations with . It is therefore numerically advantageous to use instead in such situations a FENE-inspired nonlinear force that has no sharp maximum extension,

(0) |

To the lowest order in the extension, this nonlinear force is the same as the FENE force (DPD model for a visco-elastic fluid) if we identify with , but unlike Eq.(DPD model for a visco-elastic fluid) it remains finite for all . We will use this soft force in the last section of this Letter. Once again, we stress that the dumbbells with the force law (DPD model for a visco-elastic fluid) or (DPD model for a visco-elastic fluid) do not represent individual polymers but rather collective viscoelastic degrees of freedom.

The equations of motion are solved with a version of the velocity-Verlet algorithm [19]. Unless noted otherwise, we take , particle density and timestep . Simulations are performed in the plane Couette flow geometry with periodic boundary conditions in the (streamwise) and (spanwise) direction, and no-slip walls perpendicular to the (gradient) axis. The typical dimensions of our simulation box are so that when all particles are dumbbells we have dumbbells at density . The walls are implemented as a soft repulsion potential in the normal direction — see [18] for details.

## Rheological properties of the dumbbell DPD model

We now show that this dumbbell DPD fluid exhibits the main characteristics of polymer rheology. In a simple shear flow, the stress tensor for an incompressible polymer solution has 3 independent components: there is one independent off-diagonal component while the diagonal components are linear combinations of the so-called first and second normal stress differences and . While these both vanish for a Newtonian fluid, the normal stress differences characterise the rheological properties of a viscoelastic fluid, in particular its relaxation time . More precisely, for our plane Couette geometry we have to lowest order in ,

(0) |

where , with and being the solvent and the polymer contributions to the total viscosity. The expression for the first normal stress difference defines the effective polymer relaxation time . The second normal stress difference for most polymer solutions. For the Oldroyd-B constitutive equation, vanishes identically [1, 2].

To estimate the shear viscosity and the relaxation time for different dumbbell parameters, we have performed simulations at various shear rates and computed the components of the total stress tensor from the pairwise interactions by time-averaging the virial formula [18, 20]:

(0) |

over a time interval of length 200, which was long enough to ensure convergence; here and is the average velocity of particles at the vertical position . We have verified that only one shear stress component takes significant values (typically ), and that so that indeed , as desired. The accuracy of (Rheological properties of the dumbbell DPD model) has also been validated by comparing its results to the total momentum transferred through the walls [18].

The viscosity and the relaxation time are then obtained from linear and quadratic fits of and to as illustrated in Fig. 1(a). By varying the fraction of dumbbells from 0% to 100%, we can separate the individual contributions of and ; as shown in Fig. 1(b), the total viscosity increases essentially linearly with the concentration, so that for the choosen DPD parameters and .

At fixed polymer concentration, is approximately inversely proportional to for large , as can be inferred from Table 1. This is as expected for the relaxation of a dumbbell in a viscous medium [1]. (We have found that we can increase while keeping by using chains with more than two beads [20]). We have checked that there is little shear-thinning in our model: decreases by at most 10% at high shear rates [18]; we have also observed that the magnitude of shear-thinning increases with . We have also estimated from the decay of the shear stress autocorrelation function [11]. As shown in Fig. 1(c), it can be approximated by a single exponential that yields and for , and and , respectively. These results are within 20% of the values quoted in Table 1.

0.1 | 0.77 | 4.18 | 0.91 | 5.29 | 1.04 | 5.63 |
---|---|---|---|---|---|---|

0.125 | 0.67 | 3.07 | 0.76 | 3.55 | 0.82 | 4.22 |

0.2 | 0.49 | 2.01 | 0.55 | 2.25 | 0.58 | 2.72 |

0.5 | 0.27 | 1.2 | 0.29 | 1.23 | 0.32 | 1.23 |

1 | 0.15 | 0.64 | 0.16 | 0.66 | 0.16 | 0.65 |

The above results show that the dumbbell DPD model is a faithful mesoscopic representation of the Oldroyd-B type constitutive equation for polymer rheology.

## Polymer stretching and orientation statistics

The statistics of single-polymer stretching and orientation in shear flows has recently been studied both experimentally [3] and theoretically [21]. We now demonstrate that our DPD dumbbell model is also capable of qualitatively reproducing the single-molecule results [3, 21].

We first study the extension of the FENE dumbbells (DPD model for a visco-elastic fluid) for various values of the Weissenberg number that describes the balance between the shear flow and elasticity. Fig. 2(a) shows the Probability Distribution Function (PDF) of the dumbbell extension for three values of . Upon the increase of the Weissenberg number, the shape of changes very much in line with what has been found for single polymers in shear flows by other methods [22, 21, 4]: a peak at small and a power-law decay at large extensions for , and a rather flat distribution across a large range of extensions above the coil-stretch transition for a single polymer at . At very high the distribution becomes strongly peaked close to but we have not explored this regime.

The average orientation of a polymer in a shear flow, which in our dumbbell model translates into the average orientation of the dumbbell, is characterised by the spherical angles and . The distribution of , the angle between the dumbbell and the shear plane, is found to be well approximated by a Lorentzian, in agreement with single-dumbbell models [23, 24]. We focus here on the PDF of the angle that the projection of each dumbbell in the shear plane makes with the flow direction; in view of the symmetry of the system, we only need to consider angles .

The distribution is shown in Fig. 2(b) together with the fits to the expression

(0) |

which can be derived from the Fokker-Planck equation for [23]. In this expression plays the role of a diffusion coefficient for the orientation angle . As can be seen from the figure, the above expression fits our data very well, and this allows us to study the shear-rate dependence of ; from the fitted values of listed in the caption, one infers that decreases with increasing shear rate approximately as . Although for a FENE model one does not expect a pure scaling form, we can compare this to what is expected for a single dumbbell with a linear spring, for which , as well as to what is expected for inextensible rods, for which [22]. This comparison indicates that over the range of shear rates considered, the finite extensibility of the dumbbells already reduces the shear rate dependence of noticeably.

## Tumbling statistics

As we noted before, the orientation statistics of an isolated polymer can largely be understood by studying a dumbbell in a shear flow [21]. This also holds for the tumbling statistics [3, 25]: in a shear flow, the polymer is most of the time almost aligned with the mean flow, but due to the torque exerted on it by shear, every now and then it rapidly tumbles before coming back close to the average orientation; the tumbling time is defined as the typical time between two such events [3, 25]. In our DPD model we have many dumbbells which interact via . Nevertheless, we shall now show that our model still exhibits the same tumbling statistics as a single dumbbell in shear flow.

To compare the tumbling time with the relaxation time , we performed DPD simulations in a box and other parameters as before and followed the dynamics of dumbbells for time units. We computed their extension and azimuthal angle every . The intermittent nature of such time traces is immediately clear from the example shown in the inset of Fig. 3.

In Table 2 we list the average tumbling times for different values of and . The values of were computed from the angular dynamics as explained in the caption of Fig. 3. We have also determined the individual tumbling times from the time interval between two stretched configurations with a contracted one in between. Both methods yield comparable results, but contrary to single polymer experiments for which angular resolution is a limitation [26], analysing the angular time traces works best numerically and is least dependent on the thresholds used to define a tumbling event.

29.80 | 20.36 | 29.20 | 23.01 |

In our simulations, the distribution of tumbling times shown in Fig. 3 decays exponentially — this is known to be a robust property of a single polymer in shear flow [22, 24]. The fact that there is only one intrinsic relaxation time in our model is another important illustration that its viscoelastic behaviour is a good mesoscopic representation of the single-relaxation time in Oldroyd-B type constitutive equations. Qualitatively, the collective tumbling of DPD polymers is very much like what one expects from Brownian dynamics of a single dumbbell: the tumbling time increases with the fluid relaxation time (we almost have proportionality for ). As expected [22], we find that for not too large .

We have also performed simulations with a very few dumbbells – only 12 out of the particles were paired up with FENE springs. Except for small string stiffness (presumably because of the poor statistics and the force singularity at large stretching), calculation of the associated tumbling times gave values very close to those with dumbbells listed in Table 2. Together with the results of the previous section, this shows that the elasticity of the surrounding dumbbells does not significantly affect local dynamics. It implies that by increasing the concentration of dumbbell polymers , we can increase the total normal stress without changing the properties of the fluid on a mesoscopic scale.

## Turbulent drag reduction

To demonstrate that our approach can also be employed to study non-trivial macroscopic flows, we perform, for the first time to our knowledge, DPD simulations of high Reynolds number turbulent flows in the presence of polymers. This system is known to exhibit drag-reduction – the phenomenon that the addition of minute amounts of polymer significantly suppresses turbulence and hence the drag [7]. It has recently been successfully captured by direct numerical simulations [27, 28] and by studying the exact coherent states [29]. We have shown previously [15] that the standard DPD model reproduces the main characteristics of turbulent plane Couette flow (albeit the onset Reynolds numbers are somewhat too high perhaps due to the significant compressibility of the DPD fluid). We now study the effect of polymers on the turbulent state as a function of polymer fraction.

To do so, we now switch to the soft FENE dumbbell interaction (DPD model for a visco-elastic fluid) for the reasons we explained before. The parameters for these runs are , , , , and (the rheological properties are , and ); the simulation box has dimensions ^{1}^{1}1In these runs, the walls are modelled by a layer of continuous DPD matter (see [15]), rather than by a smooth potential..

We first drive the flow as in [15] with a force field that generates an array of streamwise vortices that are known to dominate the turbulent state at least close to the onset [30]. The driving field is turned on at time 200, then slowly ramped up to time 440, kept at a steady value till time 1160, and then slowly turned off till time 1400. By following the maximum deviation from the linear flow field, where is the spatially averaged mean streamwise velocity profile, it is then tested whether turbulence remains or whether it decays to the laminar state. The behaviour of was studied as a function of polymer fraction and the Reynolds number . We have varied the forcing amplitude over a factor 12 and found the final turbulent state (if it exists) to be independent of the forcing amplitude. Typical data for two polymer concentrations is shown in Figs. 4(ab). The saturated values of are plotted in Fig. 4(c) as a function of for various . This stability diagram illustrates the occurrence of drag reduction within our dumbbell DPD model: for increasing polymer fraction, the onset of turbulence is pushed to higher Reynolds numbers and the turbulent amplitude is reduced. A similar stability diagram was reported for the exact coherent states in plane Couette flow [29]. Finally, as Fig. 5 illustrates, not only do the dumbbells get preferentially stretched inside the streaks (blue and red regions in the velocity plot), as one might expect [28], but the fluctuations in the stretching appear to be much larger than those in the velocity.

## Conclusion

In this paper we have shown that the dumbbell DPD model has all the main characteristics of the Oldroyd-B class constitutive equations, while the statistical properties of the dumbbells are similar to those of single polymers. It can thus be used to study polymer rheology simultaneously at the macroscopic and mesoscopic scales. This makes the DPD dumbbell model an attractive candidate for simulations of polymer rheology in complex geometries and flows, in particular in situations where impact of the macroscopic flow on the mesoscopic statistical properties is of interest. We expect the model to be especially useful and easy to implement for studies of cross-slot geometries, flows past objects, extrusion, etc. Detailed information about coarse grained polymer conformation readily available in our model may even shed new light on the mechanism of drag reduction.

The authors acknowledge NWO/NCF for SARA supercomputer time (Project No. MP-06-119). The work of ES and ANM in Leiden was supported by the physics foundation FOM (ES) and the DoP program of NWO/FOM (ANM). ANM acknowledges support from the Royal Society of Edinburgh/BP Trust Personal Research Fellowship.

## References

- [1] Bird R. B., Armstrong R. C. Hassager O., Dynamics of Polymeric Liquids Vol. 1 (Wiley, New York) 1987.
- [2] Larson R. G., The Structure and Rheology of Complex Fluids (Oxford University Press, Oxford) 1999.
- [3] Smith D. E., Babcock H. P. Chu S., Science , 283 (1999) 1724.
- [4] Babcock H., Teixeira R., Hur J., Shaqfeh E. Chu S., Macromolecules , 36 (2003) 4544.
- [5] Groisman A. Steinberg V., Nature , 405 (2000) 53.
- [6] Larson R. G., Nature , 405 (2000) 27.
- [7] Virk P. S., AIChEJ , 21 (1975) 625.
- [8] Hoogerbrugge P. J. Koelman J. M. V. A., Europhys. Lett. , 19 (1992) .
- [9] Espanol P. Warren P. B., Europhys. Lett. , 30 (1995) 191.
- [10] ten Bosch B. I. M., J. of Non-Newtonian Fluid Mech. , 83 (1999) 231.
- [11] Pan G. Manke C. W., J. of Rheology , 46 (2002) 1221.
- [12] Owens R. G. Phillips T. N., Computational Rheology (Imperial College Press, London) 2002.
- [13] Hulsen M. A., Fattal R. Kupferman R., J. Non-Newtonian Fluid Mech. , 127 (2005) 27.
- [14] Succi S., The Lattice Boltzmann Equation (Clarendon Press, Oxford) 2001.
- [15] van de Meent J. W., Morozov A. N., Somfai E., Sultan E. van Saarloos W., Phys. Rev. E , 78 (2008) 015708(R).
- [16] Chen H. D., Kandasamy S., Orszag S., Shock R., Succi S. Yakhot V., Science , 301 (2003) 633.
- [17] Chen S., Phan-Thien N., Fan X.-J. Khoo B. C., J. Non-Newt. Fluid Mech. , 118 (2004) 65.
- [18] Somfai E., Morozov A. N. van Saarloos W., Physica A , 362 (2006) 93.
- [19] Groot R. D. Warren P. B., J. Chem. Phys. , 107 (1997) 4423.
- [20] Spenley N. A., Europhys. Lett. , 49 (2000) 534.
- [21] Chertkov M., Kolokolov I., Lebedev V. Turitsyn K., J. Fluid Mech. , 531 (2005) 251.
- [22] Celani A., Puliafito A. Turitsyn K., Europhys. Lett. , 70 (2005) 464.
- [23] Turitsyn K. S., JETP , 105 (2007) 655 ; nlin/0501025.
- [24] Puliafito A. Turitsyn K., Physica D , 211 (2005) 9.
- [25] Tao Y.-G., den Otter W. K. Briels W. J., Phys. Rev. Lett. , 95 (2005) 237802.
- [26] Gerashchenko S. Steinberg V., PRL , 96 (2006) 038304.
- [27] Sureshkumar R., Beris A. N. Handler R. A., Phys. of Fluids , 9 (1997) 743.
- [28] Dubief Y., White C. M., Terrapon V. E., Shaqfeh E. S. G., Moin P. Lele S. K., J. Fluid Mech. , 514 (2004) 271.
- [29] Stone P. A., Roy A., Larson R. G., Waleffe F. Graham M. D., Phys. of Fluids , 16 (2004) 3470.
- [30] Eckhardt B., Schneider T. M., Hof B. Westerweel J., Annu. Rev. Fluid. Mech. , 39 (2007) 447.