Polymers are long to very long molecules, usually built by stringing together many identical chemical repeat units (monomers). In concentrated solutions and in a melt, polymers coil around each other and hinder each other's motions: they are "entangled". Entanglements occur because bonds between two adjacent atoms in a polymer chain can never be crossed by other such bonds. The goal of my PhD study was to simulate and to understand the dynamical and rheological behaviour resulting from this entanglement effect. Because of its relative simplicity, we have chosen polyethylene (PE) as our primary system of interest.

Microscopic simulations

First, we performed molecular dynamics (MD) simulations of a melt of relatively short C120H242 chains, see the figure 1 below.

coarse-graining of C120 Fig. 1. In an atomistically detailed molecular dynamics simulation of a melt of C120H242 chains, we identify the centre-of-mass of each part of 20 carbon atoms as one "blob". In this picture one of the chains in the melt is highlighted.

Using the atomistic MD force fields, the uncrossability of bonds is automatically guaranteed because of the extremely high energy barier that has to be overcome if two bonds want to cross. Unfortunately, MD simulations are computationally very demanding. The rapid internal motions of the molecules make it necessary to calculate the interactions at femtosecond (10^-15 s) intervals. As a consequence, simulations of even the relatively small system of C120 chains took a few months to complete. From these simulations we learned that various dynamic quantities, as measured in the C120 system, can consistently be described by the Rouse model using a single set of fit parameters, provided the length scales involved are larger than the statistical segment length.

Coarse graining

The longest characteristic time scale that occurs in a system of C120 chains is about 6 nanoseconds, which is well within reach of atomistic MD simulations. However, the longest time scale increases rapidly with chain length, and already equals several microseconds for a C1000 chain. Moreover, the system size must increase as well to avoid significant overlap of a chain with itself via the periodic boundary conditions. Hence, fully atomistic simulation of the dynamics and rheology of long polymer chains is quite impossible with current-day computer power. In order to increase the time and length scales accessible in simulations, it is necessary to describe the polymers on a more coarse-grained level. If the degree of coarse-graining is larger, a particular chain is represented by less coarse particles (which we call "blobs") and the integration time step can be increased. However, as long as the coarse-grained interactions are modelled as spherical interactions, it is important that the size of a blob does not exceed the entanglement length, since otherwise no realistic entanglement effects can arise. Taking these considerations into account, we decided to represent the center of mass of 20 consecutive CH2 groups by one blob.

Thus, we derived the coarse-grained interactions between blobs from the underlying atomistic model, see figure 1 above. The resulting interactions are so soft that the uncrossability of bonds is no longer automatically met (Fig. 2).

potential of mean force Fig. 2. Potential of mean force between bonded (squares) and non-bonded (circles) blobs. Notice that kT is about 3.7 kJ/mol at the temperature of 450 K.

To prevent such unphysical bond crossings a new uncrossability constraint, the TWENTANGLEMENT algorithm, was introduced. The idea behind this constraint is to consider the bonds between consecutive blobs to be elastic bands. As soon as two of these elastic bands make contact, an "entanglement" is created at the crossing point which prevents the elastic bands from crossing (Fig. 3).

principle of TWENTANGLEMENT Fig. 3. Principle of TWENTANGLEMENT: at an earlier time, the bonds between the two depicted pairs of blobs tried to cross each other. This caused the uncrossability constraint to insert an "entanglement" at the crossing point. Since then the attractive part of the potential between bonded blobs is a function of the path length from blob i, via the entanglement at X, to blob i+1.


The advantage of coarse-graining bottom-up, from the atomistic to the mesoscopic scale, is that all time and length scales are incorporated automatically and in their right proportions. Indeed, we found very good agreement between the dynamic results of the atomistic MD simulations of C120 and the coarse-grained simulations of B6. We observed a subdiffusive exponent in the mean square displacement of the chains, a stretching of the exponential decay of the Rouse modes, and a slowing down of the relaxation of the single chain coherent dynamic structure factor.

Both the uncrossability of chains and their stiffness at smaller scales were found to be essential for these effects to occur. Interestingly, the shear relaxation modulus initially behaves Rouse-like, but after t=5 nanoseconds, the stress in the system relaxes more slowly than in a system of Rouse chains. This was attributed to a very slow interchain stress relaxation caused by the uncrossability of chains.

Next, we studied the dynamical and rheological behaviour of melts of chains ranging from C80 to C1000 (4 to 50 blobs). We found that the dynamics of chain lengths C400 - C1000 is in approximate agreement with reptation theory for large time scales, but that the approximation of a Rouse-like primitive path moving with great freedom in a tube is too strict. A better picture would be that of a primitive path that is interacting with the neighbouring chains on every length scale up to the entanglement length N_e. In fact we identified a new length scale, called the slowing down length N_s, which is smaller than the entanglement length N_e. The effective segmental friction increases rapidly around N_s leading, at constant density, to a transition in the scaling of the diffusion coefficient from D ~ N^-1 to D ~ N^-2 and a conspicious non-exponential relaxation behaviour. These effects were attributed to strong local kinetic constraints caused by both chain stiffness and interchain interactions. The onset of non-local (entanglement) effects occurs at a chain length of C120, as exemplified by deviations from Rouse behaviour of the shear relaxation modulus. Full (rheological) entanglement effects were observed only above C400, where the shear relaxation modulus displays a plateau (see Fig. 4) and the single chain coherent dynamic structure factor agrees with the reptation model. The results for the tube diameter and the plateau modulus, as well as diffusion coefficients and viscosities (Fig. 5) were found to be in good agreement with experiment.

Finally, as an application of our coarse-grained model, we studied the nonlinear flow properties of polyethylene melts by subjecting the model to a planar Couette flow. In steady state, typical effects such as shear thinning of the viscosity and a decrease of the extinction angle with shear rate were measured. Also transient effects, such as the characteristic overshoot in the shear stress and an undershoot in the transient extinction angle upon onset of shear flow, were measured and found to be in good agreement with experiments.

Shear relaxation modulus Fig. 4. Shear relaxation modulus of a melt of C80H162, C120H242, C200H402, and C800H1602, respectively. Absolute values of negative data are represented by filled circles. Dashed lines are Rouse model predictions with the measured spectrum of relaxation times. Solid lines are predictions from a mixed Rouse and reptation approach. The arrows indicate estimates of the entanglement time tau_e. Click on the picture for a larger view.

Diffusion coefficient and viscosity Fig. 5. Self-diffusion coefficient (descending curve; right scale) and viscosity (ascending curve; left scale) versus molecular weight for polyethylene melts at 450 K. Open circles are from atomistically detailed molecular dynamics simulations, closed circles are from this work, and dashed lines are fits to experimental data. Click on the picture for a larger view.

RaPiD: polymers as a single particle

We have taken the concept of coarse-graining even further and treat an entire polymer as a single particle. An important concequence is that the integrated out coordinates constitute a slowly relaxing bath, generating transient forces. We are currently applying this new method (RaPiD = Responsive Particle Dynamics) to linear polymers, star polymers, and telechelic polymers (see figure 5 below), with encouraging results. We are now also applying these methods to architecturally more complex polymers, such as H-shaped polymers and pressure sensitive adhesives (figure 6 below).

failure mode transition in transient polymer networks under shear Fig. 5. Particle-based simulations reveal a failure-mode transition in transient polymer networks under shear; from a shear banded to a fractured state. Our observations resolve apparent inconsistencies in the experimental literature on the non-linear rheology of polymer networks with reversible nodes.

Fig. 6. Cartoon of the creation of an adhesive film. Left: polymeric latex particles (white) are dispersed in water (blue) by the action of surfactants and polar groups which predominantly reside on and in the interfacial region of the latex particles (yellow). Right: when the system is dried, the water evaporates and most, but probably not all, surfactants are drained away. The polymer chains from different latex particles start to intermix, and polar groups cause a strong bonding between the particles. We developed a new simulation model that can elucidate the relation between rheology and the properties of latex particles forming such a pressure sensitive adhesive. In the model a latex particle is treated as a single, yet extensible, particle. At this level of coarse-graining transient forces are found to be essential to correctly predict qualitative features and to allow for quantitative predictions.

Continuum scale (macroscopic) simulations

polymer flow streamlines Fig. 7. Polymer flow streamlines for divergent-convergent flows at different Deborah numbers

On macroscopic scales, molecular details are often replaced by suitable constitutive relations between polymer stress and deformation history. Such so-called Computational Fluid Dynamics (CFD) simulations of viscoelastic fluids are often hampered by the severe demands on grid resolution occurring near (solid) boundaries. We have developed a combined finite-volume CFD and second order accurate immersed boundary method for viscoelastic flows that makes it possible to use relatively coarse-grids which are still sufficiently accurate. We have used the method to study the flow of viscoelastic fluids through a full three dimensional pore space (Fig. 7).

  • S. De, S. Das, E.A.J.F. Peters, J.A.M. Kuipers and J.T. Padding, A coupled finite volume immersed boundary method for simulating 3D viscoelastic flows in complex geometries, J. Non-Newtonian Fluid Mech. 232, 67 (2016).
  • S. De, E.A.J.F. Peters, J.A.M. Kuipers and J.T. Padding, Viscoelastic flow simulations in model porous media, submitted to Phys. Rev. Fluids (2017).
  • S. De, J. van der Schaaf, N.G. Deen, J.A.M. Kuipers, E.A.J.F. Peters and J.T. Padding, Elastic instabilities in flows through pillared micro channels, submitted to Microfluidics and Nanofluidics (2017).
  • S. De, E.A.J.F. Peters, J.A.M. Kuipers and J.T. Padding, Viscoelastic flow simulations in random porous media, submitted to J. Non-Newtonian Fluid Mech. (2017).


We have collaborated with Prof. Juan Colmenero from the University of Basque Country in San Sebastian, Spain. In this project, a PhD student (Roberto Perez Aparicio) was applying the TWENTANGLEMENT method to the case of polypropylene. We are collaborating with Prof. Dimitris Vlassopoulos at the FORTH institute on Crete (Greece) in applying the RaPiD model to solutions of star polymers. We are also collaborating with Dr Joris Sprakel (previously in the group of Prof. Martien Cohen-Stuart at Wageningen University and now in the group of Prof. Dave Weitz at Harvard University) in applying the RaPiD model to micelle-forming telechelic polymers. I have worked with Prof. Christian Bailly, Prof. Wim Briels, Dr Lalaso Mohite and Dr Dietmar Auhl at the Université catholique de Louvain on the rheology of pressure sensitive adhesives. Within this project we also collaborated with the experimental groups of Prof. Creton at ESPCI Paris and Dr Thomas Schweizer in the group of Prof. Öttinger at ETH Zuerich. The continuum scale simulation work is supported by a FOM-Shell Computational Sciences in Energy Research grant.