Mesoscale modeling of shear-thinning polymer solutions I. S. Santos de Oliveira, B. W. Fitzgerald, W. K. den Otter, and W. J. Briels Citation: The Journal of Chemical Physics 140, 104903 (2014); doi: 10.1063/1.4867787 View online: http://dx.doi.org/10.1063/1.4867787 View Table of Contents: http://scitation.aip.org/content/aip/journal/jcp/140/10?ver=pdfcov Published by the AIP Publishing Articles you may be interested in Extensional rheology of shear-thickening fumed silica nanoparticles dispersed in an aqueous polyethylene oxide solution J. Rheol. 58, 411 (2014); 10.1122/1.4864620 A new model for dilute polymer solutions in flows with strong extensional components J. Rheol. 46, 1057 (2002); 10.1122/1.1501963 The effects of viscoelasticity on the transient motion of a sphere in a shear-thinning fluid J. Rheol. 41, 103 (1997); 10.1122/1.550803 Nonisothermal Elongational Flow of Polymeric Fluids According to the Leonov Model J. Rheol. 28, 581 (1984); 10.1122/1.549763 Official Nomenclature for Material Functions Describing the Response of a Viscoelastic Fluid to Various Shearing and Extensional Deformations J. Rheol. 28, 181 (1984); 10.1122/1.549739

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 128.82.252.58 On: Sat, 20 Dec 2014 01:19:10

THE JOURNAL OF CHEMICAL PHYSICS 140, 104903 (2014)

Mesoscale modeling of shear-thinning polymer solutions I. S. Santos de Oliveira,a) B. W. Fitzgerald, W. K. den Otter, and W. J. Brielsb) Computational Biophysics, University of Twente, P.O. Box 217, 7500 AE, Enschede, The Netherlands

(Received 19 December 2013; accepted 25 February 2014; published online 12 March 2014) We simulate the linear and nonlinear rheology of two different viscoelastic polymer solutions, a polyisobutylene solution in pristane and an aqueous solution of hydroxypropylcellulose, using a highly coarse-grained approach known as Responsive Particle Dynamics (RaPiD) model. In RaPiD, each polymer has originally been depicted as a spherical particle with the effects of the eliminated degrees of freedom accounted for by an appropriate free energy and transient pairwise forces. Motivated by the inability of this spherical particle representation to entirely capture the nonlinear rheology of both fluids, we extended the RaPiD model by introducing a deformable particle capable of elongation. A Finite-Extensible Non-Linear Elastic potential provides a free energy penalty for particle elongation. Upon disentangling, this deformability allows more time for particles to re-entangle with neighbouring particles. We show this process to be integral towards recovering the experimental nonlinear rheology, obtaining excellent agreement. We show that the nonlinear rheology is crucially dependent upon the maximum elongation and less so on the elasticity of the particles. In addition, the description of the linear rheology has been retained in the process. © 2014 AIP Publishing LLC. [http://dx.doi.org/10.1063/1.4867787] I. INTRODUCTION

Coarse-grained approaches towards the simulation of polymer liquids allow access to large time and length scales that would otherwise be prohibited with atomistic algorithms. This approach leads to both a simplified polymer representation and a simpler description of topological constraints or entanglements. It is these constraints that control the intermediate and long-term behaviours and are most relevant in highly concentrated systems such as melts. Varying degrees of coarse graining maybe considered. Depending on the simulation requirements, a polymer entity can be approximated as a series of coarse-grained blobs or beads each representing a group of constitutive monomers,1–3 while a highly coarsegrained representation views each entity as a point particle.4, 5 One such highly coarse-grained model is the Responsive Particle Dynamics (RaPiD) algorithm5–8 which was originally introduced to simulate polymer melts. In a previous study, we adapted the RaPiD method to simulate a variety of polymer solutions in the semidilute regime.9 Simulation results were found to be in qualitative agreement with experimental linear and nonlinear measurements for both a polyisobutylene (PIB) solution and a wormlike micellar solution.10 The shear thinning behaviour for the polymers was well described at low shear rates, but for higher shear rates the flow curves were slightly steeper than for a real solution. The emergence of shear thinning in solutions is most likely due to the fact that recently disentangled polymers have insufficient time to reentangle with new neighbours. Therefore, the reason for the faster shear thinning effect observed in our previous RaPiD simulations is the absence of polymer deformation and elona) Present address: Instituto de Física, Universidade Federal de Uberlândia,

C.P. 593, 38400-902 Uberlândia, Brazil.

b) Electronic address: [email protected].

0021-9606/2014/140(10)/104903/11/$30.00

gation that would allow particles more time to re-entangle with particles in their neighbourhood after disentangling. In this work, we extend the RaPiD algorithm to allow polymers to deform and elongate when subject to shear flow, thus providing particles with more time to entangle with other particles in their vicinity. To facilitate this we depict each RaPiD entity as a deformable particle with a FENE (Finite-Extensible Non-Linear Elastic) potential providing the penalty for deformation. Entanglements between particles are updated to account for overlap between the neighbouring components of deformable particles. This paper is organized as follows. The simulation method, including the description of the RaPiD deformable particles, is outlined in Sec. II. In Sec. III, we provide details of the systems to be simulated and the methods towards calculating linear and nonlinear measures. The effect of the RaPiD parameter set on rheological measurements for one polymer solution is described in Sec. IV. In Sec. V, we use the observations from Sec. IV to obtain a quantitative agreement of the model fluid linear rheology with experiments and to predict the nonlinear rheology for two different polymer systems. Finally, in Sec. VI we draw conclusions on this work.

II. THE RAPID SIMULATION MODEL

In this work, we study the possibility to simulate the rheological properties of a large class of complex fluids. Our general aim is to develop a simulation model that can be used to simulate large amounts of these fluids, possibly in complicated geometries and complicated flows. In order to speed up the simulations as much as possible, the description of the structural properties of the molecules constituting the fluid must be limited to a bare minimum necessary to obtain the correct flow behaviour. In the RaPiD model that we use

140, 104903-1

© 2014 AIP Publishing LLC

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 128.82.252.58 On: Sat, 20 Dec 2014 01:19:10

104903-2

Santos de Oliveira et al.

J. Chem. Phys. 140, 104903 (2014)

here, each molecule is represented by its center of mass position ri and a deformation vector di . Particles are subject to two forces; the average forces depending on the configuration (r, d) of the RaPiD particles, and transient forces describing the correlation between the dynamics of the eliminated small scale degrees of freedom and the dynamics of the retained degrees of freedom, i.e., the positions and deformations of the RaPiD particles. In Subsections II A and II B, we will discuss the corresponding contributions to the potential energy. In Subsection II C, we discuss the Brownian dynamics propagator used in the model.

A. Potential of mean force

The average forces on the RaPiD particles are governed by the so called potential of mean force defined as  c (r, d) = −kB T ln

 exp{−V (r, d, q)/kB T }dq ,

(1)

where V (r, d, q) is the potential energy of the small scale system with r and d the RaPiD degrees of freedom, or coarse variables, and q is the internal degrees of freedom of all particles and possibly solvent degrees of freedom; kB is Boltzmann’s constant, and T is the temperature. A simulation of the coarse variables using the potential of mean forces as its energy will produce the exact structural and thermodynamic properties of the system. In practice, of course, the potential of mean force can only approximately be represented in a form depending very much on the characteristics of the system to be simulated.11 In the present application of the RaPiD model, we extend our previous work on polymer solutions9 by allowing the particles to deform. We therefore propose the following potential of mean force: c (r, d) = f h (r) + f ene (d).

(2)

The first term represents the free energy of the eliminated degrees of freedom in the equilibrium state, i.e., when all deformation vectors are zero. The second term provides the free energy penalty for deforming the polymers. We now briefly describe in the remaining part of this subsection the main ingredients of both contributions to the potential of mean force. For solutions, interactions between solvent and solute molecules are important and must be included in the simulation model. Instead of explicitly simulating the solvent particles, we use the Flory-Huggins (FH) theory12, 13 to incorporate interactions between polymers and solvent molecules in the potential of mean force of the polymers, and thus obtain an implicit solvent model. In the spirit of Pagonabarraga and Frenkel7, 14, 15 we write the potential of mean force as a sum of contributions from each polymer molecule.9 The actual contribution by polymer i is equal to the configurational free energy of a polymer in a homogeneous solution with polymer volume fraction equal to φ i . The latter is the local volume fraction in the simulated solution at the position of polymer i (see below). The first contribution to the potential of mean

force c is then f h (r) =

Np 

ap (φi )

i=1

= pkB T

Np   1 − φi i=1

φi

 ln(1 − φi ) − χ φi ,

(3)

where p is the number of Kuhn segments in the chain and Np is the number of polymer chains. The interaction between the polymer chains and the solvent molecules is described by the Flory-Huggins parameter χ . Details of the derivation of this expression are provided in the Appendix of Ref. 9. The reference states for measuring both energy and entropy are pure liquids. Notice that the free energy contributed by translational degrees of freedom of the coarse variables has been excluded from the potential of mean force since it will be produced by the simulation. Stated differently chain centers are considered to be fixed in the derivation of Eq. (3). χ can be used to control the quality of the solvent in relation to the polymers. The local chain volume fraction φ i around particle i is computed as φi (r) =

Np 1 

ρmax

ω(rij ),

(4)

j =1

where ρ max defines the maximum number density of polymers the system is allowed to reach. The weight function ω(rij ) represents the relative contribution of the monomers of polymer j to the density experienced by the monomers of particle i, hence ω(rij ) must be a function that smoothly decays with increasing rij until it becomes zero at some cutoff value of about twice the radius of gyration of the polymers. Depending on the quality of the solvent, besides changing χ , it may also be necessary to change ω(rij ). A random distribution of the centers of mass of all polymers should turn the sum in Eq. (4) into the overall number density of chains, hence ω(rij ) must integrate to unity. For further details we refer to Ref. 9. The penalty for deforming the polymers will be described by a sum of FENE potentials for the individual polymers   2 Np 1 2 di ln 1 − , (5) f ene (d) = − kd0 2 d 0 i=1 √ where k is the spring constant, di = di · di with di being the deformation vector of particle i, and d0 its maximum value. Similar contributions to the potential of mean force have been used by Padding et al.16, 17 in the case of pressure sensitive adhesives. B. Transient forces

We now discuss the second type of forces felt by the RaPiD particles. In order to clearly describe the concepts we shall initially ignore the possibility for the particles to deform. At the end of this subsection, we then mention the changes to be made in the case of deformable particles.

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 128.82.252.58 On: Sat, 20 Dec 2014 01:19:10

104903-3

Santos de Oliveira et al.

It is well known that in polymer systems the chain connectivity in combination with strong repulsive forces gives rise to “topological constraints.” As an example of such forces, consider two well mixed overlapping linear polymers. Topological constraints will temporarily strongly resist displacements of the centers of mass of the two polymers, while in the equilibrium states before and after the displacement the free energies of the internal degrees may not be that different. On the other hand, it will take some time before two entangled polymers mix after being brought together. By taking minus the gradient of the potential of mean force, Eq. (1), with respect to one of the coarse variables ri one obtains an expression that averages the force on the chosen variable over all internal degrees of freedom q. Topological constraints, however, cause the system to make extremely large detours in order to achieve seemingly small changes in its configuration. In the potential of mean force all configurations are taken into account, while in reality only a subset of these is sampled until the coarse configuration has changed substantially. It is obvious therefore that the forces that drive the dynamics of the coarse system are not equal to those derived from the potential of mean force. For simplicity we will refer to topological constraints by the name entanglements. The reader should not confuse “the number of entanglements” or “entanglement number,” used below, with Z = M/Me , which in reptation theory indicates that the mass of the polymer is Z times its entanglement mass.18 The effect of entanglements is included into the model by adding to the list of dynamical variables the set n = {nij |rij ≤ rc }, i.e., one dynamical variable for each pair of polymers whose distance is less than some cutoff. nij is called the number of entanglements between particles i and j. It keeps track of the thermodynamic state of the eliminated coordinates, which for the time being include all deformations. For each coarse configuration, nij is assumed to have an equilibrium value depending only on the distance between the two polymers, and denoted by n0 (rij ). For Gaussian chains, n0 (rij ) may be chosen to be proportional to the number of monomercontacts, or the overlap between two chains at a distance rij . This overlap between Gaussian chains is itself approximately a Gaussian function of the distance between the centers of mass of the two polymers.19 However, to avoid zero forces at short distances we represent n0 (rij ) by a quadratically decaying function, as follows: ⎧ 2 ⎨ rij − 1 : rij ≤ rc , n0 (rij ) = (6) ⎩ rc 0 : rij > rc , where rc is the cutoff distance. When all entanglement numbers are equal to their equilibrium values, the eliminated degrees of freedom are in equilibrium and the forces are equal to those derived from the potential of mean force. Deviations of entanglement numbers from their equilibrium values give rise to additional forces. These forces are modeled by adding to the potential of mean force a transient potential t (r, n). To lowest order, the general form of the transient potential is a quadratic function in the deviations of the entanglement numbers from their equilibrium values,5 representing the tendency of the system to always move towards its equilibrium

J. Chem. Phys. 140, 104903 (2014)

state. So, t (r, n) =

1 s α (nij − n0 (rij ))2 . 2 i,j

(7)

The parameter α s controls the strength of the fluctuations in the entanglement numbers deviations from their equilibrium values. Forces due to t are called transient forces. If we were to freeze the values of the coarse variables, i.e., we would freeze the coarse configuration of the system, the entanglement numbers would gradually relax towards their equilibrium values and the corresponding forces would relax to zero. As can be seen from Eq. (6), since n0 (r) is a monotonic function of r, the entanglement numbers drive the system towards configurations with interparticle distances such that all n0 (ri, j ) are equal to the instantaneous values of the ni, j . The effect of the transient forces is to severely slow down the dynamics of the particles, mimicking the effect of entanglements between long flexible chains. In a nutshell, one might say that the centers of mass coordinates alone do not constitute the full set of slow variables that determine the dynamics of the system, and that the remaining slow variables have been lumped together into one thermodynamic degree of freedom for every interacting pair of particles. From this point of view, making the model less coarse is to transfer a number of slow variables from the thermodynamic part, where they are incorporated implicitly, to the dynamic part where they are treated explicitly. This we will now do by making the particles deformable. The simplest way to make the particles deformable is to consider each polymer to consist of two components which can entangle with similar components of surrounding polymers. We therefore introduce entanglement centers 1 ri,1 = ri − di , 2

(8)

1 ri,2 = ri + di , (9) 2 where ri is the deformable particle centre of mass and di is the distance between the entanglement centres located at ri,1 and ri,2 . We extend the transient potential to read  2 1  1 nim,j n − n0 (rim,j n ) . (10) t (r, d, n) = α 2 i,j m,n 4 Here n and m run through the values 1 and 2, nim, jn is the number of entanglements associated with the entanglement centers i, m and j, n, and rim, jn is the distance between these two centers. Notice that we retain the same definition of n0 , as in Eq. (7), i.e., the one given in Eq. (6). The factor of one fourth in front of n0 (rim, jn ) is taken for cosmetic reasons (see the Appendix). C. Brownian dynamics

Having defined the interaction potentials between the particles, we proceed to define the propagator which maps the configuration at time t to that at time t + dt. Since the time evolution of soft matter systems is highly overdamped

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 128.82.252.58 On: Sat, 20 Dec 2014 01:19:10

104903-4

Santos de Oliveira et al.

J. Chem. Phys. 140, 104903 (2014)

we adopt the Brownian dynamics propagator where the position of particle i is updated during the simulation according to Refs. 5 and 6,

2kB T dt 1 ∂ dri = − (c + t )dt + i ξi ∂ri ξi + V (yi )ˆex dt.

(11)

The first term on the right hand side contains minus the gradient with respect to ri of the sum of the conservative and transient potential, resulting in the force applied to particle i. The force due to the transient potential is a sum of pair terms, each proportional to the derivative of n0 (rim, jn ) for some value of j, m, and n. In order to prevent such pair-forces becoming zero when the corresponding distance is zero, we have ensured that the derivative of n0 (r) for r = 0 is non-zero as expressed in Eq. (6). The force on particle i is turned into a displacement over time dt by dividing by the friction ξ i and multiplying by the time step dt. The friction felt by a particle is assumed to depend linearly on the degree of entanglement with neighbouring chains:   1 (12) |nim,j n | n0 (rim,j n ), ξi = ξ s + ξe 4 j =i n,m where ξ e is the friction per entanglement and ξ s is the background friction by the solvent. The somewhat peculiar way of writing the linear dependence on the number of entanglements by taking the product of |nim, jn | and n0 (rim, jn ) and next taking the square root is to ensure that the friction between two particles becomes zero beyond the cutoff distance. The second term in the Brownian propagator describes the random displacements characteristic for Brownian motion; i is a time-dependent Markovian random vector, with unit variance and zero mean. The last term is included when a shear flow along the x-direction is applied, with an overall velocity gradient γ˙ in the y-direction imposed by Lees-Edwards sliding boundary conditions.20 Similarly, the equation of motion for the deformation vector di connecting the polymer components is given by 1 ∂ (f ene + t )dt ξid ∂di

2kB T dt + i + [V (yi2 ) − V (yi1 )]ˆex dt. (13) ξid

ddi = −

The first term on the right hand side now describes the derivative of the FENE and transient potentials with respect to the deformation vector di . In the Appendix, we will argue that the friction on the deformation vector is equal to one fourth of the friction on the centre of mass: ξid = ξi /4. The next term in Eq. (13) is again the Brownian propagator describing random motion. The last term describes the motion due to the applied shear flow, where V (yi1 ) and V (yi2 ) are the streaming velocity at the two extremes of the polymer. For any given configuration of the coarse variables, the entanglement numbers relax towards their corresponding equilibrium values. This relaxation is assumed to be governed

by a simple linear equation augmented with a random contribution obeying the fluctuation-dissipation theorem:5, 6   1 1 dnim,j n = n0 (rim,j n ) − nim,j n dt τ 4  2kB T dt + im,j n . (14) ατ The first term gives rise to an exponential decay of nim, jn towards n0 (rim, jn ) with characteristic time τ . In most polymer systems, a single relaxation time is not enough to describe the correct dynamics of the system. In the case of linear polymers, dynamical processes occur over a wide spectrum of characteristic times, either of the Rouse or the reptation type. In other cases, a wide range of characteristic times exists as experimental polymer systems are polydisperse. From a different point of view, one may expect that Gaussian chains which hardly overlap may disentangle faster than those which overlap severely. In the former case, the number of entanglements is not only smaller than in the latter case, but the entanglements are less dense, and therefore disentangle more easily. This effect is included in the model by refining the definition of the relaxation time to be a function of the distance between two neighbouring particles as given by the expression   r im,j n , (15) τim,j n = τ0 exp − λ where τ 0 is a time constant and λ defines the decay length of the relaxation time. The second term on the right hand side of Eq. (14) describes random fluctuations in nim, jn , with im, jn being a time-dependent random Markovian scalar with zero mean, unit average, and without correlations across particle pairs. III. SYSTEMS AND METHODS

The RaPiD model, as presented in Sec. II, has a number of parameters and functions that allow the simulation of a large class of materials. Some of the parameters refer to simple experimental conditions, such as temperature and volume fraction of polymer, or molecular properties, such as the radius of gyration Rg , the number of Kuhn segments in the polymer, the Flory-Huggins parameter and the solvent friction ξ s . These parameters will be called the “set parameters.” A second class of parameters, referred to as the “simulation parameters,” which composed of α, τ , λ, and ξ e is defined with reference to the transient potential outlined in Sec. II. Together with n0 (rij ) they describe the coupling between the dynamics of the eliminated and retained coordinates. In the transient potential t , the variable α governs the strength of the fluctuations of the entanglement numbers around their equilibrium values and at the same time the strength of the transient forces. The characteristic relaxation time τ (Eq. (14)) can, in principle, be estimated from the crossing frequency of the storage G (ω) and loss G (ω) moduli measured in experiments. However, since in most polymer systems a range of frequencies is observed, which we introduce by using Eq. (15), the crossing frequency only provides a rough estimate of the dominant characteristic time τ d . Assuming that the radial distribution function between polymers is very close

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 128.82.252.58 On: Sat, 20 Dec 2014 01:19:10

104903-5

Santos de Oliveira et al.

J. Chem. Phys. 140, 104903 (2014)

TABLE I. List of experimental (set) parameters and simulation parameters used to model the linear rheology for a PIB (polyisobutylene) and a HPC (hydroxypropylcellulose) polymer solutions. Value Set parameters Temperature Radius of gyration Density Average volume fraction Flory-Huggins parameter Number of monomers Solvent friction Simulation parameters Entanglement number deviation Entanglement friction Maximum τ Decay length of τ Spring constant Maximum spring elongation

Symbol

PIB

HPC

Unit

T Rg ρp φ χ p ξs

300 40 3.5 0.11 0.5 2700 2.45 × 10−9

300 15 4.3 0.3 0.5 70 1.84 × 10−10

K nm particles/Rg3 ... ... mon./polymer kg/s

α ξe τ0 λ k d0

18 1 × 10−5 4 0.30 60 0.3

1.75 2 × 10−5 3 0.41 150 0.3

kB T kg/s s Rg kB T /Rg2 Rg

to unity, one may expect that τ d in the simulation model is the one that corresponds to a separation r that makes r2 exp {−r/λ} a maximum, i.e., τ d = τ 0 e−2 . The parameters in this second class, together with k and d0 which are associated with the FENE potential, will be adjusted to describe the experimental data as well as possible.

from which the storage and loss moduli can be computed as

A. Systems and experimental data

We have chosen two rather different polymer solutions to illustrate the viability of our model. The first is a high molecular weight PIB (polyisobutylene) solution in Pristane. The linear and nonlinear rheological properties of this system have been published in Ref. 10. From the materials properties and conditions of the experiment, we define the set parameters listed in Table I. The second system is an aqueous solution of hydroxypropylcellulose (HPC). Rheological measurements for this fluid are reported in Refs. 21 and 22. The set parameters for this system which are defined with reference to experimental conditions and material parameters are listed in Table I. Experimental data in the linear regime are the storage and loss moduli G (ω) and G (ω), respectively. They are related to the shear relaxation modulus G(t) by  ∞ sin(ωt)G(t)dt, (16) G (ω) = ω 0





n  Gi τi2 ω2 , 1 + τi2 ω2 i=1

(19)

G (ω) =

n  Gi τi ω , 1 + τi2 ω2 i=1

(20)

where Gi and τ i were obtained by adjusting them until the experimental data were reproduced. In both cases, it turned out that good representations of the data could be obtained with n = 4; the corresponding results are given in Table II. We stress that the shear relaxation moduli obtained with the above fit procedure are just guides to the eye allowing a quick comparison with our simulated results. The quality of our final simulations must be judged on their capability to reproduce the experimental storage and loss moduli.

B. Methods

All simulations are performed in cubic boxes with V = 216 Rg3 , containing 756 coarse-grained particles. TABLE II. Coefficients used to fit the experimental G (ω) and G (ω) data with the Maxwell model to obtain G(t), for a PIB and HPC solutions.

∞

G (ω) = ω

G (ω) =

cos(ωt)G(t)dt.

(17)

Polymer solution

i

Gi (Pa)

τ i (s)

PIB

1 2 3 4 1 2 3 4

848.36 314.24 65.80 2.65 5809.89 1958.30 653.68 42.79

0.008 0.079 0.541 3.707 0.011 0.033 0.205 1.744

0

Since computer simulations give quick access to the shear relaxation modulus we have inverted these equations in order to obtain G(t). To this end we have assumed that the shear relaxation modulus may be represented by a sum of Maxwell modes as G(t) =

n 

Gi exp(−t/τi )

(18)

HPC

i=1

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 128.82.252.58 On: Sat, 20 Dec 2014 01:19:10

104903-6

Santos de Oliveira et al.

J. Chem. Phys. 140, 104903 (2014)

We calculate the shear relaxation modulus according to (21)

ξe = 10

-6

0 2

1

0 2

ξe = 10

-8

(22) 1

0

-4

α = 1 kB T ξe = 10

2 1

(23)

-6

ξe = 10

1

-8

0

ξe = 10

2 1 0

α = 10 kBT ξe = 10

-4

3 2 1 0

-6

3

ξe = 10

In this section, we investigate the characteristic features of the shear relaxation modulus G(t) as a function of the RaPiD parameters for the PIB polymer solution. We restrict ourselves to the simplified non-deformable (spherical) particle method outlined in Sec. II and presented in more detail in Ref. 9, since we have found that the RaPiD parameters for the spherical particle model can be easily translated to those for the deformable particle model, using the relations outlined in the Appendix, while basically retaining the linear rheology. In Fig. 1, we plot G(t) curves for different values of α, ξ e , τ 0 , and λ subject to the “set parameters” for the PIB polymer solution listed in Table I. For constant α and at small times all G(t) curves have the same maximum value. In some cases with λ = 0.2Rg (black curves), this is not observed over the time scales shown but at lower time scales (not shown here) the curves reach this value. As α increases, the maximum value in G(t) increases, indicating an increasing contribution of transient forces to the stress tensor in Eq. (22). Two pronounced relaxation zones are evident, particularly for small τ 0 where G(t) appears to be a sum of two independent contributions. For constant values of ξ e the slowest contribution is approximately independent of the transient potential suggesting it develops from the conservative force. The faster contribution, which depends on all parameters, increases with α and is insensitive to the friction. The actual relaxation of G(t) depends strongly on τ 0 and λ. When λ = ∞ such that τ = τ 0 the entanglement contribution to G(t)

2

2 1 0 3

-8

IV. LINEAR RHEOLOGY DEPENDENCY ON RAPID PARAMETERS

0

log (G (t)) [Pa]

where Sxy is the xy-component of the stress tensor, and γ˙ is the applied shear rate.

log (G (t)) [Pa]

Sxy (γ˙ ) , γ˙

ξe = 10

1

The sums only run through pairs, with each pair counted once. The first line includes forces arising from fene and t , the second those from fh . Non-equilibrium runs were performed using LeesEdwards boundary conditions with flows in the x-direction and gradients in the y-direction. Viscosities were calculated according to η(γ˙ ) =

-4

2

where V is the volume of the simulation box and Sxy is the xy component of the stress tensor. The angled brackets indicate that Sxy at some time t0 , multiplied by its value at time t0 + t, is averaged over all times t0 during an equilibrium run. The stress tensor is calculated according to the expression 1  f ene,t (rim − rj n )Fim,j n S=− V i,j m,n 1  fh − (ri − rj )Fi,j . V i,j

τ0 = 20 s

ξe = 10

V Sxy (t)Sxy (0), kB T

log (G (t)) [Pa]

G(t) =

α = 0.1 kBT τ0 = 2 s

τ0 = 0.2 s

2 1 0 -4

-2

0

2

-4

-2

0

2

-4

-2

0

2

log (t)

FIG. 1. Stress autocorrelation function G(t) calculated for a set of simulation parameters, as a function of time. On upper panel α = 0.1 kB T, in the middle α = 1 kB T, and at the lower panel α = 10 kB T. The black curves are calculated for λ = 0.2 Rg , blue ones for λ = 0.4 Rg , and red curves are for constant τ = τ 0 . In each panel from top to bottom ξ e = 10−4 kg s−1 , 10−6 kg s−1 , and 10−8 kg s−1 , while from left to right τ 0 = 0.2 s, 2 s, and 20 s.

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 128.82.252.58 On: Sat, 20 Dec 2014 01:19:10

Santos de Oliveira et al.

J. Chem. Phys. 140, 104903 (2014)

is approximately an exponentially decaying function with a decay time determined by τ 0 . We observe that the decay of G(t) becomes less abrupt with decreasing values of λ indicating that the spectrum of decay times is tuned by λ. Finally, we consider the effect of ξ e on G(t). The friction variable ξ e controls the time taken for the system to relax such that decreasing ξ e leads to a reduction in the time taken, i.e., the curves shift to the left. As the friction decreases particles undergo larger displacements, thus decreasing the time to decorrelate stress-values. Since the friction does not change the maximum values of the stress tensor, the G(t) values at small times remain the same.

(a)

3

10

2

[Pa]

104903-7

Gexp(t) Gs(t) Gd(t) Gc(t) Gt(t) Gf(t)

10

1

10

0

10 -4 10

-3

-2

10

10

-1

0

10

t [s]

10

(b)

3

10

V. TWO EXPERIMENTAL POLYMER SOLUTIONS 2

10

[Pa]

In this section, we will demonstrate the effectiveness of the RaPiD simulation model to describe the linear and nonlinear rheology of viscoelastic polymer solutions.

(ω)

1

10

0

10

A. PIB polymer solution with deformable particles -1

η0 η(γ˙ ) = , 1 + (τcr γ˙ )

(24)

10 -2 10

-1

10

0

10

1

2

10

3

10

ω [rad/s]

10

FIG. 2. Linear rheology of a PIB polymer solution. (a) G(t), (b) G (ω), and G (ω). Red lines are experimental results by Snijkers et al.,10 where the experimental G(t) is obtained by applying the Maxwell model. The black circles are simulation results using the spherical particle representation and open blue squares using deformable particles to represent the RaPiD particles. In (a), the dotted-dashed line is the shear relaxation including only conservative forces Gc (t), the dashed line accounts for transient forces Gt (t), and the dotted line the FENE forces.

where η0 is the zero shear viscosity, τ cr is a constant and the exponent  represents the shear-thinning slope. Fig. 4 shows  as a function of the maximum extension d0 . For d0 = 0.1 Rg the slope obtained is approximately that for a spherical particle while it decays with increasing d0 to a value of about 0.78 for d0 = 0.7 Rg . As d0 increases from 0.1 to 0.7 Rg , η0 increases by 17%, while τ cr varies around 0.215 s, with a variation of only ±0.025 s. The effect of the spring constant on 2

10

. η (γ) [Pa.s]

We consider the PIB polymer solution discussed in Ref. 10. First, we estimate the “simulation parameters” using the spherical particle model. The experimental shear relaxation modulus for this system is shown in Fig. 2(a). To recover this curve we have followed the trends observed in Sec. IV where we identified the effect of each parameter on G(t). First, we tuned α to obtain G(t) at small times. To approximate the relaxation in G(t) both τ 0 and λ were varied concurrently without giving much attention to shifts along the time axis. Finally, we adjust ξ e to get the correct position on the time axis. For the PIB polymer solution the RaPiD parameters found are α = 6 kB T, τ 0 = 4 s, λ = 0.3 Rg , and ξ e = 6 × 10−6 kg/s. The calculated shear relaxation modulus, Gs (t), is shown in Fig. 2(a), and the corresponding storage and loss moduli are presented in Fig. 2(b). Having reproduced the linear rheology remarkably well, we now consider the nonlinear rheological response by calculating the viscosity (Eq. (23)) as a function of the shear rate. As can be seen in Fig. 3, the simulated shear viscosity (closed circles) shows good agreement for low shear rates and the shear-thinning behaviour is qualitatively recovered over the range of shear rates studied. For higher shear rates however we note that the decay in η(γ˙ ) is much steeper than with the real fluid. It is our aim to ameliorate on this by applying the deformable particle model. Using the relations in the Appendix, we translate the transient parameters for the spherical particle model to those for the deformable particle model. In this model, we have two further parameters; the maximum extension or elongation of the polymers d0 and the spring constant k. To test the effect of d0 on the rheology we set k = 60 kB T /Rg2 and have simulated the shear viscosity at a range of shear rates for many d0 values. These shear viscosities curves were fitted with the Cross model23

Experiment Spherical particles Deformable particles

1

10

0

10

-1

10

0

10

1

. -1 γ [s ]

10

2

10

FIG. 3. Comparison between the simulated and experimental shear viscosity for a PIB polymer solution. Simulation results for the RaPiD particles represented as spherical particles (circles) and deformable particles (squares), while the solid red line represents experimental measurements.10

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 128.82.252.58 On: Sat, 20 Dec 2014 01:19:10

Santos de Oliveira et al.

104903-8

J. Chem. Phys. 140, 104903 (2014)

0.1

1.02 Spherical particles - 0.12 ε (d0) = 0.76 d 0 Deformable particles

0.08

P (d)

ε (d0) [ - ]

0.96

0.9

0.06

. γ=0 . -1 γ=1s . -1 γ = 10 s . -1 γ = 32 s . -1 γ = 100 s

0.04

0.84 0.02 0.78 0.1

0.2

0.3

0.4

d0 [Rg]

0.5

0.6

0.7

FIG. 4. Shear-thinning slope , defined in Eq. (24), as a function of the maximum elongation length of the FENE springs (d0 ). Each circle is a result from the fit of the simulated η(γ˙ ) for shear rates ranging from 1 to 100 s−1 . The dashed line is a fit to the resulting (d0 ), and the red dotted line represents the slope for systems including spherical particles only.

the shear viscosity was investigated by fixing the maximum spring deformation to d0 = 0.25 Rg while k was varied from 40 to 100 kB T /Rg2 . For the range of k values studied,  was found to be constant at 0.89, η0 decreased by only 4% from its maximum value at k = 40 kB T /Rg2 to a minimum at k = 100 kB T /Rg2 and τ cr remained approximately constant at 0.215 s. We conclude that the results are relatively insensitive to variations in the spring constant k, and that the maximum deformation d0 determines the final slope in the flow curves. Given the experimental value of  we read from Fig. 4 that d0 must be 0.3 Rg . In order to compensate for the effects of this choice on η0 and the overall time scale of G(t), we had to fine tune the original estimate of the transient parameters. The resulting parameters are given in Table I and the linear rheology obtained with these parameters is shown in Fig. 2 where the subscript “d” refers to the deformable particle model. Also shown in Fig. 2 are the contributions to G(t) from each of the individual forces experienced by the deformable particles. Cross correlation of the various stress contributions turned out to be very small. At low to intermediate times the shear relaxation modulus is mostly described by the transient force while the FENE contribution, when added to the transient contribution, adjusts the total G(t) to the experimental value. The Flory-Huggins conservative potential is responsible for the description of the curve at longer times. In Fig. 3 it is seen that, with the present model, the shearthinning behaviour of the system is fully recovered. The increase of the shear viscosity at higher shear rates when compared to the spherical particles representation is attributed to the deformation of the FENE springs. In Fig. 5, we plot the probability distribution of the deformable particle elongation for various shear rates. At γ˙ = 0 the correct equilibrium distribution is recovered. For increasing γ˙ the peak shifts to larger d indicating increased elongation. Due to this elongation, the entangling of two neighbouring particles lasts for longer and this prolongs the effect of the transient potential. It is tempting to compare Fig. 5 with information available in the literature on the deformation of polymers in shear flow, in particular Fig. 1 in the work of Baig et al.24 Direct comparison however is difficult since our systems are poly-

0 0

0.1

0.2

0.3

d [Rg]

FIG. 5. Probability distribution of the deformable particle elongation (d) for a range of shear rates with d0 = 0.3 Rg and k = 60 kB T /Rg2 . The black dashed line is obtained with the Boltzmann factor with the FENE potential.

mer solutions of rather low concentration (although obviously above c*), while those in Ref. 24 are polymer melts in which stretching of polymers will be much more severe. Moreover, it is difficult to guarantee the same definition of Weissenberg numbers in both cases. Since our overall relaxation time is in the order of a few tenths of a second, our Weissenberg numbers range from 0 to about 50. In this range, our results are roughly similar to those in Ref. 24. At the large Weissenberg numbers in Ref. 24 chains are severely stretched and the radius of gyration increases appreciably. We expect however that this effect will be much less pronounced in polymer solutions of concentrations as low as ours. As a consequence, it must be expected that d0 as it has been introduced here depends on the concentration of the polymers. It is conceivable that a less restrictive potential than the FENE potential will be more generic. We also computed the orientational order parameters for the deformable particles at various shear rates. The orientational ordering is defined through25 Pαβ =

 N  1  3ˆ ˆ 1 diα diβ − δαβ , N i=1 2 2

(25)

where α, β = x, y, z, dˆ i is the particles’ deformation unit vector, and δ αβ is the Kronecker delta function. Pαβ is symmetric and traceless, while diagonalization of Pαβ results in three different eigenvalues, P2 , P1 , and P0 . We use the largest, P2 , to define a system order parameter which is 0 for a completely disordered phase and 1 for a completely ordered phase. The average P2  as a function of the applied shear rate is shown in Fig. 6. From γ˙ = 0.1 to 1 s−1 P2  remains close to zero indicating that the system is effectively disordered. From γ˙ = 1 s−1 P2 increases and reaches 0.6 for γ˙ = 100 s−1 , demonstrating increased ordering. To calculate the preferred direction of the particles we take the eigenvector, D, associated with P2 and project it on the x-direction, D · eˆ x = cos θ,

(26)

where θ is the angle of D in respect to the x-direction. For γ˙ < 1 s−1 , θ  is close to the equilibrium value of 45◦ . With increasing γ˙ , the angle decreases to approximately 10◦ at

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 128.82.252.58 On: Sat, 20 Dec 2014 01:19:10

Santos de Oliveira et al.

104903-9

J. Chem. Phys. 140, 104903 (2014) 3

1.0

10

0.8 < θ > / 45

. η (γ) [Pa.s]

< P2 > [ - ] 0.6

o

0.4

2

10

Experiment Spherical particles Deformable particles

0.2 0.0 0.1

1

. -1 γ [s ]

10

1

100

10 -1 10

FIG. 6. Average order parameter P2  and average orientation angle θ  scaled by 45o of the deformable particle along the x-direction for a range of shear rate values.

0

10

1

. -1 γ [s ]

10

2

10

FIG. 8. Comparison between the simulated and experimental shear viscosity for a HPC solution. Results for both the spherical and deformable particle coarse graining are shown. The solid green line is a fit of the simulation data using the Cross model (Eq. (24)) for the deformable particles approach.

γ˙ = 100 s−1 and hence demonstrates shear-induced alignment of deformable particles along the flow direction. B. HPC aqueous solution with deformable particles

As a second example we have used the deformable particle model to reproduce the linear and nonlinear rheology of an aqueous solution of HPC.21, 22 The calculated G(t) and storage and loss moduli are presented in Fig. 7 where we demonstrate good agreement between experiments and simulations. We also show the shear relaxation modulus separately for each force contributing to G(t) in Fig. 7(a). The as4

10

(a) 3

G (t)

10

Gexp(t) G(t) Gc(t) Gt(t) Gf(t)

2

10

1

10

-4

10

-3

-2

10

10

t [s]

-1

10

0

10

4

10

(b)

sociated simulation parameters for the HPC fluid are listed in Table I. As for the PIB polymer solution, we find that the Flory-Huggins conservative potential is mostly responsible for the decay of the stress correlation at higher times. However unlike the PIB polymer solution, the FENE contribution is dominant at lower times, while the transient contribution describes intermediate times. As for the PIB solution, the sum of the relaxation curves from each potential effectively yields the full G(t) curve demonstrating that cross-correlations between different force types do not significantly contribute to G(t). As for the PIB solution, we computed the shear viscosity for varying shear rates and compared the flow curve to experimental measurements as shown in Fig. 8 where we observe good agreement between simulations and experiments. Unlike for the PIB polymer solution, this nonlinear rheology is a prediction from the algorithm since the parameters were tuned using the linear rheology only. Also in this case, the spherical particle model was not able to reproduce the correct shear flow curves at higher shear rates. The weak shearthinning effect for the HPC solution could only be captured by the simulations using the FENE deformable particles. For the spherical particle model a higher slope, or shear viscosity, at shear rates beyond 3 s−1 is observed. The difference, as for the PIB solution, can be attributed to the ability of the deformable particles to elongate and hence resist shear deformation.

3

[Pa]

10

2

VI. DISCUSSION AND CONCLUSIONS

G (ω)

10

1

10

G (ω)

0

10 -1 10

0

10

1

ω [rad/s]

10

2

10

FIG. 7. (a) G(t) for a HPC polymer solution. Also shown are the contributions to G(t) computed for conservative, transient or FENE stresses. (b) Comparison of the loss and storage moduli from experiments (solid line) and simulations (symbols).

We have simulated shear-thinning viscoelastic fluids with an extended Responsive Particle Dynamics method, which is a coarse-grained mesoscale model that can recover both the linear and nonlinear rheology of viscoelastic fluids. We have extended the RaPiD representation of coarse-grained polymer chains such that each polymer is depicted as a deformable particle. For these deformable particles, the conservative potential used in the spherical particle representation is adjusted by including a FENE potential which imposes a free energy penalty upon deformation. The description of entanglements between particles is updated, thus prolonging the effect of the

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 128.82.252.58 On: Sat, 20 Dec 2014 01:19:10

104903-10

Santos de Oliveira et al.

J. Chem. Phys. 140, 104903 (2014)

ACKNOWLEDGMENTS

4

. N1 (γ) [Pa]

10

The authors acknowledge financial support from the EU through the FP7 European Soft Matter Infrastructure, ESMI Grant Agreement No. 262348. I. S. Santos de Oliveira acknowledges R. Pasquino for providing experimental data for the HPC solution.

Experiment Spherical particles Deformable particles

3

10

2

10

1

10

0

10

1

10

. -1 γ [s ]

2

10

FIG. 9. Comparison of the first normal stress difference N1 from experiments and simulations for the PIB polymer solution. Results for both the spherical and deformable particle coarse graining are shown.

transient forces between the particles during shear flow. This strategy mimics the stretching of the chains when subjected to shear flow, an observation expected for real fluids. Notice that large shear flows tend to lower entangling of neighbouring polymers but that at the same time forces between stretched polymers are transmitted over larger distances. It is the absence of the latter effect which causes too strong shear thinning in the original, spherical RaPiD model. We have applied the deformable particle model to simulate two different polymer solutions, and for both cases have recovered to good agreement the experimental moduli and shear-thinning behaviours. We also considered the contributions from each force to the G(t) curves for both solutions. In both cases, the Flory-Huggins conservative potential describes the curve for longer times. However, at low to intermediate times the shear relaxation modulus is mostly accounted for by the transient force for PIB polymer solution but by the FENE contribution for the HPC solution. Through comparison with experimental data we calculated the shear thinning exponent for both solutions confirming the intuition that it depends on the maximum elongation d0 more so than the spring constant k of the FENE. Since the changes made to the original representation of the individual particles addresses their elongational properties, it is interesting to check how these changes influence the predicted first normal stress difference. In Fig. 9, we have plotted the experimental first normal stress differences for the PIB polymer solution together with the ones predicted with the spherical and the deformable model. Similar to the shear viscosities, with the spherical model the thinning behaviour is too pronounced, whereas the deformable model does a much better job. Unfortunately, the absolute values are still not perfect. In conclusion, we have shown that representing RaPiD particles as deformable particles can provide a better control over the simulated nonlinear rheology for viscoelastic polymer solutions. The observations obtained here can be used as a reference to model other complex fluids such as polydisperse fluids and long chain polymer melts using the RaPiD simulation method.

APPENDIX: COMPARISON BETWEEN RAPID SPHERICAL PARTICLES AND DEFORMABLE PARTICLES

In the main text, we have seen that the configuration of each polymer may be described either by the center of mass position vector and the deformation vector, (ri , di ), or by the two position vectors of the entanglement centers, (ri,1 , ri,2 ). The latter are related to the first by 1 ri,m = ri + (−1)m di . 2

(A1)

Using the equation of motion given in Sec. II C, we find for the displacements of the entanglement centers   ∂ 1 1 + d dt dri,m = − ξi 4ξi ∂ri,m   ∂ 1 1 − d dt − ξi 4ξi ∂ri,m

2 kB T dt 2 kB T dt d + i + i . (A2) ξi 4ξid Here m is the complement of m, i.e., m = 2 − 12 (1 + (−1)m ). If we now assume that hydrodynamic interactions should not be important, the second term on the right hand side vanishes, so ξid = ξi /4. As a result

dri,m

2 ∂ =− dt + i ξi ∂ri,m

4 kB T dt . ξi

(A3)

√ Here we have used i + di = 2i . Next, let us compare the spherical particle model with the deformable particle model. The transient force on particle i, due to the presence of particle j, in the spherical particle model and the deformable model read Fsij = α s (nij − n0 (rij ))

dn0 ri − rj drij rij

(A4)

and   1 1 dn0 rim − rj n nim,j n − n0 (rim,j n ) 4 4 drim,j n rim,j n m,n   1  = nim,j n − n0 (rim,j n ) 4 m,n

Ftij = α

×

dn0 rim − rj n , drim,j n rim,j n

(A5)

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 128.82.252.58 On: Sat, 20 Dec 2014 01:19:10

104903-11

Santos de Oliveira et al.

J. Chem. Phys. 140, 104903 (2014)

respectively. Taking the limit of the latter for d going to 0 we find    dn0 ri − rj 1 t nim,j n − n0 (rij ) . (A6) Fij = α 4 drij rij m,n If we now assume that in equilibrium mostly d = 0, and we set  m, n nim, jn = nij we find that the spherical particle and deformable particle models yield the same transient forces in the case α = 4α s . Under these assumptions all other forces are the same in both models as well, since for d = 0 the FENE force is zero and the Flory-Huggins forces are independent of d by construction. In order for the dynamics in both models to be the same we put the frictions in both models equal to each other    1 |nim,j n | n0 (rim,j n ) = ξes |nij |n0 (rij ). (A7) ξe 4 m,n In the limit of d = 0 and taking the typical situation when nij = n0 (rij ) and nim,j n = 14 n0 (rij ), we obtain ξe = ξes .

(A8)

In the same way, we must make the dynamics of the entanglements equal in both models. Consider, in the limit when d = 0,     1 dnim,j n = nim,j n dt n0 (rij ) − τ m,n m,n +

 m,n

 im,j n

2 kB T dt . ατ

(A9)

Since this must be the equation for the change of nij , we find √ τ = τ s . Notice that with m,n im,j n = 4 ij and α = 4α s

the fluctuation-dissipation theorem is also recovered for the spherical particle model. 1 G.

Zhang, K. Daoulas, and K. Kremer, Macromol. Chem. Phys. 214, 214 (2013). 2 L. Liu, J. Padding, W. den Otter, and W. Briels, J. Chem. Phys. 138, 244912 (2013). 3 A. Clark, J. McCarty, and M. Guenza, J. Chem. Phys. 139, 124906 (2013). 4 D. M. Sussman and K. S. Schweizer, Phys. Rev. Lett. 109, 168306 (2012). 5 W. J. Briels, Soft Matter 5, 4401 (2009). 6 A. van den Noort, W. K. den Otter, and W. J. Briels, Europhys. Lett. 80, 28003 (2007). 7 P. Kindt and W. J. Briels, J. Chem. Phys. 127, 134901 (2007). 8 J. T. Padding, E. van Ruymbeke, D. Vlassopoulos, and W. J. Briels, Rheol. Acta 49, 473 (2010). 9 I. S. Santos de Oliveira, A. van den Noort, J. T. Padding, W. K. den Otter, and W. J. Briels, J. Chem. Phys. 135, 104902 (2011). 10 F. Snijkers, G. D’Avino, P. L. Maffetone, F. Greco, M. Hulsen, and J. Vermant, J. Rheol. 53, 459 (2009). 11 C. N. Likos, Phys. Rep. 348, 267 (2001). 12 P. J. Flory, J. Chem. Phys. 10, 51 (1942). 13 M. L. Huggins, J. Chem. Phys. 9, 440 (1941). 14 I. Pagonabarraga and D. Frenkel, J. Chem. Phys. 115, 5015 (2001). 15 S. Y. Trofimov, E. L. F. Nies, and M. A. J. Michels, J. Chem. Phys. 117, 9383 (2002). 16 J. T. Padding, L. V. Mohite, D. Auhl, W. J. Briels, and C. Bailly, Soft Matter 7, 5036 (2011). 17 J. T. Padding, L. V. Mohite, D. Auhl, T. Schweizer, W. J. Briels, and C. Bailly, Soft Matter 8, 7967 (2012). 18 M. Doi and S. F. Edwards, The Theory of Polymers Dynamics (Oxford Science Publications, Oxford, England, 1986). 19 R. L. C. Akkermans and W. J. Briels, J. Chem. Phys. 115, 6210 (2001). 20 M. P. Allen and D. J. Tildesley, Computer Simulations of Liquids (Oxford University Press, Oxford, England, 1987). 21 R. Pasquino, F. Snijkers, N. Grizzutti, and J. Vermant, Rheol. Acta 49, 993 (2010). 22 R. Pasquino, D. Panariello, and N. Grizzuti, J. Colloid Interface Sci. 394, 49 (2013). 23 M. M. Cross, J. Colloid Sci. 20, 417 (1965). 24 C. Baig, V. Mavrantzas, and M. Kröger, Macromolecules 43, 6886 (2010). 25 P. G. Gennes and J. Prost, The Physics of Liquid Crystals (Oxford University Press, Oxford, England, 1993).

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 128.82.252.58 On: Sat, 20 Dec 2014 01:19:10