Open boundary molecular dynamics of sheared star-polymer melts

Abstract

Open boundary molecular dynamics (OBMD) simulations of a sheared star polymer melt under isothermal conditions are performed to study the rheology and molecular structure of the melt under a fixed normal load. Comparison is made with the standard molecular dynamics (MD) in periodic (closed) boxes at a fixed shear rate (using the SLLOD dynamics). The OBMD system exchanges mass and momentum with adjacent reservoirs (buffers) where the external pressure tensor is imposed. Insertion of molecules in the buffers is made feasible by implementing there a low resolution model (blob-molecules with soft effective interactions) and then using the adaptive resolution scheme (AdResS) to connect with the bulk MD. Straining with increasing shear stress induces melt expansion and a significantly different redistribution of pressure compared with the closed case. In the open sample, the shear viscosity is also a bit lowered but more stable against the viscous heating. At a given Weissenberg number, molecular deformations and material properties (recoverable shear strain and normal stress ratio) are found to be similar in both setups. We also study the modelling effect of normal and tangential friction between monomers implemented in a dissipative particle dynamics (DPD) thermostat. Interestingly, the tangential friction substantially enhances the elastic response of the melt due to a reduction of the kinetic stress viscous contribution.

---

1. Introduction

Die swelling¹ is a well known phenomenon in polymer melts and most viscoelastic liquids which consist of the sudden expansion of the liquids after exiting a slit or orifice. The most frequent explanation has a microscopic origin: molecules elongate in the stream direction and compress perpendicularly exerting extra elastic pressure in the normal planes. This leads to the so-called normal stress differences which are the landmark of viscoelasticity. Despite this accepted view, the devil is in the details and although considerable effort has been carried out since the middle of the last century (see e.g.ref. 1 for historic details), accurate modelling of polymer melts is a very difficult task. Any of the many constitutive relationships¹ for continuum models cannot generally predict the rheology of a new molecular polymer design. In turn, the huge span in time scales in any standard polymer melt limits the scope of molecular dynamics (MD) to simple rheological tests with extremely small samples under simple (usually steady) flows. However, a detailed account of bonded and non-bonded interactions in atomistic (AT) simulations (see e.g.ref. 2–4) is able to grasp relevant information, maybe then to feed continuum (fluid dynamics) models. In between, coarse grained (CG) molecular modelling is useful due to many reasons. We give at least a couple of reasons: first, polymer science has some degree of universality which benefits the use of simplified models, quite often able to provide insight and valuable predictions.^4,5 Second, the theory of coarse graining to extract precise coarse potential interactions pertaining to the atomistic model at hand is now advancing at a relatively fast rate.^4,6,7 More recently, the community has started to recognize the relevance of the dynamic aspects of coarse graining either based on GENERIC⁸ or the Mori–Zwanzig formalism.⁹ The idea is to perform short atomistic simulations⁸ to extract the quantities determining the reversible and irreversible dynamics of the slow variables, such as the friction kernels.⁹ These friction kernels (between polymer “blobs”) are naturally implemented in the dissipative particle dynamics (DPD) method as, notably, Español et al.⁹ showed (under the Markovian and pairwise approximations) that DPD can be formally derived from the Mori–Zwanzig coarse graining route. Another relatively newer route is the use of hybrid models concurrently combining continuum and molecular simulations. Yasuda et al.¹⁰ and other (mostly Japanese) groups have exploited a version of these hybrids to model polymer melts. We have also worked in this field using essentially the same technique used in this work for simple fluids.^11,12 Contrary to Yasuda et al.¹⁰ and other hybrid schemes, the present method is designed to open up the simulation box so as to consider melt expansions. Here, we apply the idea to study the rheology of much larger, polymeric molecules. At present, we restrict to the open boundary molecular dynamics (OBMD) simulation,¹³ without connection to the continuum side. However, the present method naturally connects with continuum fluid dynamics via hybrid schemes.^14–20

Another aspect of this work is the use of DPD as a tool to check the effect of monomer friction on the rheological behaviour of these CG models. In particular, the tangential friction between blobs naturally arises while performing dynamic coarse graining⁹ but is seldom included in these sort of analyses (see for example the DPD study for melts by Fedosov et al.²¹). To conclude our comments on methodological aspects, we also note that the present simulations in open domains use a useful trick which consists of using an even lower detailed molecular model to feed the molecular reservoir close to the open system boundaries. The idea is taken from the so-called “adaptive resolution scheme” (AdResS)^22–34 and permits to generalize the use of AdResS in standard periodic (closed) boxes (see e.g.ref. 11, 12 and 35–37).

The main physical question we pose here is what are the rheological consequences of imposing a fixed pressure load to a sheared sample of polymer melt, compared to the case of shearing at a fixed volume. This question was also raised by a bunch of groups spread over the last two decades^38–43 with contradictory results. It is indeed a particularly relevant question for molecular simulations because the vast majority of numerical studies on melts have considered closed (usually periodic) systems, while many rheological experiments are carried out under a normal load: melts across slabs in ref. 44 (cited experiments therein) or a cone and plate rheometry under a fixed load.¹

From the fundamental side, this question connects with the already mentioned die swelling phenomenon whose details are still not completely understood. It also connects with another interesting question which is, what consequences do boundary constraints have on flowing (far-from-equilibrium) polymer melts (see ref. 45 for a recent study). The OBMD is a flexible tool for these questions on boundary constraints because it can be tuned to permit different ensembles⁴⁶ such as the grand canonical, isoenthalpic, isothermal, constant stress (Neumann-like) or constant shear (Dirichlet-like).¹⁴ In particular, it could be useful to validate theories for non-equilibrium thermodynamics (such as extended thermodynamics^47,48), or the far from trivial fate of fluctuations of mass (related to sound modes) and momentum in sheared complex fluids, which often lead to undesirable instabilities in sheared or extruded melts, like the shark skin.^1,49

Even in unentangled melts, the influence of boundary or global constraints on the density expansion of sheared melts is still poorly understood with the studies presenting contradicting results. At least there is a consensus on the fact that for a given shear, the shear viscosity η is larger in the isochoric (NVT) constraint than under either a constant pressure^39,40 or a constant load.^43,44 A clear manifestation of this effect was presented in a numerical study of Thompson et al.,⁴⁴ in which a slab of lubricating liquid (20-mers) flowing between two solid walls presented a shear thinning exponent β_η ≃ 2/3 under a constant normal load, while just 0.5 under a constant volume (here η ∼ [small gamma, Greek, dot above]^−β_η with [small gamma, Greek, dot above] the shear rate). This effect has been ascribed to the shear dilatancy manifesting in a larger hydrostatic pressure under a constant volume.^39,40 Indeed, the viscosity of polymer melts often increases with the pressure⁴⁹ as also observed in ref. 44. In Section 6, we offer a more precise analysis showing that the spring stress is reduced in the open system due to the smaller intermolecular friction in the expanded melt.

There are several kinds of “isobaric” conditions: a constant hydrostatic pressure (P_iso) is usually termed isobaric, while a constant normal load (here P₂₂) is closer to industrial processes like slit extrusion¹ or lubricants.⁴⁴ Both can show substantially different rheological behaviours when compared with constant volume studies (see e.g.ref. 44). Experimental studies of the Couette flow under a constant load are also scarce and indicate shear expansion and a measurable increase of the melt viscosity with increasing external pressure.^50,51

The most striking differences in the computational literature are found in the density variation with shear. Dlugogorski³⁸ (for a FENE dumbbell) and Daivis and Evans (modelling decane)³⁹ seem to be the first to perform a molecular simulation showing the density decrease under shear (using isobaric conditions). Daivis and Evans use the term “shear dilatancy”, following the term used by Reynolds on the same phenomenon (see the quotation in ref. 39). Note that “shear dilatancy” has been later often used instead of “shear thinning”, but here it is not. By contrast, Xu et al.⁴⁰ (attractive linear chains up to 50-mers) reported just the opposite result (compression under shear). For (purely repulsive) branched chains (under NVT), they also reported a reduction in the hydrostatic pressure with shear. A subsequent study by Matin⁴³ for linear chains at a constant load (and chain lengths up to 50 monomers) found shear dilatancy and also a non-monotonous trend for the hydrostatic pressure (as Xu et al.⁴⁰ and others^2,3 did). Shear dilatancy was later also found by Bosko et al.⁴¹ while analysing dendrimer melts under isobaric conditions. Consistently, they found a pressure increase under a constant volume. Shear dilatancy is the trend we also observe but, under our imposed constant normal load (which is closer to the experimental setup^1,39,44) we find that the density (and viscosity) is controlled by the load and not by the hydrostatic pressure. However, there is a lack of studies and the question remains about universality of shear dilatancy of polymer melts under a Couette flow (consistent with the die swelling phenomenon under a Poisson flow¹).

To conclude the introduction some words should be said about the technological relevance of star molecules (see discussion in ref. 52–54) which certainly arises from their unique dynamic features. Star molecules present a broad range of relaxation times associated with different molecules pulsations (rotation, elastic deformation, and arm entanglements) analysed in ref. 53. Each relaxation time triggers a change in the rheological regime once the external perturbation (shear) exceeds the corresponding threshold rate. The present results also illustrate this phenomenon. Moreover, star molecules bridge the gap between linear polymers and colloids^54,55 and can present interesting (colloidal-like) ordering effects, sometime enhanced due to their softer character.^53,56–59 In this context, a suggesting observation in these simulations is the onset of some ordering in the neutral direction at large shear rates whose origin (hydrodynamic or entropic) remains to be established.

In Section 2, we briefly present the OBMD method which was otherwise more thoroughly explained in ref. 13. The star molecule melt model presented in Sections 3 and 4 shows that the OBMD correctly reproduces thermodynamic equilibrium according to the grand canonical ensemble. In Section 5, we present the results for sheared melts in the absence of tangential friction between monomers and analyze the results in Section 6 according to the pressure balance. This serves to enlighten the discussion on the effect of tangential friction in Section 8. Finally, Section 9 discusses some interesting results obtained for melts presenting severe viscous heating that depends on the characteristics of the DPD thermostats and their friction kernels. Comparison with the previous results is made in Section 10, while conclusions and future outlook are given in Section 11.

2 Open boundary molecular dynamics

We begin by briefly explaining the OBMD method which combines features of the open MD⁴⁶ and adaptive resolution.¹¹ The reader is referred to a review on open MD¹⁴ and ref. 13 for a more detailed presentation of the present OBMD implementation to star polymers. The OBMD simulation is carried out in an open rectangular box which, in the present setup, permits the MD domain to exchange mass and momentum through two of its boundaries (along the x₂ direction) with a reservoir (called buffer) which is maintained at some desired thermo-mechanic state. The two domains of finite extent of the buffer embed the central part of the box (the MD domain). They allow for molecular insertion or deletion so as to keep their average molecular density fixed (typically to a fraction between 0.5 and 0.7 of the bulk density). The OBMD is therefore not periodic in the coupling direction. Molecules are free to enter or leave the buffer from or to the MD domain, but in doing so they cross another layer where they gradually change their atomistic resolution (73 monomers for the star molecule considered hereby) to a reduced CG model, comprising one only spherical “blob” per molecule. Obviously the CG layer is placed inside the buffer domain (which here also contains a smaller atomistic part). This strategy permits to perform an otherwise impossible task: the insertion of new polymeric molecules into the melt. New molecules are inserted into the low resolution layer of the buffer, where soft CG interactions govern the dynamics of the blob-model polymers. Soft effective interactions can be obtained from the Boltzmann iteration procedure,^60,61 although we shall see that in principle, the consistency of OBMD (in terms of pressure balance across the layers) does not depend on the CG potential chosen. The insertion of these blob molecules is carried out by the USHER scheme⁶² and the change from CG to monomer molecular resolution (usually termed atomistic resolution, AT) is carried out by the AdResS.^22–34

The dynamics of the monomers can be described by the following equations of motion:

Here r_i denotes the position of the i-th particle, v_i its velocity, and m_i its mass. The total force acting on this particle has three contributions: the external force F^ext_i acting only on the particles at the buffer (to impose the desired momentum flux); the adaptive resolution force F^ad_i, which accounts for all types of particle–particle interactions, and the thermostat F^th_i contribution (here, it is applied to the whole system).

The adaptive resolution force F^ad_i is constructed to allow for a momentum conserving “alchemic” transformation of the molecules, which takes place gradually along the transition layer (where 0 < w < 1, see below). The transition is achieved by the following linear combination of the AT and CG forces (F^AT_αβ and F^CG_αβ, respectively),²²

F^ad_αβ = w(x_α)w(x_β)F^AT_αβ + (1 − w(x_α)w(x_β))F^CG_αβ. (3)

Both expressions are correspondingly weighted by a position dependent function w(X), whose value equals 0 in the CG region and 1 in the AT one and gradually changes in between (transition layer). The adaptive resolution force provided by the AdResS is not derived from a Hamiltonian and does not conserve the energy.²³ It is however constructed to obey Newton's third law, which ensures the conservation of total linear momentum of the system and can thus be used to study fluid flows. This fact is unimportant for the present study which targets isothermal sheared systems (not-conserving energy anyhow). The OBMD model might be generalized by using the recent Hamiltonian AdResS (H-AdResS)^63,64 which is based on a global Hamiltonian (i.e. also allows Monte Carlo simulations). In such a hypothetical case, extra care should be taken with the momentum conservation because of the presence of drift-forces in H-AdResS coming from the free energy difference between the AT and CG models.⁶⁵

An essential function of the buffer region is the imposition of boundary conditions to the open MD box. This is done by adding an extra “external” force at the buffer regions, F^ext, calculated from eqn (4) (see e.g.ref. 66).

Here, P_out and P_in represent the total linear momenta of the particles that were removed and inserted into the simulation in the last time step of integration δt. is the momentum flux tensor, while n is the outward normal vector of an open-end plane of the box.¹¹ In general, the pressure tensor contains normal and tangential contributions, i.e.. The external force is designed to exactly conserve the linear momentum over the whole particle system (buffer + MD) and it is distributed among the buffer particles according to F^ext_i = G(x_i)F^ext. To allow a different distribution of the normal and tangential forces, the distribution function G is chosen to be a tensor defined by eqn (5),¹³with g_∥ determining the spatial distribution of the normal force and g_⊥ the distribution of shear stress. Both functions are depicted in Fig. 2.

Many OBMD applications (as those illustrated hereby) involve transfer of momentum (the pressure tensor) from outside the MD domain. This requires a momentum conserving thermostat. In production runs, we used the DPD thermostats (while a strong damping Langevin for equilibration purposes). Our choice for the DPD thermostat is not only to conserve momentum (in principle one could use the Lowe-Andersen⁶⁷) but also because of modelling purposes. The Mori–Zwanzig formalism, under Markovian conditions leads to coarse-graining dynamics with DPD-like equations of motion.⁹ The solid theoretical results show that the friction kernels introduce an important modelling aspect. The friction kernels in a CG model of some real melt, should ideally be measured from force–force correlations of the detailed all-atom model.⁹ Here, we adopt a simpler but yet useful route, which is to study how friction affects the rheology of the model melt. The generic form of the DPD thermostat used is

where R̃_ij is the fluctuating force constructed to satisfy the fluctuation-dissipation under equilibrium conditions. We refer the reader to ref. 68–70 for details on the DPD implementation. As in ref. 9, 69 and 70 the friction kernel has normal and tangential components,

Γ = γ_∥n_ijn_ij + γ_⊥t_ijt_ij,

where n_ij is the vector joining two monomers i and j and t_ij determines the directions in the perpendicular plane. The kernels γ_∥ and γ_⊥ are distance dependent. Its shape (for a particular all-atom model) can be obtained from dynamic CG,⁹ here we will use Heaviside functions with a certain cutoff distance (see Table 1 for details). Most thermostats require or introduce some form of friction, albeit, the great majority of simulations of polymer melts do not consider thermostating (and its added friction) as part of the molecular model, but just a way to remove the heat dissipated under shear. Also most DPD simulations, such as the relatively recent work on sheared melts²¹ do not introduce the tangential friction between blobs, but rather take the form of the most standard DPD kernel (normal friction alone). However, as also pointed out by Padding and Briels,⁵ friction should be considered as a part of the CG model. Indeed, as shown by Hijon⁹et al. (CG of a star molecule as the unit blob), tangential and normal frictions can be quite different from each other. Here, we start to explore how tangential and normal frictions affect the rheology of a star molecule under far-from-equilibrium conditions.

Table 1 Thermostats used in simulations. Standard means a standard DPD thermostat and the adaptive DPD is explained in eqn (7). The transverse DPD thermostat from ref. 70 is denoted as “tdpd”

Label	Kernel cutoff (R^DPDcut)	γ ∥	γ ⊥
Sdpdshort	2^1/6σ	[1.0–20.0]	0
Sdpdlong	1.5 × 2^1/6σ	1.0, 5.0	0
Adpd	2^1/6σ	1.0, 5.0	0
Tdpd	1.5 × 2^1/6σ	1.0	1.0

---

In the following sections, we present and analyse the results obtained for the DPD thermostat with no tangential friction between monomers (blobs). This analysis is then used to understand the effect of the tangential friction, considered in Section 8. In the case of the normal friction alone, the substantial heat dissipated by the sheared melt requires from us to implement a slight modification of the DPD thermostat to maintain a constant (kinetic) temperature at the largest shear rates. Details of this modified DPD thermostat (we call it “adaptive”) and friction kernels are given in Section 3.

Importantly, simulations in closed isothermal (periodic) boxes (NTV) are used as a reference to investigate the effect of open boundaries. Shear flow in the closed system has been simulated using Lees–Edwards boundary conditions⁷¹ and the SLLOD algorithm.^72,73

3 Setup and melt models

The simulation setup is illustrated in Fig. 1 and 2. The polymer melt is exposed to a Couette flow in the x₁ direction being sheared along the x₂ direction (gradient direction). The vorticity or neutral direction is x₃. In the closed box, we use the SLLOD^72,73 dynamics to impose the desired shear rate [small gamma, Greek, dot above] in a closed periodic box (with a constant particle number and volume). In the open setup, the surfaces located at x₂ = ±L₂/2 are subjected to equal normal pressures p^ext₂₂ and (opposite sign) tangential stresses ±p^ext₁₂, in such a way that the rotational part of the shear flow turns counter-clockwise in the flow–gradient (x₁ − x₂) plane. No constraint is imposed to the remaining component of the pressure tensor, resulting in 〈p₁₃〉 = 0 and a self determined 〈p₃₃〉. The box is periodic in the other two directions x₁ and x₃ so this setup corresponds to an slice of polymer melt with a fixed load at two of its boundaries (at x₂ = ±L₂/2).

Fig. 1 Schematic representation of the open system at equilibrium along the longitudinal direction. All monomers and polymers in the system are of the same kind. Different colors are used for the sake of clarity of the picture. In the region of interest polymers are represented in the highest, i.e. atomistic (AT) resolution. Buffer regions (see text) are heterogeneous, i.e. containing regions of different resolutions. The change of resolution from AT (dots) to coarse-grained, i.e. CG, (spheres) occurs in the hybrid region (HY) of each buffer, carried out by the AdResS. New molecules are inserted into the CG part of buffers. The system is open at both ends of the box. The upper part of the figure depicts schematic representations of a molecule inside regions of different resolutions.

Fig. 2 Buffer distribution function. The force F^ext, that acts on buffer regions in order to impose the boundary conditions at open ends of the box to the region of interest, is distributed among molecules inside buffers via the depicted functions, forming the distribution function tensor G, given by eqn (5). The force acting on each molecule equals F^ext_i = G(x_i)F^ext.

The melt is made of the star polymer model already presented in ref. 9 and 74. Each polymer consists of 73 monomers, i.e. 12 arms of 6 monomers attached to the central monomer. In what follows, we use m₀, σ₀ and ε₀ for mass, length and energy units and we will arbitrarily set these units to m₀ = 1, σ₀ = 1, and ε₀ = 1, respectively. The resulting time unit is τ₀ = σ₀(m₀/ε₀)^1/2 = 1. Excluded volume interactions of monomers are modelled by the repulsive Weeks–Chandler–Anderson interaction with diameter parameter σ = 2.415σ₀ and energy parameter ε = 1. The interactions between two adjacent bonded monomers are harmonic with a linear spring of stiffness constant K = 20.0ε/σ₀². The equilibrium distance between non-central monomers is r^eq_ij = 2.77σ₀, while the equilibrium distance between the central monomer and the first monomer of an arm is r^eq_ij = 3.9σ₀. The size of the simulation box is 390 × 117 × 117 (in units of σ₀).

Simulations were performed at a fixed constant (monomer's kinetic) temperature (T = 4.00 ± 0.01). The results presented in the following sections correspond to simulations obtained using a modified DPD thermostat which we label as “adpd”. This adpd thermostat has no tangential friction γ_⊥ = 0 while γ_∥ = 1 (some results also with γ_∥ = 5). The cutoff of the Heaviside friction kernel is R_dpd = 2^1/6σ. We refer to Table 1 for thermostat and kernel details. We recall that the effect of the tangential friction is analysed in Section 8. In Section 9, we illustrate the heat dissipation and temperature increase observed upon using standard thermostats DPD with normal friction.

Thermostating sheared polymer melts is a delicate issue due to the large amount of heat they dissipate. At a large enough shear rate or large shear stress, the temperature of the melts increases. The same phenomenon is also observed in experiments and industrial processes (extrusion) at high shear rates (typically above 500 Hz).^47,75,76 Phenomenological temperature “corrections” for the melt viscosity are often used in industry and experiments.⁷⁷ Although this problem goes beyond the present manuscript, the temperature of a system under a non-equilibrium steady state is also a fundamental problem because equipartition is lost and the different temperature definitions present slight variations (kinetic versus configurational temperature⁷⁸). In the literature, few studies report on the problem of viscous heating in molecular simulations (see the exceptions in ref. 79–81) and many published materials elude reporting on possible temperature variation in their sheared thermostated systems. After the present experience, we believe that some of the data presented in previous papers might be somehow biased by temperature. We will show an indication later in Section 9.

We considered four different thermostats: the “standard” (sdpd)⁶⁹ and “transverse” (tdpd)⁷⁰ DPD thermostats and also a modified DPD thermostat. The latter is able to extract larger amounts of heat by self-adapting its temperature parameter T_DPD which controls its random force term. This “adaptive DPD thermostat” (adpd) as we call it, dynamically adjusts T_DPD according to a sort of coupled heat equation,

where T_tgt = 4 is the temperature of the target system and T_MD is the kinetic temperature obtained from the variance of the peculiar velocities of the monomer r_i, u_i = v_i − v_f(r_i). Here, v_f(r_i) is the flow velocity at the position of the monomer r_i evaluated on-the-fly from (time averaged) the binned x₂ coordinate. Eqn (7) resembles the characteristic equation of the Berendsen thermostat,⁸² where the linear differential equation in time is solved for the current temperature of the system. In our case, on the other hand, T_DPD is just the temperature incorporated in the equations for the DPD thermostat and not the actual temperature of the system. The basic idea is simple: if the system under non-equilibrium sheared state is producing substantial heat due to the friction, the noise term (or in physical terms, the hypothetical reservoir temperature T_DPD) of the thermostat should be made colder to rapidly extract heat. The adpd thermostat does not alter the equilibrium state (indeed we first checked this fact), but reduces the noise term under substantial shear.

The thermostat time τ_DPD (acting like a coefficient of heat transfer to the “reservoir”) was set to τ_DPD = 100δt. The thermostat nominal temperature T_DPD was updated using a simple explicit Euler scheme for eqn (7) with a time step of 100δt. All the thermostats used are applied to the relative velocities of monomers. The friction kernels of the DPD thermostats (damping and noise terms are constructed using the same kernel⁷⁰) are chosen to be Heaviside functions with a cutoff distance R^DPD_cut, i.e. γ(R) = 1 for R ≤ R^DPD_cut and zero otherwise.

All the results presented here correspond to a normal friction kernel γ_∥ = 1 (some results for the adpd thermostat were also carried out for γ_∥ = 5 to test sensitivity). To test the effect of tangential friction we also run simulations with γ_∥ = γ_⊥ = 1 in the tdpd thermostat. The thermostat details and labels used are given in Table 1.

The integration step ranges from δt = 0.01τ₀ to δt = 0.005τ₀ for the highest shear rates. Note that τ₀ = 1 is smaller than the standard Lennard-Jones time (monomer–monomer interaction), τ = σ(m/ε)^1/2 = 2.415τ₀ where m = 1 is the monomer mass.

The equilibration of the lump of melt in the OBMD simulation is conducted by a modified version of the AdResS, whereby the weighting function w in eqn (3) is gradually increased in time starting from w = 0 (CG model).⁸³ The weighting function is therefore switched from a position-dependent to a time-dependent one. The resolution is thus gradually sharpened from CG to AT. The procedure is described in detail in ref. 13. After equilibration, each simulation is run for 10 000τ₀.

4 Melt at equilibrium

4.1 Characteristic times

In view of the close relationship between the structure and dynamics that takes place in sheared melts (and complex fluids in general), it is interesting to present the range of physical times of the melts before analysing their structural transformation with shear. In star polymers, one can observe three types of relaxation phenomena.⁸⁴ First is elastic deformation of the overall shape of polymers, second relaxation occurs via the rotational diffusion, and the third one regards disentanglement of arms of every star polymer. Each of the relaxation processes can be estimated from the integral of the corresponding normalized autocorrelation function (ACF), via, and these are given by eqn (8)–(10), respectively.R_i represents the center-end vector of arm i, R_i its length, t time, and f the number of arms of each polymer. And i and j are indices of different arms within the same polymer. Each autocorrelation function decays with its characteristic time of the relaxation process^84,85 and these are given in Table 2 for the different DPD thermostat friction kernels considered. We checked that under an equilibrium state all thermostats produce consistent results, in terms of pressure and density, while the correlations and the characteristic times differ, as they should. The adpd and the sdpdshort correspond to the same friction kernel and produce similar relaxation times at equilibrium (within statistical uncertainty) indicating that the adpd modification does not essentially alter the dynamics. We observe that the disentanglement of the arms occurs more rapidly than rotational diffusion and that is due to the short length of the arms, as each contains only 6 monomers. The longest time is the diffusion time for the center of mass (CoM) of molecules τ_dif = R_g²/D, where R_g = 7.6σ₀ is the radius of gyration of polymers and D the diffusion constant of CoM, which are different for every friction kernel. The differences in the relaxation times between the sdpdshort and sdpdlong thermostats are due to the larger thermostat cutoff distance R^DPD_cut leading to a larger interparticle friction in the sdpdlong case.^70,86 Interestingly, the ratio between relaxation times and viscosity is similar for all the thermostats in Table 1 (values coincide within error bars), regardless of the kernel and R^DPD_cut.⁸⁷ This indicates that the translational and orientational dynamics are affected in a similar way by the thermostats.

Table 2 Polymer time scales obtained at equilibrium for different setups and thermostats (see text). All times are in units of τ₀ = τ/2.415, where τ = σ(m/ε)^1/2 is the standard Lennard-Jones time scale for monomers. In all standard and adaptive DPD cases, the thermostat damping constant is γ_DPD = 1m₀/τ₀, while for the tdpd thermostat γ = 1m₀/τ₀ and γ_⊥ = 1m₀/τ₀. τ_elas, τ_rot, and τ_arm are defined as characteristic decay times of the ACF given by eqn (8)–(10), respectively. Diffusion time is defined as τ_dif =R_g²/D, with R_g = 7.6σ₀ the average radius of gyration in equilibrium

Simulation	τ arm	τ elas	τ dif	τ rot
Open sdpdshort	33 ± 2	6 ± 1	800 ± 100	55 ± 5
Closed sdpdshort	33 ± 2	6 ± 1	800 ± 100	55 ± 5
Open sdpdlong	55 ± 2	10 ± 1	1200 ± 100	100 ± 5
Closed sdpdlong	55 ± 2	10 ± 1	1100 ± 100	100 ± 5
Open adpd	33 ± 2	3 ± 1	700 ± 100	59 ± 5
Closed adpd	33 ± 2	3 ± 1	800 ± 100	59 ± 5
Open tdpd	195 ± 5	35 ± 5	4800 ± 400	332 ± 5

---

Using the adpd thermostat (normal friction alone) we get τ_elas = 3 ± 1, τ_arm = 33 ± 2, τ_rot = 59 ± 5 and τ_dif = 700 ± 100, which illustrates the wide range in time scales involved in these sort of simulations. These times are similar in the open and closed systems (in equilibrium) and are compared with those obtained with other thermostats and kernels in Table 2. However, as shown below the rheology of the melt is not determined by the molecular diffusion, but rather by the molecular relaxation times. For the star molecule under study with f = 12 arms of length l_a = 6σ (i.e. 6 monomers per arm), molecular rotation is slower than the elastic relaxation of arms, the ratio being τ_rot/τ_arm > 10 in all the cases considered (see Table 2). This ratio determines the type of rheological behaviour of the melt according to a theoretical approach based on solving the Fokker–Planck equation for the bond distribution.^88,89 We will come back to this issue in Section 5.3.

4.2 Equation of state

As explained in ref. 13, the OBMD simulations of the melt at equilibrium provide the correct average thermodynamic variables. The pressure equation of state obtained in OBMD agrees with that obtained from a standard NVT simulation for all the volume fractions studied.¹³ More precisely, the average equilibrium density 〈ρ〉 (p^ext₂₂) obtained with OBMD at a fixed normal pressure p^ext₂₂ is consistent with the equilibrium pressure p calculated in a closed MD simulation at a fixed density ρ = 〈ρ〉 (p^ext₂₂). In passing, it should be highlighted that the same CG potential was used for all the melt densities considered, indicating that the pressure consistency at the MD domain is independent of the coarse-grained potential used in the buffer.

In the following, we work with the polymer volume fraction, defined as Φ = ρ_Nπσ³/6, where ρ_N = N/V is the number density (N and V respectively represent the number of monomers and the volume of the region of interest and the mass density is ρ = mρ_N with monomer mass m = 1, thus ρ = 0.136Φ). For the NVT ensemble, the study conducted here corresponds to p = (0.093 ± 0.001) and Φ = 0.20, fixed for any shear rate. In the open box, we fix p^ext₂₂ = 0.093 and for zero shear rate (equilibrium) we get 〈Φ〉 = (0.20 ± 0.01). Fig. 3 shows the normalised density profile of polymers (ρ^slab_M/ρ^avg_M). Here ρ^slab_M denotes mass density of polymers in each slab, where it is measured, and ρ^avg_M its average value. The latter is constant in closed simulations, i.e. ρ^avg_M = 0.0271, while its calculated value in open cases equals ρ^avg_M = 0.0271 ± 0.0001. We observe that the obtained density profile is flat in the region of interest with some minor artefacts, which are due to the lack of statistics. Along the buffer zones the density gradually decreases as a consequence of the application of the external pressure p^ext₂₂. As it has been explained in previous related studies on hybrid atomistic-continuum schemes,^14,66 this inhomogeneity does not affect the transfer of pressure and stress from the exterior, provided that the density profile is flat around the hybrid interface. This is indeed the case, as it can be seen in Fig. 3 (see the interface between the “region of interest” and “buffer”). The radial distribution function (RDF) of CoMs of molecules is in perfect agreement with NVT simulations.

Fig. 3 Normalised density profile (NDP – top panel) and comparison of RDFs (bottom). The NDP is depicted in the direction, in which the box is open. It is normalised to the desired value of the density, which corresponds to the occupational factor Φ = 0.2. RDFs are calculated from molecules inside the region of interest in open simulation (red) and closed systems (green).

4.3 Mass fluctuations

In the grand canonical ensemble, mass fluctuations are related to the integral of the RDF,⁹⁰ so the excellent agreement between the RDFs obtained from NVT and open boxes suggest that the fluctuations in the number of molecules inside the domain of interest should be thermodynamically consistent with the grand canonical prescription: namely, the relative fluctuation in polymer mass M = ρ_MV should be Std[M]/M = [k_BT/(Mc_T²)]^1/2, where c_T is the isothermal sound velocity c_T² = (∂p/∂ρ)_T. c_T, which is evaluated using the pressure equation of state from ref. 13, is shown in Fig. 4. The agreement between OBMD and the NVT ensemble is excellent. In the same figure, we plot the results for the relative mass fluctuations at different pressures (average densities) and compare with the grand canonical prediction. We find that the agreement is very good, notably because of the tiny relative mass fluctuations in the (not small V = 3 504 384σ₀³ = 248 805σ³) volume considered for our open box: which range between 3.6% and 0.26% at the largest density considered (Φ = 0.2). In terms of mass density variance σ_ρ² = ρk_BT/(Vc_T²), for the state we considered hereafter under shear (p = 0.093ε₀/σ₀³, Φ = 0.2 and ρ = 0.0271m₀/σ₀³) the OBMD result is σ_ρ² = (5.0 ± 0.3) × 10⁻⁹m₀/σ₀³ in excellent agreement with the variance predicted for the μVT ensemble 4.79 × 10⁻⁹. The conclusion of this study strongly supports our claim that the OBMD equilibrium simulation samples the grand canonical ensemble without any (or negligible) bias. It has to be said that the value of the external chemical potential μ^ext = μ(p^ext₂₂,T) however cannot be imposed in OBMD, although it could be reconstructed following the standard Gibbs–Duhem route with varying external pressure. New implementations of the AdResS^29,65 might also be used to evaluate μ.

Fig. 4 Top panel: The relative fluctuation of the polymer mass (standard deviation over average polymer mass M) in the interest MD-domain of the open setup (see Fig. 2). Comparison is made with the grand canonical theoretical result (green dashed line) Std[M]/M = [k_BT/(Mc_T²)]^1/2 and, also included, the ideal gas limit (red dashed) obtained with the isothermal sound velocity c^id_T = (k_BT/M_m)^1/2 with M_m = 73m₀ the molecular mass. The volume of the interest domain is V = 156 × 117 × 117σ₀³. Recall that the monomer LJ diameter is σ = 2.415σ₀ and its mass is m₀ = 1. The right ordinate axis shows the values of the isothermal sound velocity c_T = (∂p/∂ρ)_T (from the pressure equation of state) that is compared with ω/k_eff obtained from the oscillation frequency ω of the total mass autocorrelation function in the MD domain (see bottom panel). The effective wavenumber is k_eff = (π/370)σ₀⁻¹ and the total open box is L₂ = 390σ₀. Bottom panel: the time autocorrelation function of the mass in the MD domain at equilibrium with imposed external pressure p^ext₂₂ = 0.001ε₀/σ₀³ (dashed line is the fit to extract ω, see text); the inset illustrates the evolution of the total mass in the MD domain for this case.

A typical outcome for the time evolution of the total mass of polymers in the MD domain is shown in Fig. 4. It presents oscillations, suggesting that it might contain some information about the sound velocity of the system.

Sound propagation can indeed be studied in our simulations because we work with momentum preserving (DPD) dynamics (by contrast, sound is damped in the standard Langevin dynamics). It is noted that an “ideal” open system should be transparent to all waves, meaning that all waves, either created by inner mass fluctuations of any wavelength or by external waves travelling across, should leave the system and do not reflect back. This implies, in particular, that in the absence of external longitudinal forces, fluctuations of the total mass should have no memory, being a white noise (or at least a broad-band signal).⁹¹ The presence of correlations in the total mass of our “region of interest” is in fact due to partial reflections of waves at the rarefied buffers, where density fluctuations are reduced (recall we fix the average total mass of these reservoirs, see e.g.ref. 66). In the autocorrelation (ACF) of the total mass fluctuation, the dominant wavelength should obviously be the largest possible one compatible with this condition and the total simulation box length. A density mode with wavenumber k = 2π/λ decays like 〈ρ(k,t)ρ(k,0)〉 = 〈δρ²〉exp[−ν_Lk²t]cos[ωt] with ν_L being the sound attenuation coefficient⁹² and ω = c_Tk the oscillation frequency. Inspection of Fig. 3 indicates that the density profile of the whole system (MD + buffer) roughly conforms to Dirichlet boundary conditions with a fixed density at the end of the buffer domains ρ(x₂ = ±L₂/2)<ρ₀. In such a case, the longest wavelength available to mass fluctuations of the system would be λ_eff ≃ 2L₂. We fitted the time autocorrelation function of the MD mass to extract ω and compare it with the ansatz ω = c_Tk_eff using ω = c_Tk_eff with k_eff = 2π/λ_eff. The best fit to the simulation results corresponds to λ_eff = 760σ₀, while 2L₂ = 780σ₀ and it is shown in Fig. 4 in terms of ω/k_eff and compared with the sound velocity of melt c_T. The excellent match confirms that the mass in the MD domain has memory induced by reflections of sound waves against the low density domains (buffers). Such reflections could be reduced by coupling the MD with a continuum hydrodynamic field outside,⁶⁶ or by imposing a non-reflecting boundary condition⁹¹ (still to be generalized to MD, see also ref. 92). However, in the present scenario we find that these results are quite interesting because they suggest the possibility of measuring the sound velocity c_T from the fluctuations of the total MD mass. In particular, it might allow to measure how the sound velocity is modified in a sheared melt, c_T(Wi). Although a detailed study would be certainly required, just as an indication, we find that c_T(Wi) (estimated in this way, i.e. from mass fluctuations) decreases.

5 Sheared melt with normal friction

This section presents results for the star molecule model with zero tangential friction between monomers. Kernel and thermostat details (adpd) are given in Sections 3 and 4. We decided to first focus on this model to avoid embarrassing the discussions with details of different cases (open versus closed, normal versus tangential friction) and also because this model presents a richer dynamical behaviour, whose analysis will be useful to understand the effect of tangential friction in Section 8.

5.1 Weissenberg number

The Weissenberg number Wi is a useful number to compare and dissect different regimes in polymer rheology. It is defined as

Wi = τ_rel[small gamma, Greek, dot above], (11)

where τ_rel is the longest molecular relaxation time and [small gamma, Greek, dot above]⁻¹ is the “shear flow time” needed to affinely deform a square box of sheared fluids into a parallelepiped with an angle of 45° between adjacent planes. As stated, in eqn (9) and Table 2, the longest relaxation of our star polymer is related to the molecular rotation so τ_rel = τ_rot in eqn (11). In fact, the diffusion of the CoM of the molecules is much slower (see Table 2). The CoM diffusion does not directly sample the structure of molecules, whose modification under flow is related to the non-Newtonian character of polymers. The Peclet number determines the ratio between CoM diffusion and flow advection Pe = τ_dif[small gamma, Greek, dot above] and for our setup it is about 10 times larger than Wi. In colloidal suspensions, the shear thinning typically starts for Pe > 1 due to shear banding. Interestingly, star polymers constitute a sort of bridge between the open structures of linear polymers and the compactness of colloids. Thus, one might elucubrate that the onset of shear thinning in compact stars could well be due to shear banding (collective molecular ordering in lanes), rather than by (individual) polymer elongations. We shall see later on that both (collective and individual) effects take place in our sheared system, although we advance that the transition to the non-Newtonian regime takes place at Wi = τ_rot[small gamma, Greek, dot above] > 1. Hence at least for the (not so compact) star polymer studied here, shear thinning is not determined by collective ordering at straining rates faster than the CoM diffusion. The hybrid character of star molecules (between colloids and polymers) is the subject of current research.^56–58

5.2 Density and hydrostatic pressure

Denoting the symmetric pressure tensor (exerted by the melt) as P_α,β, the hydrostatic pressure is defined as We find that our model for a star polymer melt expands when sheared under the normal load (see Fig. 17 which collects the results from different cases). This behaviour is consistent with the die swelling phenomenon of polymer melts¹ and the increase in the normal pressure P₂₂ observed in all cases. By contrast, the hydrostatic pressure decreases with Wi in the open domain, while in the NVT box presents a non-monotonous trend (with not large variations). As explained in Section 6 P_iso is influenced by several molecular mechanisms, whose relevance changes with [small gamma, Greek, dot above] along with the molecular structure. The density (in the open domain) and pressures (P₂₂ in NVT and also P_iso in both cases) are quite sensible to temperature changes, as shown in Section 9.

5.3 Rheology

We start by determining the zero shear viscosity, η₀ = η([small gamma, Greek, dot above] = 0). A standard way is the Green–Kubo equilibrium route (via the integral of the time autocorrelation of shear stress).⁹³ To maintain consistency with the non-equilibrium route we follow the experimental approach which fits the trend for η([small gamma, Greek, dot above]) for a range of values of the shear rate with some suitable expression such as the popular Carreau fit.^50,94 We also checked that Green–Kubo viscosities⁹⁵ are similar to the Carreau-fitted ones within statistical uncertainty (about 10%). The Carreau fit has the following expression,

η = η₀[1 + (τ_η[small gamma, Greek, dot above])²]^−β_η/2, (13)

shown in Fig. 5 along with simulation data. Fig. 5 shows the shear viscosity η obtained under open and closed setups for several models with different thermostats and friction kernels (see Table 1). η([small gamma, Greek, dot above]) was calculated from the off-diagonal pressure tensor component P₁₂ = −η[small gamma, Greek, dot above], which was measured in simulations. From eqn (13), the viscosity shear thinning exponent β_η, is such that,

η → [small gamma, Greek, dot above]^−β_η for large [small gamma, Greek, dot above]

and for polymer melts β_η ranges between 0.4 and 1.^1,88 This fit of eqn (13) also provides an estimation of the zero-shear viscosity η₀ and a characteristic time τ_η related to the onset of the shear thinning regime. η₀ differs for different friction kernels of the DPD thermostat (see Table 1). The adpd and sdpdshort thermostats (identical friction kernels) consistently provide the same zero shear rate viscosity η₀ = 0.5 (that does not increase largely for γ_∥ ≤ 10). The sdpdlong model with an increased kernel cutoff R^DPD_cut = 1.5 × 2^1/6σ presents η₀ = (0.60 ± 0.1) + 0.29γ_∥ an increase consistent with the increasing relaxation time with friction† as deduced in the analysis of Kindt and Briels.⁹⁶ Tangential friction (tdpd) leads to η₀ = 2.6 (γ_∥ = γ_⊥ = 1). The meaning of τ_η becomes clear when it is compared with the estimated star rotational relaxation time τ_rot. Here are their values for different thermostats: τ_η = 61 and τ_rot = 59 for adpd; τ_η = 58 and τ_rot = 55 for sdpdshort; τ_η = 125 and τ_rot = 100 for sdpdlong; and τ_η = 287, τ_rot = 332 for tdpd. The error bar of the given values is approximately ± 5. We observe τ_η ≃ τ_rot indicating that the onset of shear thinning, which takes place at [small gamma, Greek, dot above]τ_rot > 1, coincides with the molecular deformation altering the equilibrium rotational diffusion.

Fig. 5 Normalized shear viscosity obtained for several models under open and closed conditions. Details of the model (varying in thermostats and kernels) are given in Table 1 and 2. Lines are best Carreau fits.

We now focus on the adpd model and defer the discussion on the tangential friction to Section 8. The system temperature is fixed to T = 4.0 ± 0.01. In this case, we get a zero shear viscosity η₀ = 0.50 ± 0.05. At larger Wi, the shear thinning exponents β obtained from Carreau fits (eqn (13)) are found to be slightly steeper in the open domain β_η = 0.41(2) than in the NVT box β_η = 0.35(4). Thus at a fixed shear rate, the open system is slightly less viscous than the closed sample. This is in agreement with previous studies for linear and branched melts carried out at a constant and unconstrained density (see Section 10).

The first and second normal stress coefficients Ψ₁ = N₁/[small gamma, Greek, dot above]² and Ψ₂ = N₂/[small gamma, Greek, dot above]² (for the adpd model, discussed in this section) are shown in Fig. 6. For Wi < 20 we find a decrease in Ψ₁ consistent with the Carreau standard behaviour (eqn (13)), providing Ψ₁([small gamma, Greek, dot above] = 0) = 21 ± 1 and an exponent β_Ψ1 ≃ 1.30 ± 0.04 (i.e. N₁ ∼ [small gamma, Greek, dot above]^0.7) quite similar for both ensembles. The relaxation time for Ψ₁ obtained from the fit is also consistent with τ_η = τ_rot within error bars. At Wi > 20, we find a measurable decrease of Ψ₁ with respect the Carreau trend (see Fig. 6), which takes place at slightly smaller Wi in the open case. This corresponds to a loss of the elastic component of the melt at large shear rates. The second normal stress coefficient Ψ₂ = N₂/[small gamma, Greek, dot above]², shown in Fig. 6 is also quite similar in both environments and at large [small gamma, Greek, dot above] scales like Ψ₂ ∼ −[small gamma, Greek, dot above]⁻¹. The similar behaviour for N₂ under open and closed boxes might be due to the fact that we just fix the normal load (in the x₂ direction) and not the hydrostatic pressure as in some other studies,^39,40,42 presenting different trends for N₂ in NVT and NP_isoT constraints.

Fig. 6 First and second normal stress coefficients of the star polymer melt under open and closed conditions and the adpd thermostat for γ_∥ = 1 (blue) and γ_∥ = 5 (purple). Dashed line for Ψ₁ is the Carreau fit (see text).

The normal stress ratio VR ≡ −Ψ₂/Ψ₁ and the recoverable shear strain⁸⁷

SR ≡ Ψ₁/(2η[small gamma, Greek, dot above]) = (P₂₂ − P₁₁)/(2P₁₂) (14)

are standard indicators of viscoelasticity⁸⁷ (e.g. SR vanishes for a Newtonian fluid). As shown in Fig. 7 in our melt model, SR increases with the shear rate, as expected; however, both indicators (SR and VR) clearly show that a change in the elastic component of the model takes place at Wi > 20. Notably, SR decreases, so the melt becomes less compliant with shear strain and stores less elastic energy across the flow–gradient plane, with a jump in the VR. We advance that the amount of the elastic energy loss at large Wi depends on friction forces (notably on the presence of tangential friction as shown in Section 8). This is a clear indication that calibration of friction from detailed all-atom models⁹ is crucial to represent or simulate some particular real melt.

Fig. 7 (top pannel) Normal stress ratio and (bottom pannel) recoverable shear strain index for some model melts under open and closed environments. The results for the case with tangential friction (tdpd) are also included.

5.4 Gyration tensor

The thermodynamic and rheological properties of any polymeric system are intimately related to the deformation of the molecules induced by flow. It is therefore convenient to start by presenting the results for the average gyration tensor of our star polymer model under shear, whose components are shown in Fig. 8. At Wi > 1 the flow induces the alignment of the star molecules in the flow direction and at the same time a compression in the gradient direction. The width of the stars in the neutral direction slightly increases up to Wi ∼ 20 and AT larger shear rates, they also start to contract in this direction. It is also instructive to observe the shape of molecules in their principal deformation axes, obtained from the diagonalization of the gyration tensor. The molecular axes are s₁ = cos(θ_G)x₁ − sin(θ_G)x₂, s₂ = −sin(θ_G)x₁ + cos(θ_G)x₂, and s₃ = x₃ (since G₁₃ = G₂₃ = 0). The ith eigenvalue is noted as G̃_i. The angle θ_G is the average molecular tilt over the flow direction that satisfies:

Fig. 8 Eigenvalues of the gyration tensor scaled with the monomer diameter σ². The results include cases with the adpd friction kernel γ_∥ = 1 (no tangential friction) and the tdpd model γ_∥ = γ_⊥ = 1 discussed in Section 8. Estimated error bars of the data are approximately 5%.

The results for the eigenvalues of gyration tensor are shown in Fig. 8 and reveal the strong contraction of the molecules in the s₂ direction. In particular, at Wi ∼ 50 we observe that G̃₂^1/2 ∼ 1.12σ so the chain width in the contracted axis reaches the monomer (or blob) diameter. In a real polymer at this shear rate, the flow starts altering the structure of the chain “blobs” smaller than σ. But this is precisely our modelling unit length so we will take this as the maximum Wi to be explored here (roughly Wi < 100). Above this (size dependent) shear rate, one should also expect a shear rate dependence of the blob–blob friction kernel (see ref. 97 for a related study on this issue). It is also interesting to note that the total volume of the molecules (evaluated from the product of the eigenvalues of the gyration tensor) first increases (molecules expand) up to Wi ∼ 10 and due to contraction in the neutral direction and then tends to contract (in the form of a highly extended ellipsoid) at larger Wi. In this sense, the behaviour of the arms in the neutral direction of the star have some similarities with what is found in tethered polymers,^98,99 which also present a maximum volume at intermediate Wi.

6 Analysis of the stress–structure and dynamics relationship

The present analysis provides the relationship between stress components and the molecular structure and dynamics. It focuses on the results presented in the previous section, which corresponds to the adpd star model, with zero tangential friction. The conclusions will be useful for the comparison done in Section 8.

6.1 Pressure balance in closed system

Pressure is the leading mechanical variable in rheology as it directly connects the microscopic world with material properties. We now elaborate on the pressure balance in an attempt to connect the molecular structure with the observed rheology.

The pressure exerted by the melt‡ is calculated using the Irving-Kirkwood method and the method of planes,¹⁰⁰ which is particularly suitable for the open setup as it was designed for non periodic boundaries. The pressure can be first decomposed in virial and kinetic parts. The virial pressure includes contributions from spring forces and intramolecular forces (both acting amongst pairs of monomers of the same molecule) and intermolecular forces (between monomers of different molecules). The pressure balance allows to analyze the different molecular contributions,

P_α,β = P^kin_α,β + P^spring_α,β + P^intra_α,β + P^inter_α,β. (16)

The pressure tensor is symmetric for the type of molecules considered hereby (in fact for most polymers).¹ The “average” pressure on the system is given by the hydrostatic pressure which is just the third part of the trace of P (see eqn (12)). The hydrostatic pressure is not involved in (although indirectly affects) the main rheological quantities such as the shear stress P₁₂ and the first and second normal stress differences, defined as

N₁ = P₂₂ − P₁₁, (17)

N₂ = P₃₃ − P₂₂. (18)

It is customary to introduce the traceless pressure tensor

P̂ ≡ P − P_isoI,

with I = δ_α,β the unit second rank tensor. Indeed, the hydrostatic pressure can be decomposed in different molecular contributions with A = kinetic, springs, etc. Also, due to the linearity of the trace operator we can decompose the traceless stress tensor P in a sum of contributions of traceless tensors, i.e., with P̂^A = P̂^A − (1/3)Tr[P^A]I. This is useful to simplify the analysis. In fact, before discussing the differences between open and closed ensembles, we inspect some generalities of the pressure balance focusing on the results obtained for the closed (NVT) ensemble. To that end, we present in Fig. 9 the contributions to the pressure components of the melt under shear to the hydrostatic pressure and the different traceless stress tensors (component-wise: gradient P̂₂₂, flow P̂₁₁, neutral P̂₃₃ and shear P₁₂ obviously P₁₃ = P₂₃ = 0).

Fig. 9 Components of the different contributions to the traceless tensor (, with A= kinetic, springs, etc.) of the traceless stress tensors P̂^A = P^A − (1/3)Tr[P^A]. The results correspond to the closed setup using the adpd thermostat with γ_DPD = 5. P₀ is the (isotropic) pressure at [small gamma, Greek, dot above] = 0, the rest of terms are explained in the text. The error bars of the data are around 0.002–0.005ε/σ₀³.

6.1.1 Spring stress. Bonded interactions are hereby modelled by linear springs with a non zero equilibrium distance. The contribution of the springs of star molecules to the traceless stress tensor is illustrated in Fig. 9 (again, for the case of closed systems). As expected, as the shear rate is increased, molecules tend to be stretched in the flow direction and compressed in the gradient direction. Springs tend to restore the equilibrium (spherical) shape of the molecules by producing a negative stress (compressing) in the x₁ direction and a positive (expanding) in the gradient direction. As in any material with elastic properties, both effects contribute to enlarge the first normal stress difference N₁ = P₂₂ − P₁₁. In turn, we find that elastic contribution to the second normal stress difference N₂ = P₃₃ − P₂₂ is negative, as it happens in real polymeric fluids.¹Fig. 9 shows that the off-diagonal component of springs' contribution to the pressure tensor P^spring₁₂ is the most important contribution to the shear stress, although as discussed below, the kinetic pressure becomes significant at large Wi (above Wi > 20 for the adpd model discussed in this section).

6.1.2 Kinetic pressure. The kinetic pressure tensor was obtained from the peculiar velocities δv_i = v_i − u(r_i), where u(r) = u₁(x₂)e₁ is the average velocity field obtained by binning the velocity gradient direction x₂. Indeed, the ideal part of the hydrostatic pressure is one third of the trace of the kinetic pressure tensor, (1/3)Tr[P^kin] = ρT. However, its traceless part, P̂^kin ≡ P^kin − ρTI, contains relevant information about the correlations in velocity fluctuations over different axes. In particular, velocity fluctuations of monomers in the x₁ and x₂ directions become correlated at [small gamma, Greek, dot above] > 0 and contribute to the shear stress and viscosity. In our setup, this correlation becomes negative (owing to the imposed counter-clockwise shear) and leads to a negative contribution to P^kin₁₂ with a non-negligible viscous contribution to the shear viscosity η = −P₁₂/[small gamma, Greek, dot above] which tends to reduce the shear thinning exponent (see the comparison with the case with tangential friction in Section 8). Although seldom pointed out in the literature (see ref. 4), at high enough shear rates the kinetic pressure of monomers becomes significant in polymer melts. This is in fact what we clearly observe in Fig. 9 for the star model with normal friction (adpd). The anisotropy in the kinetic pressure observed in the covariance of peculiar velocities of monomers is certainly much larger than what would be observed in a simple fluid.^47,78 It is also of a different nature because in the melt, the monomer motions are constrained and they fluctuate and rotate being tethered to the center of their star molecule. This leads to different monomer peculiar velocities in the flow, gradient, and neutral directions which are reflected in the kinetic pressure. In Fig. 9, we show P̂^kinversus Wi (recall that P^kin = P̂^kin + ρTI and diagonal components of P^kin are positive). At Wi > 1, the kinetic pressure contributes to the first and second normal stress N^kin₁ and N^kin₂. As soon as the molecules become stretched in the flow direction (for Wi > 1), velocity fluctuations in the flow direction are enhanced with respect to those along the gradient and neutral directions: P̂^kin₁₁ > 0 and P̂^kin₂₂ < 0, P̂^kin₃₃ < 0. We find that N^kin₂ < 0 (see Fig. 10 below), so the kinetic stress acts in the same way as the elastic (springs) components. However, the kinetic contribution to the first normal stress differences is negative N^kin₁<0 so it is just opposite of the elastic contribution of the chain. Our conclusion is that the kinetic pressure of monomers tends to reduce the elasticity of the melt in the flow–gradient direction at large shear rates.

Fig. 10 Contributions of different mechanisms to the summing up the total normal stress differences, shear stress, and hydrostatic pressure. The results for the adpd thermostat under open and closed conditions. The error bars of the data in graphs are around 0.002–0.005ε/σ₀³.

6.1.3 Intramolecular pressure. Intramolecular pressure gives an indication of excluded volume effects and molecular collisions within one molecule. For our moderate-size molecules, it has a minor contribution to the total stress. Intramolecular collisions induce a viscous (not-restoring) stress in the melt which slightly contributes to the shear stress and tends to counterbalance the elastic first normal stress difference N^intra₁ < 0 (just like the kinetic pressure – see Fig. 10). It has the same effect in the gradient–neutral plane (N^intra₂ > 0) where due to the molecule stretching, monomers tend to collide less in the gradient direction. However, above Wi ∼ 20 a shallow maximum is observed in P^intra₃₃ (and a minimum for P^intra₂₂ – see Fig. 9) suggesting that the reduction of the arm extension in the neutral direction has the consequence of increasing monomer collisions in the gradient direction (thus N^intra₂ diminishes – see Fig. 10). Evidence of arm retraction in the neutral direction is shown in Fig. 8 as a clear the reduction of G̃₃ = G₃₃ for Wi > 20 (also observed in “sdpdshort” and “sdpdlong” thermostats). A somewhat similar rebound in the intramolecular pressure at a large shear rate has been observed in simulations of entangled and disentangled linear polymer melts (see ref. 3 and reference therein). However, in these studies the monomer–monomer interaction is attractive, and this fact modifies (reverses) the intramolecular contributions (polymer compression then reduces the intramolecular pressure, unlike for our purely repulsive potential).

6.1.4 Intermolecular pressure. According to Fig. 9, the contribution of the intermolecular pressure to the melt is minor. The seemingly irrelevant contribution of intermolecular stress is assumed in many theoretical models⁸⁷ to explain viscoelasticity (stress-optic rule). However, this is a simplistic view because intermolecular forces are responsible for spreading the external momentum introduced through the system boundaries (note that, in this respect, the open boundary setup behaves like a real experiment). In fact, internal forces (between monomers of the same molecule) sum up to zero so they cannot modify the CoM velocity of the molecule. For the present star molecule (with relatively short arms), momentum in the flow direction is transferred and maintained across x₂ by intermolecular collisions.⁹ These friction forces gradually build up the elastic stress in the molecules, before finally collapsing to a stationary value in the steady state. The central role of intermolecular forces can be also seen by considering the alternative molecular formulation of the pressure tensor¹⁰¹ (based on molecules (μ) CoM and virial pressure proportional to ). This “molecular pressure” formulation is however less informative because the effects of all internal forces, like the springs, are hidden in the spatial CoM distribution g({R}). Nevertheless, it serves to illustrate the central role of intermolecular forces and shed some light on the apparently striking similarity between the intermolecular pressure and the total pressure dependence with Wi, which can be seen by comparison between the corresponding (inter and total) panels of Fig. 9 (note the difference in values). This similarity between intermolecular (IM) and global magnitudes (such as IM energy and hydrostatic pressure, see Fig. 9) has also been reported as “striking” in previous sheared melt simulations.^2,3 It is interesting to note that the dominant velocity gradient component of the IM force P^inter₂₂ reaches a plateau around Wi ∼ 20 and slightly decreases at larger Wi (as the total P₂₂ does – see Fig. 9). This is indicative of a change in the dynamics of polymers, which according to the concomitant increase observed in P^inter₁₁, most probably start to rotate and collide more often in the flow direction.

6.1.5 Hydrostatic pressure. The hydrostatic pressure is key in shear induced polymeric phenomena, such as shear induced crystallization⁸⁷ or separation of blends.⁴⁷ Its dependence on the shear rate is not well understood. P_iso depends on the molecule size^40,43 and their architecture^3,40,42,88. However, there have been reports of both increase and decrease in different (sometimes contradicting) studies. Fig. 9 indicates that P_iso can present a non-monotonous trend with Wi due to changes in the molecular structure under shear. Non-monotonous trends for P_iso([small gamma, Greek, dot above]) have also been observed for different polymers in ref. 2, 3, 40 and 43. The present analysis (see Fig. 9) reveals in fact that P_iso depends on a competition of several mechanisms. Below Wi ≃ 10 the hydrostatic pressure varies little but presents a slight descend, probably due to the chain expansion (revealed by the analysis of the gyration tensor) which decreases the intramolecular pressure. At larger shear rates Wi > 10, two opposite mechanisms enter into the play: First, an increase in P₁₁ due to the strong increase of kinetic pressure (and to a lesser extent to intermolecular collisions). This is probably due to molecular rotations similar to the tank-thread motion reported for a star molecule solution¹⁰² and also pointed out in ref. 2 and 3. And second, a decrease in P₃₃ due to the contraction of the stars in the neutral direction (see Fig. 8) and consequent reduction of the kinetic pressure P^kin₃₃. Both effects nearly counterbalance each other (see Fig. 9) leaving small variations in P_iso. The relevance of these effects depends on the boundary conditions (the open case is analysed in next section) and type of friction (Section 8). More generally variations of P_iso depend on the presence of attractive monomer interactions and the molecular structure.

It is interesting to note that the decrease of neutral kinetic pressure in our adpd star model starts to take place around Wi ≃ 20, which is precisely the ratio τ_rot/τ_elas ≃ 20. Above this shear rate the flow strains faster than the elastic relaxation of the molecules thus reducing the fluctuations of the arms in the neutral direction. This effect finally induces a net decrease in P_iso above Wi > 30. We shall see (Section 10) that this effect is absent (or at least delayed) while introducing the tangential friction, again indicating that friction should be an essential part of any CG model.⁹

6.1.6 Intermolecular forces and the Hookean limit. To summarize, at low shear rates the elastic energy stored by the melt grows in response to intramolecular (non-bonded) interactions (here mainly friction forces). At large enough shear rates viscous (Newtonian) effects coming from the kinetic and intra molecular pressures, tend to modify (normally reduce) the elastic response of the melt (notably first normal stress differences). Our findings are in agreement with the conclusions of Kroger et al.⁴ clearly and succinctly summarized in their book (page 144) in relationship with the breakdown of the linear stress-optic-response (SOR) due to the Newtonian viscous transport at large Wi. In our simulations we find a linear relationship between G₂₂ − G₁₁ and P₂₂ − P₁₁ (essentially similar in open and closed setups), indicating the validity of the linear SOR up to Wi ≃ 10. In the approach by Kroger et al. (applied to linear multibead FENE chains), both, kinetic and intramolecular pressures, are collected in what they call the “simple fluid” stress contribution.^4,89 We find here that both contributions can have different roles, which more generally should probably depend on the molecular shape (kinetic pressure) and monomer (intramolecular) interactions (attraction/repulsion).

To highlight the relevancy of the intermolecular stress in the linear SOR regime, we plot in Fig. 11 (top panel) the total and elastic contributions (springs) to the normal and shear stress against the intermolecular counterparts. At low shear rates, about Wi < 10, the same linear relation is found for the first normal stress difference and shear stress, N₁ ≃ 8N^inter₁ and P₁₂ ≃ 8P^inter₁₂, while we found N₂ ≃ 4N^inter₁. Approximately the same linear relation holds for the elastic stress and also for the normal molecular strains evaluated with the gyration tensor (G₁₁ − G₂₂ and G₂₂ − G₃₃, scaled in Fig. 11(top panel)). This provides the Hookean limit of the melt N₁ ≃ C(G₁₁ − G₂₂) and N₂ ≃ (C/2)(G₃₃ − G₂₂) with C = 62.5 (we find that P₁₂ ∼ CG₁₂ holds only for smaller Wi < 5). In the open (adpd) setup, the linear regime for N₁ and P₁₂ perfectly agrees with closed simulations; however we found deviation from linearity in the case of N₂, an issue which deserves further investigation. The bottom panel of Fig. 11 presents the results for stars with the tangential friction, analysed in Section 10.

Fig. 11 Total stress, elastic stress components, and molecular strain differences in normal planes measured from the average gyration tensor (G₁₁ − G₂₂ and G₃₃ − G₂₂) plotted against the intermolecular stress component corresponding to each case (i.e., N^intra₁, N^intra₂, etc.). The results for adpd (top panel) and tdpd models under open conditions (bottom panel). Estimated error bars of the depicted quantities are approximately 5%.

7 Effect of open environment

7.1 Density and hydrostatic pressure

In our open domain, we fix the load of the melt in the gradient direction (P^ext₂₂ = P₂₂) and this produces a redistribution of the pressure tensor, reducing its component in the flow and neutral directions and also, indirectly, its shear stress. This is deduced from Fig. 12 where we compare the traceless stress tensor and the hydrostatic pressure in the open and closed setups at the fixed T = 4 temperature. Note that P^ext₂₂ = P₂₂ = P̂₂₂ + (1/3)P_iso is constant in the open setup (within statistical uncertainties). In the open domain, the sheared melt expands in the gradient direction; a phenomenon similar to the die swelling observed in polymer extrusion at the pipe orifice and related to other viscoelastic phenomena.¹⁰³ In the open domain, this corresponds to a decreasing melt density (at faster shear rates) and brings about a smaller hydrostatic pressure than in the closed environment at similar Wi (see Fig. 12). However, the relative decrease of P_iso is larger than the density jump. This fact is due to several related effects we now analyze from the inspection of Fig. 10. Indeed, at the fixed temperature, a lower density leads to a lower kinetic pressure P^kin_iso = ρT found in the open domain (this trend also applies here to the intramolecular pressure P^intra_iso because our model considers purely repulsive nonbonded interactions, the opposite effect could arise for attracting chains^2,3). However, an even larger reduction in P_iso with respect to the closed box is observed due to the smaller intermolecular pressure in the open box (see Fig. 10 for Wi > 10). Indeed, at a high shear, a less dense melt presents less molecular collisions, less intermolecular friction and thus less elastic load. As stated, at the fixed temperature, the elastic strain is essentially activated by intermolecular friction in the melt. Notably, at Wi ∼ 20, these intertwined effects induce a reduction of about 25% of hydrostatic pressure of the open domain mainly arising from the decrease of elastic stress. At this Wi, density has only decreased about 10% (see Fig. 17). In agreement with this comment, we note that this depressurizing effect is doubled when the tangential friction is added as commented in Section 10.

Fig. 12 Hydrostatic pressure and components of the traceless stress obtained in closed and open domains at fixed T = 4 and under the adpd thermostat (normal friction). Error bars of the measured quantities are approximately 0.005ε/σ₀³.

7.2 Normal stress differences and shear stress

Fig. 10 compares the different contributions to the normal stress differences and the shear rate in the open and closed systems. Much of what has been already said in the previous section applies here. The results for the shear stress nicely corroborate what we pointed out before about the close relationship between intermolecular and elastic stresses. Fig. 10 (right bottom) clearly shows that the kinetic and intra molecular shear stresses are essentially equal in the open and closed domains. The decrease in elastic shear stress found in the open case is due to the reduction in the intermolecular friction at a lower density (although seemingly paradoxical, the “nominal” contribution to P^intra₁₂ in the monomer pressure balance is minor). In all instances, at low enough shear (Wi < 10) the elastic stress is close to the total stress (and proportional to the intramolecular stress). This is the regime of validity of the stress optic rule which is broken at larger shear due to viscous (and compressible) effects related to the “simple fluid” of monomers. Remarkably, for the present model of the star polymer melt, the kinetic normal stress becomes the dominant “viscous” contribution and at large Wi > 20 it even induces a decrease in the first normal stress difference N₁.

7.3 Molecular ordering under shear

Fig. 13 presents the angle of the largest eigenvector of different contributions to the pressure tensor with the flow direction (measured according to eqn (15)). We also include the molecular orientation, measured from the angle associated with the gyration tensor (eqn (15)). This plot condenses what has been already mentioned in previous sections. Recall that a spherical molecular structure provides θ_G = 45° and similarly from eqn (15) a Newtonian fluid without elastic component, necessarily presents θ = 45°. This is what is observed for the total pressure tensor angle θ_P and the molecular orientation θ_G at Wi → 0 in Fig. 13. As the straining rate is made faster both angles decrease in a similar fashion, however, at Wi > 10, the pressure tensor angle θ_P presents a minimum and starts increasing towards 45°. By contrast, the molecular orientation θ_G keeps aligning with the flow direction. This in an indication of the loss of the Hookean behaviour of the melt which here is mainly due to the kinetic pressure (see its principal direction in Fig. 13). The stress direction of springs also aligns with the flow, although their angle is larger than the molecular orientation (a similar outcome was observed in ref. 3). Finally, note the close match between the direction of intermolecular forces and the total pressure tensor. As stated, intermolecular forces are the driving mechanism of transformation between the viscous flow and elastic energy. Lastly, as observed by other authors,^39,42,43 the open boundary does not modify the molecular structure or orientation with Wi, when compared with the closed case. Here, Fig. 13 shows another remarkable result: the significant redistribution of pressure in the open case (see Fig. 12) does not alter the orientation of the different contributions of the pressure tensor at the increasing shear rate. The orientation of pressure eigenvectors is a function of their eigenvalues which in turn determine the material properties of the polymer (its viscoelasticity⁸⁷). A relevant example is the recoverable shear strain (SR) given by eqn (14). Material properties should not depend on the constraints used to perturb the polymer and this is precisely what our analyses provide.

Fig. 13 Top panel: angle between the principal eigenvectors of the contributions to the stress tensor and the principal eigenvector of the average gyration tensor. Bottom panel: comparison with the model including transversal friction (tdpd model). Angles are with the same color coding.

To observe the collective order of the star molecular, we calculated the CoM pair distribution function g(R_ij). Fig. 14 illustrates the marginal distributions for different planes and at increasing shear rates. It is illustrative to draw the directions of the principal components of the pressure and the gyration tensor to observe the departure from the Hookean (linear SOR) regime. Above Wi > 20 molecules start to orient in lanes in the flow–gradient plane, as indicated by the elongated shape of the CoM distribution in Fig. 14. In linear polymers at large shear rates, this effect creates order and even crystallization (ref. 87). This collective order is also clearly visible in one snapshot of the system, in Fig. 15. Molecules with different relative velocities (closeby in the gradient direction) slide over creating a stress which is larger along a direction differing from the orientation of individual molecules. This direction of maximal stress is correlated with the CoM distribution (see Fig. 14) which shows bright spots where molecules slide over and depleted regions in the “wake” of each molecule. At the largest shear rates considered, we also observe some collective ordering in the neutral (vorticity) direction. In sheared colloids, lanes of particles in the neutral direction appear due to hydrodynamic interactions.¹⁰⁴ This could be a plausible hypothesis, in view of the hybrid colloid–polymer nature of the star molecule direction.

Fig. 14 The pair distribution of molecules, g(R) providing the probability of finding the CoM of a molecule at a distance (vector) R from the target molecule CoM. Left: Marginal probability in the flow–gradient and (right) gradient–neutral planes. The results for increasing Wi in a closed domain (open simulations at similar Wi are visually indistinguishable). Green line denotes the direction of the deviatoric stress and the blue line the molecular orientation obtained from the gyration tensor eqn (15).

Fig. 15 Snapshot of a closed box simulation under shear rate (Wi = 42). For the sake clarity only 10% of randomly chosen molecules are shown. Colors indicate the position in the neutral direction. The snapshot clearly shows the formation of lanes tilted in the flow–gradient plane created due to sliding over rows of oriented molecules. Compared to Fig. 14, the snapshot is rotated for 180° with respect to the gradient axis.

8 Effect of tangential friction

To investigate how the tangential friction between monomer blobs alters the rheology of the melt model, we use a standard (i.e. not adaptive) DPD thermostat⁷⁰ with γ_∥ = 1 and γ_⊥ = 1 (see Section 3, 4, and 5.3 for details). As stated the friction kernels are Heaviside functions, in this case with a cutoff distance R_dpd = 1.5 × 2^1/6σ. In the following, we label this tangential friction thermostat as “tdpd”. This tdpd thermostat was found to be strong enough to keep the relative increase of the system temperature smaller than 5% at the largest shear rates considered.

In a real system friction acts by reducing the relative velocities of interacting monomers, generally in the normal and also in the tangential directions. Under Markovian and pairwise interaction assumptions, this form of friction leads to the DPD equations as shown by the Mori–Zwanzing dynamic coarse graining applied to the microscopic Liovillian dynamics.⁹ The same effect is properly captured by the tdpd thermostat, although here at a qualitatively level. Reducing the monomer relative velocities immediately leads to a reduction in kinetic pressure which has large consequences on the system rheology. In particular, the behaviour of the melt is essentially ruled by its elastic component, activated by the more effective intermolecular friction. Just to illustrate this point, we plot in Fig. 16 the contributions of the first normal stress difference (a direct measure of viscoelasticity) for the tdpd case. Compared with Fig. 10 (for the γ_⊥ = 0 adpd case) the tdpd model has a much smaller kinetic pressure and N₁ is essentially determined by the elastic stress (particularly as Wi increases). The same conclusion applies to P₁₂ in Fig. 16.

Fig. 16 Balances for the first normal stress difference and the shear stress for the tdpd model including tangential friction. Comparison with Fig. 10 for the zero tangential friction case. Error bars of the measured quantities are approximately 0.002–0.005ε/σ₀³.

An interesting difference related to the presence of tangential friction concerns shear dilatancy. Fig. 17(a) presents the relative density expansion δρ/ρ₀ = 1 − ρ/ρ₀ for different models. Let us now focus on the adpd and tdpd models which are kept isothermal (non-isothermal cases are discussed in Section 9). The density expansion of the adpd model (without tangential friction) scales like δρ/ρ₀ ∼ Wi, while tangential friction (tdpd) leads to a much softer trend δρ/ρ₀ ∼ Wi^0.5 (although it expands relatively more at moderate shear rates). Under a constant normal load, shear dilatancy is a consequence of the growth of pressure in the velocity–gradient direction. In the case of small kinetic effects (tdpd) this growth is controlled by the expansion force arising from the compressed springs. This elastic pressure appears as soon as molecules start to align with the flow and compress in the gradient direction. Under enough tangential friction, the elastic stress is dominant and also controls the hydrostatic pressure, which in the absence of kinetic pressure effects, presents a faster decrease at large Wi compared to the adpd case (see Fig. 18). Of course, this decrease is also related to the fact that the tpdp simulations were done in the open system; notably for tdpd we get about 50% reduction in hydrostatic pressure for less than 10% reduction in density (see Fig. 17).

Fig. 17 Top panel: Relative density variation for the different cases studied (open and closed, and different models and thermostats, see Table 1). Bottom panel: relative temperature increase observed in simulations with standard DPD thermostats with vanishing tangential friction. The relative variation of T for the adpd thermostat is smaller than 0.01 and less than 0.045 for the tdpd case (both of them not shown in the graph). Details of the thermostats are given in Table 1. The lines correspond to eqn (21), with the characteristic constant A defined there. Error bar estimations of the data are approximately 0.005 and 0.01 for the upper and lower graph, respectively.

Fig. 18 Hydrostatic pressure and components of the traceless stress for the adpd model (normal friction between blobs) compared with the tdpd model (including tangential friction). Both cases in open boxes at fixed T = 4. Error bars of the measured quantities are approximately 0.002–0.005ε/σ₀³.

The eigenvalues of the gyration tensor shown in Fig. 8 also indicate that adding tangential friction makes star molecules “stiffer”, in the sense that one needs larger values of the Weissemberg number to deform them. This observation is however somewhat misleading because for a fixed Wi, the real (physical) shear rate [small gamma, Greek, dot above] = Wi/τ_rot is now smaller due to the increase in τ_rot with the friction. In any case, the tangential friction is expected to alter the stress–strain relations in the Hookean regime (related to the linear stress-optic rule coefficient). This is (indirectly) seen in Fig. 11 where we plot the normal strain differences G₁₁ − G₂₂ and G₂₂ − G₃₃ (also G₁₂) against the corresponding intermolecular stress differences (against P^intra₁₂). We choose this plot to illustrate two facts: first, if the kinetic pressure is minor, the intermolecular friction is the leading mechanism driving the molecular deformation (and its elastic response). Second, molecular strains (and elastic stresses, not shown in the figure) in the flow–gradient and gradient–neutral planes (corresponding to the first and second normal stress differences) present a quite similar linear scaling (the shear deformation as well) with the intermolecular stress. This is to be contrasted with the top pannel of Fig. 11 (model in the absence of tangential friction), where the second normal stress (and G₂₂ − G₃₃) is half of the first N₁ counterpart (also G₂₂ − G₁₁). A perfect alignment between intermolecular forces, elastic stress and molecular strain was also found by Kroger et al.⁸⁹ in linear FENE chains. It thus seems that tangential friction (tdpd) helps to reduce the second normal stress in such a way that N₁ and N₂ present similar scaling laws N₁ ∼ [small gamma, Greek, dot above]^0.68±0.02 with N₂ ∼ −0.3N₁. This is to be compared with the adpd case in Fig. 6 (N₁ ∼ [small gamma, Greek, dot above]^0.70 and N₂ ∼ [small gamma, Greek, dot above]^1.0). As shown in Fig. 7 the tdpd model yields −N₂/N₁ ≃ 0.3 at any Wi < 100. This value is characteristic of disentangled melts (being −N₂/N₁ = 2/7 the theoretical prediction for small shear rates^1,88).

The monotonous increase of elastic storage with [small gamma, Greek, dot above] found in the tdpd model is reflected in the recoverable shear strain (SR) shown in Fig. 7. Somewhat paradoxically, adding tangential friction increases the melt elasticity. Albeit, this reinforces the conclusions in Section 6: the intermolecular friction is the principal mechanism loading elastic stress into an disentangled melt. In passing, we note that in the tpdp model the orientational resistance parameter m_G = Wi tan(2θ_G) grows like m_G = 3.7Wi^0.65 (at least for Wi < 100), a scaling which agrees with that reported for stars in solution¹⁰⁵ (the prefactor being however about twice larger in our melt).

Not unexpectedly, the zero shear viscosity for the tdpd star model is larger η₀ = 2.6 than the adpd case η₀ = 0.5. The relaxation time is also larger τ_η = 287 (compared with 60). However, the tdpd shear viscosity thinning is faster; we find β_η ≃ 0.5 compared with 0.4 for the adpd case (recall η ∼ [small gamma, Greek, dot above]^−β_η). Again, this is also a consequence of a much smaller contribution of the viscous stress coming from kinetic effects. If the tangential friction is absent, the kinetic (and intramolecular) contribution increases the shear stress and the viscosity at any Wi leading to a softer shear thinning exponent.

9 Thermostats and heat dissipation

9.1 Density and temperature

We now briefly analyze the effect of the temperature increase due to heat dissipation in the sample. In all cases, the temperature reaches a steady state, but as shown in Fig. 17 (bottom panel) plotting δT/T₀ = T/T₀ − 1, we face severe viscous heating while performing the first row of simulations with the sdpdshort and (stronger) sdpdlong thermostats (see Table 1). Heating is observed at Wi > 10 and irrespective of the damping parameter (we tried up to γ_DPD = 50m₀τ₀⁻¹). Tangential friction drastically reduces heating (δT/T₀ < 0.045 for Wi < 70) however the adpd thermostat enabled us to simulate zero tangential friction at a fixed temperature, providing δT/T₀ < 0.01.

Let us focus on the “heated” runs at the increasing shear rate to illustrate the effect of uncontrolled temperature. Fig. 17 shows that heating introduces further melt expansion under shear and this kinetic energy induces larger hydrostatic pressure (which it is seen to increase with shear in the sdpdshort case). Pressurization due to viscous heating can also alter the rheology response. This is seen in Fig. 5 where the sdpdshort case presents shear thickening at Wi > 10, but only under closed conditions. The shear thickening reported in some of the published studies on polymer melts (closed box simulations) might be due to viscous heating (see e.g.ref. 39). More interesting than this elucubration is the results of Fig. 5 for the sdpdshort-open. The viscosity obtained for the same thermostat in an open environment is not affected by the temperature increase with [small gamma, Greek, dot above]. In fact it presents the very same trend as the adpd case. This observation indicates that shear viscosity is dominated by the normal load. In fact, the same outcome was also observed for the sdpdlong model (with shear exponent β_η = 0.39) and for all γ_∥ ≤ 10 considered. Thus, this insensitivity of viscosity to temperature found only under normal load is not probably due to a cancellation of effects sometimes observed in experiments⁵⁰ (viscosity decreasing with T and increasing with P₂₂). Rather it should be due to the viscosity dominated by pressure as happens in highly pressurized melts.⁵⁰ Here, monomers interact via purely repulsive forces (WCA) which might map a sample under large pressure, in fact, adding attractive interactions would probably trigger temperature effects on viscosity.

9.2 Viscous heating

The rate of heat production per unit volume due to viscous dissipation is Q̇_η = η[small gamma, Greek, dot above]² leading to a larger steady temperature whose value depends on the heat extraction rate. The onset of temperature increase is usually determined by a non-dimensional parameter which depends on Wi and on the rate of cooling Q̇_c (see the recent computational study in ref. 10). Although in this work we shall not focus on heat and entropy production, we believe it is interesting to share our observations on this phenomenon, partly because of the relatively small simulation literature accurately reporting heating effects in sheared, thermostatted melts. A simple equation for the heat produced in the sheared melt includes frictional gain and cooling,where c_X is the specific heat capacity (molar) at a constant pressure (X = p) or volume (X = V) and dQ_c = ρ_nc_XdT (here ρ_n is the monomer number density). The DPD thermostat extracts (kinetic) energy upon pair collisions, at a rate which is proportional to the temperature difference T − T₀, where T is the temperature of the system (kinetic) and T₀ the nominal temperature of the thermostat.§ The value of cooling rate α (which has units of the number of colliding pairs divided by volume and time) should scale aswhere we note that w(r) is the DPD kernel, which is simply a Heaviside function in our thermostat. The heat cooling rate of the thermostat increases with its damping coefficient γ_DPD, the kernel cutoff R^DPD_cut, and with the square of the monomer local density, which determines the number of thermalizing monomer collisions. In the steady state Q̇ = 0 we thus have, Using the viscosity trend η = η(Wi) obtained from simulations (see Section 5.3 above), we plot prediction eqn (21) in Fig. 17 (lines) for different cases considered. The agreement is quite reasonable, indicating that the temperature increase in the melt can be forecast using a simple thermodynamic argument. Best fits to eqn (21) provide A ≡ αc_P = 0.011 in the open-sdpdshort thermostat, while somewhat larger A = 0.017 for the closed-sdpdshort thermostat. For the closed sdpdlong thermostat, which has about twice as much colliding partners within the dpd kernel, we consistently get A = 0.03. A preliminary calculation of c_X from the variance of the system energy in a closed (NVT) equilibrium simulation c_V = 〈δU²〉_NVT/(NT²) provides c_V ≃ 1.7. From eqn (21) we get the same order of magnitude (A ≃ 0.05 for sdpdshort and A = 0.09 for sdpdlong). A better agreement is found while comparing with the tdpd (sdpdlong) thermostat whose best fit provides A = 0.06(5) against the prediction 0.09. This indicates that eqn (20) should also depend on the number of degrees of freedom the thermostat acts upon (3 in the case of tdpd, 1 otherwise). A more refined calculation would also require including the dependence of c_X and g(r) on the shear rate.

10 Comparison with previous studies

It is interesting to compare our results with previous rheological studies, some of them were carried out at isobaric (constant P_iso) or constant load (P₂₂) constraints. As stated in the Introduction, the number of studies of flowing melts under a constant pressure (either hydrostatic pressure or normal load) is not large. However, they present significant discrepancies in the density and pressure variations with shear. For instance, ref. 42 presents the results for dendrimer melts under the isobaric conditions (constant hydrostatic pressure) revealing a decrease in the melt density under shear. For linear chains, ref. 40 presents just the opposite effect (contraction under shear), while ref. 43 (constant load) reports shear expansion (density increase).

The shear thinning exponents found here for a star polymer melt are consistent with those found in other simulations for somewhat similar systems such as hyperbranched and dendrimer polymers.^41,42 Recall that the shear thinning exponent of any quantity Φ is β_Φ with Φ ∼ [small gamma, Greek, dot above]^−β_Φ at large [small gamma, Greek, dot above]. For viscosity we find β_η ≃ 0.4 (adpd) and β_η ≃ 0.5 (tdpd), for the first normal stress coefficient β_Ψ1 ∼ 1.30 (adpd) and 1.31 (tdpd), while for the second one the values are β_Ψ2 ∼ 1.0 (adpd) and 1.31 (tdpd). For dendrimers, Bosko et al.⁴¹ reported shear thinning exponents increasing with M under NVT (β_η ∈ [0.28–0.36]) while, for NP_isoT the were roughly independent of M (β_η ∈ [0.37–0.39]). The same work reported β_Ψ1 ≃ 1.27 and β_Ψ2 ≃ 1.23 under NP_isoT, while β_Ψ1 ≃ 1.1 β_Ψ2 ≃ 1.0 under NVT. Closed simulations for hyperbranched polymer melts⁴² predicted β_η ≃ 0.3 (slightly increasing with molecular mass M), while β_Ψ2 ≃ 0.95 and β_Ψ1 ≃ 1 (both roughly independent on M). The scaling of (first and second) normal coefficients are probably sensitive to the type of external constraint (either isobaric P_iso or constant load P₂₂). As an indication, a numerical study at constant load for linear chains⁴³ reported the same exponents found in this work (β_Ψ1 = 1.35 and β_Ψ2 ≃ 1 respectively). Concerning shear viscosity, if the density is allowed to relax (under isobaric or constant load conditions) the shear thinning exponents reported in the literature show much less variation with the molecular weight than under a fixed density. This has been shown for linear chains under a constant load⁴³ and also for dendrimers.⁴² In both cases, the shear thinning exponent in the variable density case was found to be close to 0.4, while it ranged from 0.25 to 0.4 in the closed system (increasing with the molecular weight in both cases). This larger insensitivity of shear thinning exponents under fixed load conditions agrees with our observations in simulations presenting viscous heating mentioned in Section 9.

More recently, very similar papers^106,107 studied the rheology and dynamics¹⁰⁸ of star polymers with different functionalities in solution. They reported that the contribution of the stars to the sheared solution viscosity scales like Wi^−0.4 which is quite close to the shear thinning exponents we find here for the melt. According to these simulations,^106,107 the first normal stress coefficient scales like Ψ₁ ∼ Wi⁻¹ in solution, although the authors claim an exponent of −4/3 at very high shear rates Wi > 100. In the melt, we observe the −4/3 power law at smaller shear rates. The second normal coefficient seems also to scale slightly differently in solution Ψ₂ ∼ Wi^−4/3 than in the melt Ψ₂ ∼ Wi⁻¹. The strong similarities in solution and melt indicate that rheological properties are mainly ruled by conformational changes of the chains and maybe also that the hydrodynamic effects in the melt are somewhat similar to those in a liquid solvent. Finally, it is interesting to note that the range of values for β_η and β_Ψ1 for stars, dendrimers and hyperbranched molecules are consistent with the theoretical calculation of Kroger et al.⁸⁸ based on the Fokker–Planck equation for the bond vector distribution of multibead linear chains having slow rotational diffusion (compared with the entanglement relaxation) and a finite deformation energy.⁸⁸ These two features are indeed consistent with the nature of short-armed stars, dendrimers and hyperbranched molecules with internal excluded volume interactions.

The number of experimental studies on rheology of sheared star polymers is not abundant either, but it is growing fast due to the foreseen technological applications.^52,53 Most of these studies consider long armed stars which present entanglements. However, the shear thinning exponent we found for the tangential friction case β_η = 0.5 is quite close to the star case of the experiments by Tezel et al.¹⁰⁹ (4 arms stars with M_a = 140 kg mol⁻¹; the lowest molecular weight considered in these experiments). Although entanglements in star polymers are not yet fully understood in the non-Newtonian regime,⁵² they should be responsible for the substantially increase in shear thinning (exponents close to 0.8) observed in the recent experiments with star molecule melts of Snijkers et al.⁵² Shear thinning exponents found in melts of linear chains are around β_η = 1 (predicted by reptation theory). An open question is why star molecules significantly reduce the shear thinning.⁵²

The analyses presented in Section 6 based on the (exact) balance of the pressure coming from the different contributions (kinetic, intramolecular, intermolecular and bonds) provide some clues which can be useful in a broader study. These types of analysis were performed by Baig et al.,^2,3 Matin⁴³ and Tschopp.⁴⁵ However these balances were based on the energy budget which has less direct rheological consequences than the stress. Also, surprisingly, the kinetic contribution was not considered (leading to unbalanced analyses). The sole exception we found is the work of Kroger et al.⁸⁹ for linear FENE chains which, like ours, is based on the exact stress budget. In elongational flow of linear chain melts, Kroger et al.⁸⁹ found that viscous effects came mainly from intramolecular collisions (they called it “simple fluid” stress offset). Here, we find that in star molecules the kinetic pressure can be the dominant viscous contribution to the melt rheology at large shear rates (a possibility in fact recognized by Kroger et al. in their book⁴). The kinetic pressure is however reduced with the tangential friction and this warns about the need of dynamic coarse-graining⁹ to represent a realistic model. While linear chains tumble in shear flow,¹¹⁰ star molecules perform a quite different motion called tank-threading^102,105 (whereby arms rotate around the molecule center). One can speculate that the reduction of shear thinning observed in star melts is due to an enhanced kinetic pressure related to the tank-threading. Above Wi > 10 we observe the tank-threading motion in our melt with monomer angular momentum growing much faster (ω ∼ [small gamma, Greek, dot above]^0.6 from the preliminary results) than it does in stars in solution (where ω is constant^102,111). Concerning intermolecular interactions, quite often neglected in the literature (see e.g.ref. 87), we find that these are the key to transfer the externally imposed stress into molecular strain and elastic stress. This observation agrees with that of Kroger et al.⁸⁹ Intermolecular friction is the sole possible mechanism if entanglements are not relevant (like in our short arm star molecules).

11 Conclusions

We have conducted OBMD simulations of the melt of star polymers (73 monomers, 12 arms of 6 monomers per arm) using a relatively new modelling technique which combines the adaptive resolution and open-MD (respectively, introduced in ref. 22 and 46 and used in many other studies). From a technical point of view, this work substantially enlarges the size of molecules exchanged between the open system and the reservoir.¹³ In this work, we however focus on what happens if a melt is sheared under a constant normal load instead of a constant volume. The constant load is in fact the condition in many experiments (see ref. 40 and 44) although the simulation literature on the subject is not abundant. We have also presented some conclusions on how the tangential and normal monomer–monomer frictions affect these CG molecular models.

The OBMD permits to perform several new features in a molecular simulation. At equilibrium, the OBMD correctly represents the mass fluctuations of the grand canonical ensemble and interestingly, it could permit to study how fluctuations and sound velocity are altered under non-equilibrium (e.g. sheared) states. OBMD also allows to impose the external shear stress P^ext₁₂ as required for (flux based) hybrid continuum-MD simulations.^15,66 It could also enable the validation of theories like extended thermodynamics⁴⁷ predicting different outcomes for the “conjugate” constant stress and constant shear rate [small gamma, Greek, dot above] non-equilibrium ensembles.

Concerning the present study, we observe that under a constant normal pressure, the melt expands when sheared (shear dilatancy) leading to substantial depressurization and slightly decreasing the shear viscosity. This behaviour was observed in most previous simulations on sheared melts, but surprisingly it is still unclear if it is the general (universal) trend (e.g. see ref. 40). This study provides new information on the rheology of sheared melts: notably, we see that the type of monomer friction is a key aspect for the system rheology. From a theoretical standpoint,⁹ friction between monomers (or rather “polymer blobs”) arises as soon as one consider a CG view of the complex molecule (which is the standard case in polymer science). The pressure balance analysis reveals that in the absence of tangential friction, the monomer kinetic stress becomes significant at large shear rates, increasing the system viscosity (reducing the shear thinning exponent) and diminishing its elastic response (e.g. normal stress difference and stress recovery). We also observed viscous heating in some simulations (e.g. using the sdpdshort and sdpdlong thermostats) revealing a viscosity jump (shear thickening) in closed systems. By contrast, the viscosity of the heated (and less dense) open samples did not change in trend (shear thinning) under a constant load. This indicates that the melt viscosity is controlled by the normal pressure, at least for the present type of molecules with purely repulsive interactions and no significant entanglements. We expect the OBMD to become useful in other studies, such as adding the energy transfer⁶³ or extending the incompressible coupling in ref. 10 to perform hybrid simulations of compressible melts including the transfer of dissipated heat through/across system boundaries. Thus, the properties related to heat conduction could be investigated.¹¹²

Acknowledgements

We thank J. J. Freire for initial configurations and K. Kremer for useful discussions. J. S. and M. P. acknowledge financial support through grants P1-0002 and J1-4134 from the Slovenian Research Agency. J. S. acknowledges financial support from Slovene Human Resources Development and Scholarship Fund (186. JR). R. D-B acknowledges support from the Spanish government under national MINECO project FIS2013-47350-C5-1-R. Partial support from the COST Action MP1305 is kindly acknowledged.

References

1. R. Larson and H. Brenner, Constitutive Equations for Polymer Melts and Solutions: Butterworths Series in Chemical Engineering, Elsevier Science, 2013.
2. C. Baig and V. G. Mavrantzas, J. Chem. Phys., 2010, 132, 014904.
3. C. Baig, V. G. Mavrantzas and M. Kröger, Macromolecules, 2010, 43, 6886–6902.
4. M. Kröger, Models for Polymeric and Anisotropic Liquids, Springer, Berlin Heidelberg, 2005.
5. J. Padding, E. van Ruymbeke, D. Vlassopoulos and W. Briels, Rheol. Acta, 2010, 49, 473–484.
6. V. A. Harmandaris, N. P. Adhikari, N. F. A. van der Vegt and K. Kremer, Macromolecules, 2006, 39, 6708–6719.
7. C. F. Abrams, J. Chem. Phys., 2005, 123, 234101.
8. P. Ilg, H. C. Öttinger and M. Kröger, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2009, 79, 011802.
9. C. Hijon, P. Espanol, E. Vanden-Eijnden and R. Delgado-Buscalioni, Faraday Discuss., 2010, 144, 301–322.
10. S. Yasuda and R. Yamamoto, Phys. Rev. X, 2014, 4, 041011.
11. R. Delgado-Buscalioni, K. Kremer and M. Praprotnik, J. Chem. Phys., 2008, 128, 114110.
12. R. Delgado-Buscalioni, K. Kremer and M. Praprotnik, J. Chem. Phys., 2009, 131, 244107.
13. R. Delgado-Buscalioni, J. Sablić and M. Praprotnik, Eur. Phys. J.: Spec. Top., 2015, 224, 2331–2349.
14. R. Delgado-Buscalioni, Numerical Analysis of Multiscale Computations, 2012, 145–166.
15. G. De Fabritiis, R. Delgado-Buscalioni and P. V. Coveney, Phys. Rev. Lett., 2006, 97, 134501.
16. K. Mohamed and A. Mohamad, Microfluid. Nanofluid., 2010, 8, 283–302.
17. J. H. Walther, M. Praprotnik, E. M. Kotsalis and P. Koumoutsakos, J. Comput. Phys., 2012, 231, 2677–2681.
18. S. Yasuda and R. Yamamoto, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2010, 81, 036308.
19. W. E, B. Engquist, X. Li, W. Ren and E. Vanden-Eijnden, Commun. Comput. Phys., 2007, 2, 367–450.
20. N. D. Petsev, L. G. Leal and M. S. Shell, J. Chem. Phys., 2015, 142, 044101.
21. D. A. Fedosov and G. E. K. B. Caswell, J. Chem. Phys., 2010, 132, 144103.
22. M. Praprotnik, L. Delle Site and K. Kremer, J. Chem. Phys., 2005, 123, 224106.
23. M. Praprotnik, L. Delle Site and K. Kremer, Annu. Rev. Phys. Chem., 2008, 59, 545–571.
24. M. Praprotnik, S. Poblete and K. Kremer, J. Stat. Phys., 2011, 145, 946–966.
25. S. Bevc, C. Junghans, K. Kremer and M. Praprotnik, New J. Phys., 2013, 15, 105007.
26. H. Wang, C. Hartmann, C. Schütte and L. Delle Site, Phys. Rev. X, 2013, 3, 011018.
27. J. Zavadlav, M. N. Melo, A. V. Cunha, A. H. de Vries, S. J. Marrink and M. Praprotnik, J. Chem. Theory Comput., 2014, 10, 2591–2598.
28. J. Zavadlav, M. N. Melo, S. J. Marrink and M. Praprotnik, J. Chem. Phys., 2014, 140, 054114.
29. A. Agarwal, H. Wang, C. Schütte and L. Delle Site, J. Chem. Phys., 2014, 141, 034102.
30. J. Zavadlav, M. N. Melo, S. J. Marrink and M. Praprotnik, J. Chem. Phys., 2015, 142, 244118.
31. J. Zavadlav, R. Podgornik and M. Praprotnik, J. Chem. Theory Comput., 2015, 11, 5035–5044.
32. H. Wang and A. Agarwal, Eur. Phys. J.: Spec. Top., 2015, 224, 2269–2287.
33. A. Agarwal, J. Zhu, C. Hartmann, H. Wang and L. D. Site, New J. Phys., 2015, 17, 083042.
34. K. Kreis, A. C. Fogarty, K. Kremer and R. Potestio, Eur. Phys. J.: Spec. Top., 2015, 224, 2289–2304.
35. D. Mukherji and K. Kremer, Macromolecules, 2013, 46, 9158–9163.
36. P. Szymczak and M. Cieplak, J. Chem. Phys., 2007, 127, 155106.
37. P. Szymczak and M. Cieplak, J. Phys.: Condens. Matter, 2011, 23, 033102.
38. B. Z. Dlugogorski, M. Grmela and P. J. Carreau, J. Non-Newtonian Fluid Mech., 1993, 48, 303–335.
39. P. J. Daivis and D. J. Evans, J. Chem. Phys., 1994, 100, 541–547.
40. Z. Xu, J. J. de Pablo and S. Kim, J. Chem. Phys., 1995, 102, 5836–5844.
41. J. T. Bosko, B. D. Todd and R. J. Sadus, J. Chem. Phys., 2005, 123, 034905.
42. T. C. Le, B. D. Todd, P. J. Daivis and A. Uhlherr, J. Chem. Phys., 2009, 131, 044902.
43. M. L. Matin, PhD thesis, Royal Melbourne Institute of Technology, Melbourne, 2001.
44. P. A. Thompson, G. S. Grest and M. O. Robbins, Phys. Rev. Lett., 1992, 68, 3448–3451.
45. M. Tschopp, J. L. Bouvard, D. Ward, D. Bammann and M. F. Horstemeyer, 2013, arXiv:1310.0728, 1–28.
46. E. G. Flekkoy, R. Delgado-Buscalioni and P. V. Coveney, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2005, 72, 026703.
47. D. Jou, J. Casas-Vázquez and M. Criado-Sancho, Thermodynamics of Fluids Under Flow, Springer, Netherlands, 2010.
48. P. Daivis, J. Non-Newtonian Fluid Mech., 2008, 152, 120–128.
49. J. Dealy and J. Wang, Melt Rheology and its Applications in the Plastics Industry, Springer, Netherlands, 2013.
50. J. Aho, Rheological Characterization of Polymer Melts in Shear and Extension: Measurement Reliability and Data for Practical Processing, Tampere University of Technology, 2011.
51. W. Friesenbichler, I. Duretek, J. Rajganesh and S. R. Kumar, Polimery, 2011, 56(1), 58–62.
52. F. Snijkers, K. Ratkanthwar, D. Vlassopoulos and N. Hadjichristidis, Macromolecules, 2013, 46, 5702–5713.
53. D. Vlassopoulos, G. Fytas, T. Pakula and J. Roovers, J. Phys.: Condens. Matter, 2001, 13, R855–R876.
54. C. N. Likos, Phys. Rep., 2001, 348, 267–439.
55. C. N. Likos, H. Löwen, M. Watzlawek, B. Abbas, O. Jucknischke, J. Allgaier and D. Richter, Phys. Rev. Lett., 1998, 80, 4450–4453.
56. S. M. T. Divoux, M. A. Fardin and S. Lerouge, 2015, arXiv:1503.04130, 1–16.
57. S.-Q. Wang, S. Ravindranath and P. E. Boukany, Macromolecules, 2011, 44, 183–190.
58. S. T. Milner, T. C. B. McLeish and A. E. Likhtman, J. Rheol., 2001, 45, 539–563.
59. C. Koch, A. Z. Panagiotopoulos, F. Lo Verso and C. N. Likos, Soft Matter, 2013, 9, 7424–7436.
60. D. Reith, M. Pütz and F. Müller-Plathe, J. Comput. Chem., 2003, 24, 1624–1636.
61. S. Bevc, C. Junghans and M. Praprotnik, J. Comput. Chem., 2015, 36, 467–477.
62. R. Delgado-Buscalioni and P. V. Coveney, J. Chem. Phys., 2003, 119, 978–987.
63. R. Potestio, S. Fritsch, P. Español, R. Delgado-Buscalioni, K. Kremer, R. Everaers and D. Donadio, Phys. Rev. Lett., 2013, 110, 108301.
64. R. Potestio, P. Español, R. Delgado-Buscalioni, R. Everaers, K. Kremer and D. Donadio, Phys. Rev. Lett., 2013, 111, 060601.
65. P. Español, R. Delgado-Buscalioni, R. Everaers, R. Potestio, D. Donadio and K. Kremer, J. Chem. Phys., 2015, 142, 064115.
66. R. Delgado-Buscalioni and G. De Fabritiis, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2007, 76, 036709.
67. E. A. Koopman and C. P. Lowe, J. Chem. Phys., 2006, 124, 204103.
68. T. Soddemann, B. Dünweg and K. Kremer, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2003, 68, 046702.
69. P. Espanol and P. Warren, Europhys. Lett., 1995, 30, 191.
70. C. Junghans, M. Praprotnik and K. Kremer, Soft Matter, 2008, 4, 156–161.
71. A. W. Lees and S. F. Edwards, J. Phys. C: Solid State Phys., 1972, 5, 1921.
72. D. J. Evans and G. P. Morriss, Phys. Rev. A: At., Mol., Opt. Phys., 1984, 30, 1528–1530.
73. A. J. Ladd, Mol. Phys., 1984, 53, 459–463.
74. J. Freire, in Branched Polymers II, ed. J. Roovers, Springer, Berlin Heidelberg, 1999, vol. 143, pp. 35–112.
75. M. SJ, Korea-Aust. Rheol. J., 2008, 20, 117–125.
76. F. Legrand and J.-M. Piau, J. Non-Newtonian Fluid Mech., 1998, 77, 123–150.
77. R. A. Mendelson, Trans. Soc. Rheol., 1965, 9, 53–63.
78. J. Casas-Vázquez and D. Jou, Rep. Prog. Phys., 2003, 66, 1937.
79. S. Bernardi, B. D. Todd and D. J. Searles, J. Chem. Phys., 2010, 132, 244706.
80. P. Padilla and S. Toxvaerd, J. Chem. Phys., 1996, 104, 5956–5963.
81. R. Khare, J. de Pablo and A. Yethiraj, J. Chem. Phys., 1997, 107, 2589–2596.
82. H. J. C. Berendsen, J. P. M. Postma, W. F. van Gunsteren, A. DiNola and J. R. Haak, J. Chem. Phys., 1984, 81, 3684–3690.
83. M. Christen and W. F. van Gunsteren, J. Chem. Phys., 2006, 124, 154106.
84. G. S. Grest, K. Kremer, S. T. Milner and T. A. Witten, Macromolecules, 1989, 22, 1904–1910.
85. G. S. Grest, K. Kremer and T. A. Witten, Macromolecules, 1987, 20, 1376–1383.
86. A. Eriksson, M. N. Jacobi, J. Nyström and K. Tunstrom, J. Chem. Phys., 2008, 129, 024106.
87. G. Strobl, The Physics of Polymers, Springer-Verlag, NY, 1997.
88. M. Kröger, W. Loose and S. Hess, J. Rheol., 1993, 37, 1057–1079.
89. M. Kröger, C. Luap and R. Muller, Macromolecules, 1997, 30, 526–539.
90. J.-P. Hansen and I. R. McDonald, in Theory of Simple Liquids, ed. J.-P. Hansen and I. R. McDonald, Academic Press, Burlington, 3rd edn, 2006, pp. 11–45.
91. R. Delgado-Buscalioni and A. Dejoan, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2008, 78, 046708.
92. R. Delgado-Buscalioni and P. V. Coveney, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2003, 67, 046704.
93. P. Allen and D. Tildesley, Computer simulation of liquids, Clarendon Press, 1987.
94. K. Yasuda, J. Text. Sci. Eng., 2006, 52, 171–173.
95. J.-P. Ryckaert, A. Bellemans, G. Ciccotti and G. V. Paolini, Phys. Rev. A: At., Mol., Opt. Phys., 1989, 39, 259–267.
96. P. Kindt and W. J. Briels, J. Chem. Phys., 2005, 123, 224903.
97. C. Baig and V. A. Harmandaris, Macromolecules, 2010, 43, 3156–3160.
98. S. Barsky, R. Delgado-Buscalioni and P. V. Coveney, J. Chem. Phys., 2004, 121, 2403–2411.
99. R. Delgado-Buscalioni, AIP Conf. Proceedings, 2007, 913, 114.
100. B. D. Todd, D. J. Evans and P. J. Daivis, Phys. Rev. E: Stat. Phys., Plasmas, Fluids, Relat. Interdiscip. Top., 1995, 52, 1627–1638.
101. R. L. C. Akkermans and G. Ciccotti, J. Phys. Chem. B, 2004, 108, 6866.
102. M. Ripoll, R. G. Winkler and G. Gompper, Phys. Rev. Lett., 2006, 96, 188302.
103. J. Brandao, E. Spieth and C. Lekakou, Polym. Eng. Sci., 1996, 36, 49–55.
104. X. Cheng, X. Xu, S. A. Rice, A. R. Dinner and I. Cohen, Proc. Natl. Acad. Sci. U. S. A., 2012, 109, 63–67.
105. M. Ripoll, R. Winkler and G. Gompper, Eur. Phys. J. E: Soft Matter Biol. Phys., 2007, 23, 349–354.
106. D. A. Fedosov, S. P. Singh, A. Chatterji, R. G. Winkler and G. Gompper, Soft Matter, 2012, 8, 4109–4120.
107. S. P. Singh, A. Chatterji, G. Gompper and R. G. Winkler, Macromolecules, 2013, 46, 8026–8036.
108. S. P. Singh, C.-C. Huang, E. Westphal, G. Gompper and R. G. Winkler, J. Chem. Phys., 2014, 141, 084901.
109. A. K. Tezel, J. P. Oberhauser, R. S. Graham, K. Jagannathan, T. C. B. McLeish and L. G. Leal, J. Rheol., 2009, 53, 1193–1214.
110. R. Delgado-Buscalioni, Phys. Rev. Lett., 2006, 96, 088303.
111. S. P. Singh, D. A. Fedosov, A. Chatterji, R. G. Winkler and G. Gompper, J. Phys.: Condens. Matter, 2012, 24, 464103.
112. J. Li and J. Lee, Acta Mech., 2014, 225, 1223–1233.

---

Footnotes

† An interesting remark is that the recoverable shear compliance⁸⁷ obtained from J_e = τ_rot/η₀ results to be J_e ≃ 110 independent of the friction kernel.

‡ Traditionally, following the experiment standpoint, rheological magnitudes are expressed in terms of the pressure on the system (here given by −P). For instance, the normal stress difference is then N₁ = −(P₁₁ − P₂₂) and this explains the opposite sign used in our (and other) simulation studies.

§ This can be proved from the equation for the time dependence of the covariance 〈dv_i/dt(t)dv_i/dt(t′)〉 of the DPD Langevin's equation and can be easily checked (and the cooling rate α measured) upon observation of an exponential convergence of T towards T₀ after an instantaneous change (increase or decrease) of T₀.

---